Ramesh Kolluru et al

*Ramesh Kolluru

Post Doctoral Fellow, Department of Aerospace Engineering, Indian Institute of
Science, Bangalore, Karnataka, India

# Robust and accurate central algorithms for Multi-Component mixture equations
with Stiffened gas EOS

Ramesh Kolluru  S V Raghurama Rao  G N Sekhar  Department of Aerospace
Engineering, Indian Institute of Science, Bangalore,Karnataka, India
Department of Aerospace Engineering, Indian Institute of Science,
Bangalore,Karnataka, India  Department of Mathematics, BMS College of
Engineering, Bangalore,Karnataka, India  Department of Mechanical Engineering,
BMS College of Engineering, Bangalore,Karnataka, India<EMAIL_ADDRESS>Ramesh Kolluru  S V Raghurama Rao  G N Sekhar

(26 April 2016; 6 June 2016; 6 June 2016)

###### Abstract

[Summary]Simple and robust algorithms are developed for compressible Euler
equations with stiffened gas equation of state (EOS), representing gaseous
mixtures in thermal equilibrium and without chemical reactions. These
algorithms use fully conservative approach in finite volume frame work for
approximating the governing equations. Also these algorithms used central
schemes with controlled numerical diffusion for this purpose. Both Mass
fraction ($Y$) and $\gamma$ based models are used with RICCA and MOVERS+
algorithms to resolve the basic features of the flow fields. These numerical
schemes are tested thoroughly for pressure oscillations and preservation of
the positivity of mass fraction at least in the first order numerical methods.
Several test cases in both 1D and 2D are presented to demonstrate the
robustness and accuracy of the numerical schemes.

, , and , (2019), Simple,Accurate and Robust Algorithms for Multi-Component
mixture equations with Stiffened gas EOS, , .

###### keywords:

MOVERS, MOVERS+, RICCA, Contact-discontinuity,$\gamma-$ based approach

††articletype: Article Type00footnotetext: Abbreviations: MOVERS, Method of
Optimum Viscosity for Enhanced Resolution of Shocks; RICCA Riemann Invariants
based Contact-discontinuity Capturing Algorithm;

##  1 Introduction

Atmospheric air is a mixture of gases which are compressible in nature. Each
of the components in the mixture have different physical and thermodynamical
properties, and very often in modelling the flow, air is assumed to be a
single component gas with constant properties. There are many applications
where due consideration should be given to each of the components in the
mixture such as gasoline and air mixture entering the combustion chamber and
combustion products exhausting from the engines. There are instances where
liquids and gases exist together like bubbles moving in the liquid, spray of
paint facilitated through a nozzle. In all the situations mentioned above
fluids exist as a mixture or as different components separated by interfaces.

Many times, the contribution of the individual components are negligible or
the variation in the properties of the components do not contribute
significantly to the flow field and hence they can be neglected. If the
properties of the components vary at large, then the individual effect of the
components are resolved or their combined effect on the mixture has to be
studied. In these cases, the classic model of single component compressible
fluid may not be appropriate.

Broadly the flow of these fluids can be classified into two categories: a)
pure interface problems, and b) multicomponent flows. For pure interface
problems, the thermodynamic properties of the fluids change only across the
interface whereas in multicomponent flows the properties vary throughout the
flow field. In pure interface problems, apart from solving for the dynamics of
each component, the interface is also tracked by a specific method like level
set method. In multicomponent flows the modelling is done without tracking any
interface. Nature of multicomponent flows can vary from low subsonic flows to
hypersonic reacting flows. Low subsonic flows often are coupled with
combustion related phenomena and therefore are not easily amenable to
numerical modelling. Modelling of supersonic and hypersonic flows can take
advantage of the sophisticated numerical methods developed for hyperbolic
systems in the past few decades, though treating supersonic and hypersonic
combustion problems are non-trivial. Some of the important contributions in
modelling multicomponent flows are due to 15, 16, 5, 8, 3, 4, 38, 6, 18, 19,
34, 35, 36.

Fernandez et al.12, aimed at constructing an efficient conservative numerical
scheme for computation of multi-species flows. The governing equations are
Euler equations and additional equations for the species with different
molecular weights and specific heats are considered. Approximate Riemann
solver of Roe has been used and modifications for evaluation of $\gamma$ in
the Roe matrix have been suggested. Donor cell approximation method for
species equations is modelled and compared with the modifications of the Roe
scheme. They conclude that the modified Roe scheme performs better than the
donor cell approximation method.

Larrouturou et al.18, have reviewed various numerical methods for
multicomponent perfect and real gas models. They have suggested modifications
for Osher, Steger-Warming, van Leer and Roe schemes for the application to
multicomponent perfect and real gases. They show clearly that, for
multicomponent flows Roe’s conditions (consistency, conservation and
hyperbolicity) get satisfied only when $\gamma(U)$ is constant and hence the
extension of Roe scheme for multicomponent mixture flows is not an easy task.

Karni 15 has carried out modelling of multicomponent fluids using Euler
equations with an additional equation for the species. Both conservative and
non-conservative form of the equations are considered and primitive form of
the equations are recommended to avoid the pressure oscillations occurring
near the material interface. Four different models of the governing equations
with variable $\gamma$ and a level set method based on distance function are
used in both conservative and primitive form. Roe linearisation method is used
in numerical simulation and compared with second order upwind methods. It is
concluded that any fully coupled conservative based numerical scheme leads to
pressure oscillations and non-preserving of positivity of mass fractions. She
recommended the use of primitive variable based approach in order to avoid the
pressure oscillations. The use of non-conservative form, however, leads to
conservative errors and incorrect shock positions.

Abgrall 3, has used a quasi-conservative approach for the calculations for
multicomponent cases and proved that an additional evolution equation of
$\gamma$ and in particular of the form $\frac{1}{\left(\gamma-1\right)}$, is
suggested to preserve mass fraction positivity and to avoid pressure
oscillations.

Abgrall and Karni 5 have reviewed numerical algorithms commonly used in the
simulations of multicomponent compressible fluid flow. They conclude that if
separate equations for individual species are solved along with the mixture
equations, then the numerical scheme developed preserves pressure equilibrium
and mass fraction positivity.

Keh-Ming Shyue in 34, 35, 36, 37, has utilized Abgrall’s model 3 for
compressible multicomponent flow problems using stiffened gas EOS, van der
Waals EOS, Mie-Gruneisen EOS, Tait EOS.

Overall the basic issues in extending the single fluid conservative numerical
schemes to multicomponent flows are

  1. 1.

preserving the positivity of mass fraction,

  2. 2.

avoiding pressure oscillations even in the first order numerical scheme,

  3. 3.

difficulties in extension to more than 2 components.

In this work novel and accurate central solvers MOVERS-n, MOVERS-1, developed
by 13 along with MOVERS+, RICCA as explained in 28 are applied to
multicomponent flows to address some of the above issues.

##  2 Governing equations for mixture with two components or species

From the literature it is observed that there are many different ways in which
the governing equations can be formulated. A simple case of non-reacting
mixture equations with two components and without diffusion is considered in
the present work. Two models based on mass fraction Y and $\gamma$ are chosen
to test the algorithms in conservative cell centered finite volume frame work.
In the following sections the governing equations for these models and basic
algorithms used to discretise them are discussed briefly.

###  2.1 Mass fraction based model

Consider the mixture of gasses consisting of two species with following
mixture properties: pressure p, density $\rho$, velocity u and temperature
$T$. The mixture pressure is given by Dalton’s Law $p=p_{1}+p_{2}$, mixture
density $\rho=\rho_{1}+\rho_{2}$, the mass fraction of the species
$Y_{k}=\frac{\rho_{k}}{\rho},\leavevmode\nobreak\ k=1,2$. Specific heat at
constant pressure and constant volume of individual species,
$c_{p_{k}},c_{v_{k}},\leavevmode\nobreak\ k=1,2$ are considered to be
constant, and the ratio of specific heats of individual species is given by
$\gamma_{k}=\frac{{c_{p}}_{k}}{{c_{v}}_{k}},\leavevmode\nobreak\ k=1,2$. The
governing equations for the mixture in conservation form are given by (2.1)

$\displaystyle\left.\begin{aligned} \frac{\partial\rho}{\partial
t}+\frac{\partial\left(\rho u\right)}{\partial x}&=0,\\\ \frac{\partial\rho
u}{\partial t}+\frac{\partial\left(\rho u^{2}+p\right)}{\partial x}&=0,\\\
\frac{\partial\rho E_{t}}{\partial t}+\frac{\partial\left[\left(\rho
E_{t}+p\right)u\right]}{\partial x}&=0,\\\ \frac{\partial\left(\rho
Y_{k}\right)}{\partial t}+\frac{\partial\left(\rho Y_{k}u\right)}{\partial
x}&=0,\qquad k=1,2\\\ Y_{1}+Y_{2}&=1\end{aligned}\right\\}$ (2.1)

The value of ratio of specific heats, $\gamma$, for the mixture, is defined as
$\gamma=\frac{{c_{p}}_{mixture}}{{c_{v}}_{mixture}}=\frac{\sum_{k}Y_{k}\gamma_{k}{c_{v}}_{k}}{\sum_{k}Y_{k}{c_{v}}_{k}}$
and the mixture pressure is given by $p=(\gamma-1)(\rho E_{t}-\frac{\rho
u^{2}}{2})$. The equation of state for each individual component can be
described by a function $p=p(\rho,e)$. These governing equations are
represented in compact notation as in the first equation of (2.2) where
enthalpy of the mixture is given by $H=E_{t}+\frac{p}{\rho}$. The ratio of
specific heats for the mixture, $\gamma$ is a function of the conserved
variable vector $U$, as
$\gamma=\frac{U_{4}\gamma_{1}{c_{v}}_{1}+(U_{1}-U_{4})\gamma_{2}{c_{v}}_{2}}{U_{4}{c_{v}}_{1}+(U_{1}-U_{4}){c_{v}}_{2}}=\gamma\left(U_{1},U_{4}\right)$.
This property of $\gamma$ for the mixture plays a role in determining the
hyperbolicity of the governing equations. It can also be observed that (2.2)
is extension of Euler equations with an additional equations for the mass
fraction of individual component gases. Hence if this set of governing
equations satisfies the hyperbolicity principle then all the algorithms which
are designed for Euler equations can be in principle extended to
multicomponent fluids.

$\displaystyle\frac{\partial U}{\partial t}+\frac{\partial
F\left(U\right)}{\partial x}$ $\displaystyle=0$ (2.2) $\displaystyle
U=\begin{bmatrix}U_{1}\\\ U_{2}\\\ U_{3}\\\
U_{4}\end{bmatrix}=\begin{bmatrix}\rho\\\ \rho u\\\ \rho E\\\ \rho
Y_{k}\end{bmatrix};F(U)=\begin{bmatrix}\rho u\\\ \rho u^{2}+p\\\ (\rho
E+p)u\\\ \rho Y_{k}u\end{bmatrix}$ $\displaystyle=\begin{bmatrix}U_{2}\\\
\frac{\left(3-\gamma\right)}{2}\frac{U_{2}^{2}}{U_{1}}+\left(\gamma-1\right)U_{3}\\\
\frac{\left(3-\gamma\right)}{2}\frac{U_{3}U_{2}}{U_{1}}-\frac{\left(\gamma-1\right)}{2}\frac{U_{2}^{3}}{U_{1}^{2}}\\\
\frac{U_{4}U_{1}}{U_{2}}\end{bmatrix}$ (2.3)

###  2.2 Hyperbolicity and eigenstructure for the mixture model

To demonstrate the hyperbolicity of equations (2.2), it is required to
evaluate the flux Jacobian matrix, it’s eigenvalues and corresponding
eigenvectors. The flux Jacobian matrix of the governing equations is given by
(2.4)

$\displaystyle A(U)$ $\displaystyle=\begin{bmatrix}0&1&0&0\\\
\frac{\left(\gamma-3\right)}{2}\frac{U_{2}^{2}}{U_{1}^{2}}+B&\left(3-\gamma\right)\frac{U_{2}}{U_{1}}&\left(\gamma-1\right)&B^{\prime}\\\
\frac{\left(\gamma-3\right)}{2}\frac{U_{3}U_{2}}{U_{1}^{2}}+uB+\left(\gamma-1\right)\frac{U_{2}^{3}}{U_{1}^{3}}&\frac{\left(3-\gamma\right)}{2}\frac{U_{3}}{U_{1}}-\frac{3\left(\gamma-1\right)}{2}\frac{U_{2}^{2}}{U_{1}^{2}}&\frac{\left(3-\gamma\right)}{2}\frac{U_{2}}{U_{1}}&uB^{\prime}\\\
-\frac{U_{4}}{U_{2}}&\frac{-U_{4}U_{1}}{U_{2}^{2}}&0&\frac{U_{2}}{U_{1}}\end{bmatrix}$
(2.4)

It can be observed that the flux Jacobian matrix is a function of $\gamma$ and
its derivatives given by (2.5)

$\displaystyle B$
$\displaystyle=\frac{p}{\left(\gamma-1\right)}\frac{\partial\gamma}{\partial
U_{1}},$ (2.5) $\displaystyle B^{\prime}$
$\displaystyle=\frac{p}{\left(\gamma-1\right)}\frac{\partial\gamma}{\partial
U_{4}}.$

The flux Jacobian matrix (2.4) in terms of specific total enthalpy of mixture
$H$ is given by (2.6).

$\displaystyle A(U)=\frac{\partial F\left(U\right)}{\partial
U}=\begin{bmatrix}0&1&0&0\\\
\frac{\left(\gamma-3\right)}{2}u^{2}+B&\left(3-\gamma\right)u&\left(\gamma-1\right)&B^{\prime}\\\
\frac{\left(\gamma-1\right)}{2}u^{3}+Bu-uH&H-\left(\gamma-1\right)u^{2}&\gamma
u&B^{\prime}u\\\ -Y_{1}u&Y_{1}&0&u\end{bmatrix}$ (2.6)

The eigenvalues of the matrix $A(U)$ 8 are $(u+a,u,u,u-a$, where
$a=\sqrt{\frac{\gamma p}{\rho}})$, and the right eigenvectors are

$\displaystyle r_{1}=\begin{bmatrix}1\\\ u-a\\\ H-ua\\\
Y_{1}\end{bmatrix},r_{2}=\begin{bmatrix}1\\\ u\\\
\frac{u^{2}}{2}-\frac{B}{\left(\gamma-1\right)}\\\
0\end{bmatrix},r_{3}=\begin{bmatrix}0\\\ 0\\\
-\frac{B^{\prime}}{\left(\gamma-1\right)}\\\
1\end{bmatrix},r_{4}=\begin{bmatrix}1\\\ u+a\\\ H+ua\\\ Y_{1}\end{bmatrix}$
(2.7)

The above system is hyperbolic as eigenvalues are real and the eigenvectors
are linearly independent.

###  2.3 Multicomponent fluid simulations with $\gamma$-based model and
stiffened gas EOS

In his fundamental work, Abgrall 3 has quoted that any numerical scheme
designed for compressible Euler equations extended to multicomponent flows
would generate pressure oscillations. It has also been suggested by the author
that use of $\frac{1}{(\gamma-1)}$ as the parameter in the quasi-conservative
approach would eliminate the pressure oscillations, this has been demonstrated
by Shyue 34, 35, 36, 37. For $\gamma$ based model the mixture equations in 1D
are given by (2.8).

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\left(\partial\rho
u\right)}{\partial x}$ $\displaystyle=0$ (2.8)
$\displaystyle\frac{\partial\rho u}{\partial t}+\frac{\partial\left(\rho
u^{2}+p\right)}{\partial x}$ $\displaystyle=0$
$\displaystyle\frac{\partial\rho E}{\partial t}+\frac{\partial\left[\left(\rho
E+p\right)u\right]}{\partial x}$ $\displaystyle=0$

Here, $Y^{i}$ represent the volume fraction of the component gases in a given
cell or control volume and $\rho,u,p,p_{\infty},\gamma$ represent the mixture
density, mixture velocity, mixture pressure, mixture stiffened pressure and
mixture gamma respectively.

$\displaystyle\frac{\left(\partial\frac{\rho}{\left(\gamma-1\right)}\right)}{\partial
t}+\frac{\left(\partial\frac{\rho u}{\left(\gamma-1\right)}\right)}{\partial
x}=0$ (2.9) $\displaystyle\frac{\left(\partial\frac{\rho\gamma
p_{\infty}}{\left(\gamma-1\right)}\right)}{\partial
t}+\frac{\left(\partial\frac{\rho\gamma
p_{\infty}u}{\left(\gamma-1\right)}\right)}{\partial x}=0$ (2.10)
$\displaystyle\frac{p+\gamma p_{\infty}}{\gamma-1}=\rho e$ (2.11)

In the above equations (2.11) refer to stiffened gas EOS. This EOS would
revert to perfect gas EOS in the limit $p_{\infty}\rightarrow 0$. These
equations (2.9) along with (2.8) are referred to as $\gamma$-based model by
Abgrall 2.

###  2.4 Numerical Methodology and Novel Algorithms

The governing equations are Euler equations in conservative form as given in
(2.12a)

$\displaystyle\frac{d\overline{U}}{dt}=-R,R$
$\displaystyle=\frac{1}{\Omega}\left[\sum_{i=1}^{N}{F_{c}\cdot\hat{n}\leavevmode\nobreak\
dS}\right],$ (2.12a) $\displaystyle\overline{U}$
$\displaystyle=\frac{1}{\Omega}\int_{\Omega}Ud\Omega\leavevmode\nobreak\ .$
(2.12b)

where $U$ is conserved variable vector, $F_{c}$ is convective flux vector on
an interface, $R$ representing net flux from a given control volume, $\Omega$
volume of control volume and $N$ representing number of control surfaces for a
given control volume.

(j-1)(j)(j+1)$F_{L}=F(U_{L})$$F_{R}=F(U_{R})$$j-\frac{1}{2}$$j+\frac{1}{2}$$F_{\mathrm{I}}$
Figure 1: Typical finite volume in 1D

The convective flux on any interface of a control volume as shown in figure(1)
for any stable scheme can be written as the sum of an average flux across the
interface and a numerical dissipative flux as given in (2.13a).

$\displaystyle
F_{\mathrm{I}}=\frac{1}{2}\left[F_{L}+F_{R}\right]-d_{\mathrm{I}}$ (2.13a)
$\displaystyle
d_{\mathrm{I}}=\frac{\mid\alpha_{\mathrm{I}}\mid}{2}\left(U_{R}-U_{L}\right)$
(2.13b)

where $\alpha_{\mathrm{I}}$ coefficient of numerical diffusion and
$d_{\mathrm{I}}$ represents the numerical dissipative flux. In the present
work the coefficient of numerical diffusion is determined by MOVERS 13, RICCA
and MOVERS+ 28 as briefed in sections(2.4.1,2.4.2 $\&$ 2.4.3) are utilised to
simulate the mixture equations based on mass fraction model with perfect gas
EOS and $\gamma$ based model with stiffened gase EOS. These algorithms are
independent of eigenstructure of the underlying hyperbolic system and can be
extended to any arbitrary EOS.

####  2.4.1 MOVERS

The central scheme of interest is due to 13 who introduced a new central
scheme named MOVERS (_Method of Optimal Viscosity for Enhanced Resolution of
Shocks_) which can capture grid aligned shocks and contact-discontinuities
accurately. This coefficient of numerical diffusion rewritten in terms of RH
conditions is given by (2.14)

$\displaystyle\lvert\alpha_{I}\rvert_{i}=\left|s_{i}\right|=\left|\frac{\Delta
F_{i}}{\Delta U_{i}}\right|,\quad
i=1,2,3,\quad\Delta(\cdot)=(\cdot)_{R}-(\cdot)_{L}$ (2.14)

In order to introduce boundedness and further stabilize the numerical scheme,
$\alpha_{I}$ is to be restricted to a physically feasible range of eigenvalues
of the flux Jacobian matrix. This process known as wave speed correction
(2.15) is incorporated such that the coefficient of numerical diffusion lies
within the eigenspecturm of the flux Jacobian i.e.,
$\alpha_{I}\in\left[\lambda_{max},\lambda_{min}\right]$.

$\displaystyle\lvert\alpha_{\mathrm{I}}\rvert=\begin{cases}\lambda_{\textit{max}},\leavevmode\nobreak\
\leavevmode\nobreak\ if\leavevmode\nobreak\ \leavevmode\nobreak\
\lvert\alpha_{\mathrm{I}}\rvert>\lambda_{\textit{max}}\\\
\lambda_{\textit{min}},\leavevmode\nobreak\ \leavevmode\nobreak\
if\leavevmode\nobreak\ \leavevmode\nobreak\
\lvert\alpha_{\mathrm{I}}\rvert<\lambda_{\textit{min}}\\\
\lvert\alpha_{\mathrm{I}}\rvert,\\\ \end{cases}$ (2.15)

This method is independent of eigenstructure of the underlying hyperbolic
systems, is simple and can capture grid-aligned stationary discontinuities
exactly. Authors 13 introduced two variations of MOVERS: $(i)$ an $n$-wave
based coefficient of numerical diffusion, corresponding to $n$ number of
conservation laws (MOVERS-n) and $(ii)$ a scalar diffusion, corresponding to
the energy equation, referred to as MOVERS-1. The robustness of the basic
scheme has been improvised through its variants by Maruthi N.H. 26 and
extended them to hyperbolic systems for magnetohydrodynamics and shallow water
flows. In this work this algorithm is chosen as the foundation to devise two
new efficient algorithms for hyperbolic systems. The simplicity and accuracy
of this algorithm make this scheme a well-suited base-line solver for further
research, apart from its independency of the eigenstrucure.

####  2.4.2 Riemann Invariant based Contact-discontinuity Capturing Algorithm
(RICCA)

The numerical diffusion evaluated using Riemann Invariant based Contact-
discontinuity Capturing Algorithm (RICCA) is given by

${\alpha}_{\emph{$I$}}=\begin{cases}\qquad\qquad\frac{|V_{nL}|+|V_{nR}|}{2},\qquad\qquad\qquad\quad\text{if
}|\Delta\mathbf{F}|<\delta\ \text{and}\ |\Delta\mathbf{U}|<\delta\\\
\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\\\
max(|V_{nL}|,|V_{nR}|)+sign(|\Delta
p_{\raisebox{-2.0pt}{\scriptsize{$\mathrm{I}$}}}|)a_{\mathrm{I}},\qquad\text{otherwise}\\\
\qquad\qquad\qquad\qquad\quad\end{cases}$ (2.16)

where $a_{\mathrm{I}}=\sqrt{\frac{\gamma
p_{\raisebox{-2.0pt}{\scriptsize\emph{$I$}}}}{\rho_{\raisebox{-2.0pt}{\scriptsize\emph{$I$}}}}}$
is the speed of sound evaluated with the values at the interface given by

$\displaystyle\rho_{\raisebox{-2.0pt}{\scriptsize\emph{$I$}}}=\frac{\rho_{L}+\rho_{R}}{2},p_{\raisebox{-2.0pt}{\scriptsize\emph{$I$}}}=\frac{p_{L}+p_{R}}{2},\Delta{p}_{\raisebox{-2.0pt}{\scriptsize\emph{$I$}}}=(p_{R}-p_{L}).$
(2.17)

####  2.4.3 MOVERS without wave speed correction - _MOVERS+_

The coefficient of numerical diffusion for MOVERS+ is given by

$\displaystyle\lvert d_{\mathrm{I}}\rvert_{j}=\Phi Sign(\Delta
U_{j})\lvert\Delta F_{j}\rvert+\left(\frac{|V_{nL}|+|V_{nR}|}{2}\right)\Delta
U_{j},\quad j=1,2,3$ (2.18)

These two new algorithms RICCA and MOVERS+

  * •

can capture steady contact-discontinuities exactly,

  * •

has sufficient numerical diffusion near shocks so as to avoid shock
instabilities, and

  * •

does not need entropy fix for at sonic points.

A similar strategy was introduced by N.Venkata Raghavendra in 42, 43 to design
an accurate contact-discontinuity capturing discrete velocity Boltzmann scheme
for inviscid compressible flows.

###  2.5 Modifications of upwind methods for multicomponent flows

As mentioned before, the application of upwind methods to multicomponent flows
is non-trivial because these methods are strongly dependent on the
eigenstructure. Larrouturou and Fezoui 18 have reviewed these modifications
needed for upwind methods, which are briefly presented here.

####  2.5.1 Extension of Steger-Warming FVS method to multicomponent gases

The Flux Vector Splitting (FVS) method of Steger-Warming method as given in
TORO_1, leads to the following split flux vectors for Euler equations.

$\displaystyle
F^{\pm}=\frac{\rho}{2\gamma}\begin{bmatrix}\left(u-a\right)\lambda_{1}^{\pm}+2\left(\gamma-1\right)\lambda_{2}^{\pm}+\lambda_{3}^{\pm}\\\
\lambda_{1}^{\pm}+2\left(\gamma-1\right)u\lambda_{2}^{\pm}+\left(u+a\right)\lambda_{3}^{\pm}\\\
\left(H-ua\right)\lambda_{1}^{\pm}+2\left(\gamma-1\right)u^{2}\lambda_{2}^{\pm}+\left(H+ua\right)\lambda_{3}^{\pm}\\\
\end{bmatrix}$ (2.19)

For the mixture equations (2.1), an additional fourth component (for the extra
mass fraction term) for the split flux vectors is given by

$\displaystyle F[4]=\rho uY_{1}=F^{\pm}[1]Y_{1}.$ (2.20)

####  2.5.2 Extension of van Leer FVS method to multicomponent gases

The details of the flux vector splitting developed by van Leer are given in
TORO_1 for Euler equations. Extension of van Leer flux splitting to
multicomponent mixture equations as a function of Mach number is given by

$\displaystyle F=F(\rho,a,M,Y)=\begin{bmatrix}\rho aM\\\ \rho
a^{2}\left(M^{2}+1\right)\\\ \rho
a^{3}M\left(\frac{M^{2}}{2}+\frac{1}{\left(\gamma-1\right)}\right)\\\ \rho
aMY_{1}\end{bmatrix}$ (2.21)

The split fluxes given in 18 are

$\displaystyle F^{\pm}=\frac{1}{4}\rho
a\left(1{\pm}M\right)^{2}\begin{bmatrix}1\\\
\frac{2a}{\gamma}\left(\frac{\left(\gamma-1\right)}{2}M{\pm}1\right)\\\
\frac{2a^{2}}{\left(\gamma^{2}-1\right)}\left(\frac{\left(\gamma-1\right)}{2}M{\pm}1\right)^{2}\\\
Y_{1}\end{bmatrix}$ (2.22)

####  2.5.3 Extension of Roe’s FDS method

Roe’s Flux Difference Splitting (FDS) method, which is an approximate Riemann
solver, cannot be directly extended to multicomponent flows in a trivial way.
In order to evaluate Roe’s numerical flux, the following information is
necessary

  * •

wave strengths $\tilde{\alpha_{i}}$,

  * •

eigenvalues of the flux Jacobian matrix $\tilde{\lambda_{i}}$,

  * •

right eigenvectors of the flux Jacobian matrix $\tilde{R^{\left(i\right)}}$.

The following basic conditions (also called as $U$ property) are to be
satisfied by Roe scheme

  1. 1.

consistency, $A\left(U_{L},U_{R}\right)=A(U)$ if $U_{L}=U_{R}=U$,

  2. 2.

hyperbolicity i.e., flux Jacobian matrix should have real eigenvalues,

  3. 3.

conservation across discontinuities,
$F\left(U_{R}\right)-F\left(U_{L}\right)=A\left(U_{R}-U_{L}\right)$.

The two component scheme for the interface flux given in 18 is

$F\left(U_{L},U_{R}\right)=\frac{1}{2}\left(F(U_{L})+F(U_{R})\right)+\frac{1}{2}|\tilde{A}|\left(U_{L}-U_{R}\right)$
(2.23)

where $U$ represents the average state between the left and right states. The
averaged state is defined as

$\displaystyle\left.\begin{aligned}
U&=\left[\tilde{\rho},\tilde{\rho}\tilde{u},\tilde{\rho}\tilde{E},\tilde{\rho}\tilde{Y}\right]^{T}\\\
\tilde{\rho}&=\frac{\rho_{L}\sqrt{\rho_{L}}+\rho_{R}\sqrt{\rho_{R}}}{\sqrt{\rho_{L}}+\sqrt{\rho_{R}}}\\\
\tilde{u}&=\frac{u_{L}\sqrt{\rho_{L}}+u_{R}\sqrt{\rho_{R}}}{\sqrt{\rho_{L}}+\sqrt{\rho_{R}}}\\\
\tilde{H}&=\frac{H_{L}\sqrt{\rho_{L}}+H_{R}\sqrt{\rho_{R}}}{\sqrt{\rho_{L}}+\sqrt{\rho_{R}}}\\\
\tilde{Y}&=\frac{Y_{L}\sqrt{\rho_{L}}+Y_{R}\sqrt{\rho_{R}}}{\sqrt{\rho_{L}}+\sqrt{\rho_{R}}}\\\
\tilde{A}(U)&=\begin{bmatrix}0&1&0&0\\\
\frac{\left(\tilde{\gamma}-3\right)}{2}\tilde{u}^{2}+\tilde{B}&\left(3-\tilde{\gamma}\right)\tilde{u}&\left(\tilde{\gamma}-1\right)&B^{\prime}\\\
\frac{\left(\tilde{\gamma}-1\right)}{2}\tilde{u}^{3}+\tilde{B}\tilde{u}-\tilde{u}\tilde{H}&\tilde{H}-\left(\tilde{\gamma}-1\right)\tilde{u}^{2}&\tilde{\gamma}\tilde{u}&\tilde{B}^{\prime}\tilde{u}\\\
-\tilde{Y_{1}}\tilde{u}&\tilde{Y_{1}}&0&\tilde{u}\end{bmatrix}\end{aligned}\right\\}$
(2.24)

The matrix $\tilde{A}$ as given in (2.24) is diagonalisable, and its
eigenvalues, as given in 8, are
$\left(\tilde{u}-\tilde{a},\tilde{u},\tilde{u},\tilde{u}+\tilde{a}\right)$
where
$\left(\tilde{a}^{2}=\left(\tilde{\gamma}-1\right)\left(\tilde{H}-\frac{\tilde{u}^{2}}{2}\right)\right)$.
The following remarks are given in 8, 6, which highlights the conditions under
which Roe scheme is not applicable.

######  Remark 2.1.

For two component fluid flow, the conservation property is satisfied only if
$\leavevmode\nobreak\ \gamma_{1}\leavevmode\nobreak\
=\gamma_{2}\leavevmode\nobreak\ =\textrm{Constant}$. If
$\gamma_{1}\neq\gamma_{2}$, then pressure oscillations are observed in the
case of steady contact discontinuities.

######  Remark 2.2.

In order to satisfy the conservation property $\tilde{A}\left(U\right)$ has to
be modified to $A\left(\tilde{U}\right)$, for which the following definitions
of
$\tilde{B}^{\prime}=\frac{c_{v1}c_{v2}\left(\gamma_{1}-\gamma_{2}\right)\tilde{T}}{\tilde{Y_{1}}c_{v1}+\left(1-Y_{1}\right)c_{v2}}$
and
$\tilde{T}=\frac{T_{L}\sqrt{\rho_{L}}+T_{R}\sqrt{\rho_{R}}}{\sqrt{\rho_{L}}+\sqrt{\rho_{R}}}\neq
T\left(\tilde{U}\right)$ are used. The Jacobian matrix $\tilde{A}$ is
diagonalisable only if this modification is incorporated.

######  Remark 2.3.

This expression
$\tilde{B}^{\prime}=\frac{c_{v1}c_{v2}\left(\gamma_{1}-\gamma_{2}\right)\tilde{T}}{\tilde{Y_{1}}c_{v1}+\left(1-Y_{1}\right)c_{v2}}$
is not easily extendable to a mixture of more than 2 components.

######  Remark 2.4.

If the two states $U_{L}$ and $U_{R}$ are supersonic and satisfy $u_{L}\geq
a_{L}$, $u_{R}\geq a_{R}$, and when using upwind schemes, the flux
$\phi\left(U_{L},U_{R}\right)=F\left(U_{L}\right)$, is satisfied in all upwind
schemes, modified for multicomponent cases, except for Roe scheme.

######  Remark 2.5.

The construction of Jacobian matrix of average the state
$\tilde{A}\left(U_{L},U_{R}\right)=A\left(\tilde{U}\right)$, is impossible for
acomplex equation of state for multicomponent cases.

It can be observed that direct extension of Roe scheme to multicomponent
fluids is not an easy task. The new central solvers introduced in
section(2.4.1,2.4.2,2.4.3), RICCA and MOVERS+, along with MOVERS-1 and
MOVERS-n, do not require any of the above modifications as they are not
dependent on eigenstructure at all.

##  3 Results for mass fraction based model

In this section the test cases used for the validation of the central solvers,
MOVERS-1, MOVERS-n, MOVERS+ and RICCA, for multicomponent flows are discussed.
Initial conditions for the 1-D shock tube are given in table (2).

Test case | $\rho_{L}$ | $H_{L}$ | $u_{L}$ | $\rho_{R}$ | $H_{R}$ | $u_{R}$ | $\gamma_{L}$ | $\gamma_{R}$ | Time
---|---|---|---|---|---|---|---|---|---
1 | 1.0 | 1.0 | 0.0 | 0.125 | 1.0 | 0.0 | 1.4 | 1.6 | steady case
Table 1: Steady contact-discontinuity test case with perfect gas EOS Test case | $\rho_{L}$ | $p_{L}$ | $u_{L}$ | $\rho_{R}$ | $p_{R}$ | $u_{R}$ | $\gamma_{L}$ | $\gamma_{R}$ | Time
---|---|---|---|---|---|---|---|---|---
1 | 1.0 | 1.0 | -1.0 | 1.0 | 5.0 | 1.0 | 1.4 | 1.4 | 0.21
2 | 1.0 | 1000.0 | 0.0 | 0.125 | 1.0 | 0.0 | 1.6 | 1.4 | 0.21
Table 2: Initial conditions for shock tube problem and mass fraction
positivity test case: data obtained from 15, 8 with Perfect gas EOS

A fully coupled approach is used for the flux evaluation approach as explained
in 8 and it is stated that the Steger-Warming scheme and van Leer Scheme will
preserve the maximum of $0\leq Y\leq 1$ as they are using a fully coupled
approach.

###  3.1 Steady contact-discontinuity

This test case refers to a contact discontinuity wherein there is a jump in
density and $\gamma$ as given in table (1). The shock tube is filled with two
different perfect gases denoted by variable $\gamma$. The initial
discontinuity is present at $x=0.5$ with $x\in[0,1]$. A total of 100 control
volumes are used in simulations. Various numerical schemes like Steger
Warming, VanLeer, MOVERS-n,MOVERS-1,RICCA and MOVERS+ are compared for this
test case. The ability of the numerical schemes to resolve the steady contact-
discontinuity is analysed here.

Figure 2: Steady contact-discontinuity Comparison of RICCA,MOVERS+,MOVERS-n,
MOVERS-1,Vanleer and Steger Warming with unequal $\gamma$

Figures (2) refer to comparison of solution obtained for the steady contact-
discontinuity case as given in table (1). It can be observed that RICCA,
MOVERS+, MOVERS-1 and MOVERS-n resolve the steady contact exactly and even the
mass fraction is resolved exactly, whereas for Steger-Warming and van Leer
methods, the contact discontinuity and the mass fractions are diffused. It can
also be observed that in Steger-Warming scheme and van Leer scheme
oscillations are present in pressure and velocity but the positivity of the
mass fraction is preserved.

###  3.2 Pressure oscillations test case

This test case (2) is used by 16, to test the positivity of the mass fraction
by the regular Godunov type conservative finite volume methods. In 16, the
authors claim that many of the numerical methods which are formulated in the
conservative finite volume method would fail to preserve the positivity of the
mass fraction. Further, the pressure and velocity have oscillations for the
regular finite volume methods.

Figure 3: Isolated material front problem using RICCA and MOVERS+

Figures ( 3) refer to isolated material front test case as described in table
(2). As commented by the authors in 5, for this test case pressure
oscillations are present for all first-order numerical schemes which are
designed based on the conservative formulation with fully coupled approach and
the positivity of the mass fractions for such scheme is doubtful. Though the
new algorithms are based on a fully coupled approach in conservative
formulation with controlled numerical diffusion, they produce mild
oscillations in pressure but they preserve the positivity of the mass fraction
unlike the other conservative schemes.

###  3.3 Sod shock tube test case

This is a standard test case whose initial conditions are given in table (2)
for variable gamma values. The data values of the test case are taken from 15
also referred in 8. The second test case is a stiffer shock tube problem with
variation of pressure, mentioned as in test case 3 in table (2). For these
shock tube problems the initial discontinuity is located at $x=0.5$, with mass
fraction to the left of the discontinuity $Y_{L}=1,Y_{R}=0$ and to the right
of the discontinuity $Y_{L}=0,Y_{R}=1$. In all of the computations 100 equally
spaced control volumes on the interval $[0,1]$ are considered with CFL of
$0.45$ till the prescribed time is reached. Numerical results are presented
for Steger-Warming scheme, van Leer scheme, MOVERS-1, MOVERS-n, MOVERS+, and
RICCA.

Figures (4) refer to the standard Sod shock tube problem whose initial
conditions are defined as test case 2 in table (2). This shock tube problem
has two different fluids with different $\gamma$ values initially separated by
a diaphragm placed at $x=0.5$. Numerical simulations are carried out with 100
control volumes for all schemes and the reference solution is generated using
10000 control volumes using Rusanov Method.

Figure 4: Sod shock tube problem using RICCA with unequal $\gamma$,
$\gamma_{L}=1.6$, $\gamma_{R}=1.2$

Figures ( 5) refer to the stiff shock tube whose initial conditions are
defined as in test case 3 in table (2). As can be seen all the numerical
schemes preserve the mass fraction positivity and no pressure oscillations are
present in the pressure.

Figure 5: Stiff shock tube problem using RICCA,MOVERS+, MOVERS-n,MOVERS-1 with
equal gamma $\gamma=1.4$

###  3.4 1D results for $\gamma$-based model

Numerical simulations have been carried out for the test cases described in 34
whose initial conditions are given in the table (3).

Sno | $\rho_{L}$ | $p_{L}$ | $u_{L}$ | $\rho_{R}$ | $p_{R}$ | $u_{R}$ | $\gamma_{L}$ | $\gamma_{R}$ | ${p_{\infty}}_{L}$ | ${p_{\infty}}_{R}$ | Time
---|---|---|---|---|---|---|---|---|---|---|---
1 | 1.0 | 1.0 | 1.0 | 0.125 | 1.0 | 1.0 | 1.4 | 1.2 | 0 | 0 | 0.12
2 | 1.0 | 1.0 | 1.0 | 0.125 | 1.0 | 1.0 | 1.4 | 4.0 | 0 | 1 | 0.12
3 | 1.241 | 2.753 | 0.0 | 1.0 | $3.059\times 10^{-4}$ | 0.0 | 1.4 | 5.5 | 0 | 1.505 | 0.1
4 | 1.0 | 1.0 | 0.0 | 5.0 | 1.0 | 0.0 | 1.4 | 4.0 | 0 | 1 | 0.2
|  |  |  | 7.093 | 10.0 | -0.7288 |  | 4.0 | 0 | 1 | 0.2
Table 3: 1D Shock tube test cases referred from 34 with stiffened gas EOS

Test case 1 is an interface only problem and consists of a single contact
discontinuity. This test case consists of two sets of data

  1. 1.

a polytropic gas with two constants states as case 1,

  2. 2.

has same states except for the changes in $\gamma$ and $P_{\infty}$ as in case
2.

Initial position of the diaphragm is located at $x=0.2$ and the length of the
shock tube, $L=1$. Results are shown for 1O RICCA and MOVERS+. 100 control
volumes are considered for computation for both cases and the results are
shown at the prescribed time of 0.12.

Figure 6: Interface only problem with $P_{\infty}=0$ simulated using MOVERS+
and RICCA using eq(10) in 34

Figures (6,LABEL:InterfaceOnlyNoStiffgaswitheq10NWSC) refer to interface only
problem with $p_{\infty}=0$, which corresponds to perfect gas EOS using RICCA
and MOVERS+. It can be observed from the figures that the pressure
oscillations are not present when gamma based model (unlike in Figure (2.8))
is used and $\frac{1}{\gamma-1}$ is used as the conservative variable as
suggested by 3 and reconfirmed by 34. Shyue has also suggested that using
eq(10) in 34 cannot be generalised to all shock interaction problems. Figures
(7,LABEL:InterfaceOnlyStiffgaswitheq10NWSC) refer to interface only problem
with $P_{\infty}\neq 0$, which corresponds to stiffened gas EOS using RICCA
and MOVERS+ schemes. It can be observed from the figure that the pressure
oscillations are not present when gamma based model. Even for a large jump in
$\gamma$ both the numerical schemes do not generate any pressure oscillations.

Figure 7: Interface only problem with stiff gas simulated using RICCA using
eq(10) of 34

Test case 2 is a two fluid gas-liquid Riemann problem with initial conditions
as given in case 3 in the table. The diaphragm position is located with gas
occupying the domain till $x\leq 0.5$ and then the liquid and the time for
computation is $t=0.1$. In these test cases the conservative formulation using
$\frac{\rho}{\gamma-1}$ is used. Data for reference solution is taken from 34.
Simulations are shown for MOVERS+ and RICCA with100 control volumes. It can be
observed form the figures (LABEL:Liquid_Gas_RP_RICCA,8) and figure
(LABEL:Liquid_Gas_RP_Pinf) that there are no pressure oscillations present in
the pressure and the velocity. The internal energy and the $p_{\infty}$ are
accurately predicted.

Figure 8: Liquid gas Riemann problem using MOVERS+ and RICCA

The third test case considered here is a shock contact-discontinuity
interaction problem with the data given as in case 4. Here the two liquids are
separated by the interface at $x=0.5$ and a shock wave located at $x=0.6$ with
the pre- and post-shock conditions as given in the table 3 and the
computations are carried out for a time $t=0.2$. For this case 200 control
volumes are considered in the computation, the results from computation are
shown in figure (9). It can be seen that the phenomenon is captured accurately
by both RICCA and MOVERS+.

Figure 9: Sock contact interaction problem using RICCA and MOVERS+

###  3.5 Extension to two dimensions

In this section, the multicomponent model described in the previous section is
extended to 2D. The governing equations for $\gamma$-based model in 2D are
given by (3.1)

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\left(\partial\rho
u\right)}{\partial x}+\frac{\left(\partial\rho v\right)}{\partial y}$
$\displaystyle=0$ (3.1) $\displaystyle\frac{\partial\rho u}{\partial
t}+\frac{\partial\left(\rho u^{2}+p\right)}{\partial
x}+\frac{\left(\partial\rho uv\right)}{\partial y}$ $\displaystyle=0$ (3.2)
$\displaystyle\frac{\partial\rho v}{\partial t}+\frac{\left(\partial\rho
uv\right)}{\partial x}+\frac{\partial\left(\rho v^{2}+p\right)}{\partial y}$
$\displaystyle=0$ (3.3) $\displaystyle\frac{\partial\rho E_{t}}{\partial
t}+\frac{\partial\left[\left(\rho E_{t}+p\right)u\right]}{\partial
x}+\frac{\partial\left[\left(\rho E_{t}+p\right)v\right]}{\partial y}$
$\displaystyle=0$ (3.4)

with the equations $\gamma$ given by

$\displaystyle\frac{\left(\partial\frac{\rho}{\left(\gamma-1\right)}\right)}{\partial
t}+\frac{\left(\partial\frac{\rho u}{\left(\gamma-1\right)}\right)}{\partial
x}+\frac{\left(\partial\frac{\rho v}{\left(\gamma-1\right)}\right)}{\partial
x}$ $\displaystyle=0$ (3.5) $\displaystyle\frac{\left(\partial\frac{\rho\gamma
p_{\infty}}{\left(\gamma-1\right)}\right)}{\partial
t}+\frac{\left(\partial\frac{\rho\gamma
P_{\infty}u}{\left(\gamma-1\right)}\right)}{\partial
x}+\frac{\left(\partial\frac{\rho\gamma
p_{\infty}v}{\left(\gamma-1\right)}\right)}{\partial y}$ $\displaystyle=0$
(3.6)

and stiffened gas EOS as given by (3.7).

$\displaystyle\frac{p+\gamma p_{\infty}}{\gamma-1}=\rho e$ (3.7)

For the 2D $\gamma$-based model given above, simulations are carried out using
RICCA and MOVERS+.

####  3.5.1 Moving Interface Problem

The first test case considered is a moving interface problem which consists of
a bubble with radius $r_{0}=0.16$ evolving in a constant velocity field
$\left(u,v\right)=(1,1)$. The initial data considered here is similar to the
1D test case as described in table (3). The pressure is uniform with value
$p=1$ while the $\rho,\gamma,p_{\infty}$ jump across the interface. Initially
the bubble is placed at $x_{c}=0.25,y_{c}=0.25$ on a domain which varies from
$x\in[0,1]$ and $y\in[0,1]$. A total of $100\times 100$ control volumes are
considered in $x$ and $y$ directions and the solution is evolved using a time
accurate scheme till $t=0.36$.

RICCA

MOVERS+

Figure 10: Surface of density at $t=0$ and $t=0.36$

Figure 11: Density contours at $t=0$ and $t=0.36$

Figure 12: Pressure distribution in the domain Figure 13: Contour and surface
view of density and pressure of interface only at $t=0.36$ on $100\times 100$
Grid

RICCA

MOVERS+

Figure 14: Surface of density at $t=0$ and $t=0.36$

Figure 15: Density contours at $t=0$ and $t=0.36$

Figure 16: Pressure distribution in the domain Figure 17: Contour and surface
view of density and pressure of interface only at $t=0.36$ on $500\times 500$
Grid Figure 18: Line plot of density across the interface Figure 19: Line plot
of pressure across the interface Figure 20: 2D Line plots of density and
pressure across the interface using RICCA Scheme on $500\times 500$ and
$100\times 100$ Grid

Figure (13) refers to the moving bubble at $t=0$ and at $t=0.36$ using RICCA
and MOVERS+ on a $100\times 100$ grid and Figure (17) refers to results
obtained on a $500\times 500$ grid. It can be observed that the position of
the bubble is captured accurately and the pressure does not have any
oscillations. Further the pressure and density plots across the bubble are
shown in the figure (20) where in it can be observed that there are no
pressure oscillations in the results generated by both RICCA and MOVERS+.

####  3.5.2 Bubble explosion problem

The second test case considered is a radially symmetric problem. It consists
of a circular bubble present initially at rest in water and suddenly explodes
due to high pressure of the water. The bubble is placed at $(x,y)=(0.5,0.5)$
and has a radius of $r_{o}=0.2$. The fluid inside the bubble has the following
initial conditions $(\rho,p,\gamma,p_{\infty})=(1.241,2.753,1.4,0)$ and the
surrounding water has the following properties
$(\rho,p,\gamma,p_{\infty})=(0.991,3.059\times 10^{-4},5.5,1.505)$. Numerical
simulations are carried out using RICCA and MOVERS+. Pressure and density
contour plots are shown in figure (25) and the variation of pressure and
density at $y=0.5$ are shown in figure (25).

MOVERS+

RICCA

Figure 21: Surface pressure of bubble explosion

Figure 22: Surface density of bubble explosion Figure 23: Contour and surface
view of density and pressure of bubble explosion under water at $t=0.058$
using MOVERS+ and RICCA on $500\times 500$ Grid

Figure 24: 2D Line and contour plots of density and ressure using MOVERS+

Figure 25: 2D Line and contour plots of density and pressure using RICCA

##  4 Conclusions

Numerical simulations of Euler equations in 1D and 2D have been carried out
using MOVERS (both scalar dissipation and vector dissipation), MOVERS+ and
RICCA, for multicomponent gases with perfect gas EOS and stiffened gas EOS
using mass fraction approach and $\gamma$-based approach in conservative form.
It can be concluded that

  1. 1.

both the numerical schemes RICCA and MOVERS+ can be extended to multicomponent
gases with different EOS without any modifications,

  2. 2.

both the numerical schemes preserve the mass fraction positivity and the
pressure positivity in the conservative approach when used in mass fraction
based model.

  3. 3.

pressure oscillations are observed in the finite volume framework for
interface only problem and when $\gamma$-based model is adopted no pressure
oscillations are observed.

  4. 4.

both these schemes can be easily extended to any number of components.

  5. 5.

both the schemes can handle large jumps in $\gamma$ without any modifications.

## References

  * 1
  * 2 R. Abgrall, An extension of Roe’s upwind scheme to algebraic equilibrium real gas models, Computers & Fluids 19(2) (1991), p.p. 171-182.
  * 3 R. Abgrall, How to prevent pressure oscillations in multicomponent flow calculations: A quasi conservative approach, J. Comput. Phys. 125 (1996), p.p. 150-160.
  * 4 R. Abgrall, Generalisation of the Roe scheme for the computation of mixture of perfect gases,Rech. Aerospat., 6 (1988), pp. 31-43 (English edition).
  * 5 R. Abgrall, S. Karni, Computations of Compressible Multifluids, Journal of Computational Physics 169 ( 2001), p.p. 594-623.
  * 6 R. Abgrall, Preliminary results on an extension of Roe approximate Riemann solver to non equilibrium flows, [Research Report] RR-0987 (1989), pp.39. <inria-00075572>
  * 7 J. Blazek, Computational Fluid Dynamics Principles and Applications, Elsevier Publications (2005).
  * 8 D. Chargy, R. Abgrall, L. Fezoui, B. Larrouturou, Conservative numerical schemes for multicomponent inviscid flows, Rech. Aerospat., 2 (1992), pp. 61-79 (English version).
  * 9 P. Colella, H. Glaz, Efficient Solution Algorithms for the Riemann Problem for Real gases, Journal of Computational Physics 59 (1983), p.p. 264-289.
  * 10 R.M.L. Coelho, P.L.C. Lage and A. Silva Telles, A Comparison of hyperbolic Solvers for Ideal and Real gas flows, Brazilian Journal of Chemical Engineering, 23, p.p. 301-318.
  * 11 T.J. Chung, Computational Fluid Dynamics, Cambridge University Press (2010).
  * 12 G. Fernandez, B. Larrouturou, Hyperbolic schemes for multicomponent EULER equations, in Nonlinear Hyperbolic Equations Theory, Computation Methods, and Applications, Vol. 24 of the series Notes on Numerical Fluid Mechanics, p.p. 128-138.
  * 13 S. Jaisankar, S.V. Raghurama Rao, A central Rankine-Hugoniot solver for hyperbolic conservation laws, Journal of Computational Physics, 228(3)(2009), p.p. 770-798.
  * 14 S.Jaisankar and S.V. Raghurama Rao, Diffusion regulation for Euler solvers, Journal of Computational Physics, vol. 221 (2007), pp. 577-599.
  * 15 S. Karni, Multicomponent flow calculations by a consistent primitive algorithm, Journal of Computational Physics., 112 (1994), pp. 31-43.
  * 16 S. Karni, S. Canic,Computations of Slowly Moving Shocks, Journal of Computational Physics 136 (1997), p.p. 132–139.
  * 17 A. Kurganov and E. Tadmor, New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations, Journal of Computational Physics 160(1) (2000), p.p. 241-282.
  * 18 B. Larrouturou, L. Fezoui, On the equations of multicomponent perfect and real gas inviscid flow, in Non-Linear Hyperbolic Problems, edited by Carasso, Charrier, and Joly, Lecture Notes in Mathematics, 1402,(Springer Verlag, Heidelberg, 1989) p.p. 69.
  * 19 B. Larrouturou, How to preserve the mass fractions positivity when computing compressible multi-component flows, [Research Report] RR-1080 (1989). <inria-00075479>
  * 20 M.S. Liou, Mass flux schemes and connection to shock instability, Journal of Computational Physics 160(2)(2000), p.p. 623-648.
  * 21 Liou, M.-S., and Steffen, C., A New Flux Splitting Scheme, Journal of Computational. Physics 107(1993), p.p 23-39.
  * 22 Liou, M.-S., A Sequel to AUSM: AUSM+ Journal of Computational Physics, 129(1996), p.p. 364-382.
  * 23 Liou, M.-S., A Sequel to AUSM, Part II: AUSM+-up, Journal of Computational Physics, 214(2006), p.p. 137- 170.
  * 24 P.D. Lax, weak solutions of nonlinear Hyperbolic Equations and their Numerical Computaataions, Comm. Pure. Appl. Math. VII (1954), p.p. 159-193.
  * 25 R.S. Lagumbay, Modeling and Simulation of Multiphase/Multicomponent Flows, Ph.D Thesis, University of Colorado, 2006.
  * 26 Maruthi N.H., Hybrid Central Solvers for Hyperbolic Conservation Laws, PhD Thesis, Indian Institute of Science, Bangalore, India, 2016.
  * 27 R. Menikoff, B.J. Plohr, The Riemann problem for fluid flow of real materials, Reviews of modern physics 61(1)(1989), p.p. 75.
  * 28 Ramesh Kolluru Novel,Robust and Accurate Central solvers for Real,Dense and Multicomponent gases, PhD Thesis, Indian Institute of Science, Bangalore, India, 2019.
  * 29 P.L. Roe, Journal of Computational Physics 43 (1981), p.p. 357.
  * 30 P.L. Roe, Characteristic-based schemes for the Euler equations, Annural Review of Fluid Mechanics, vol. 18 (1986), pp. 337-365.
  * 31 P.L. Roe, Shock Capturing, Chapter 6, Handbook of Shock Waves, , Academic Press, vol. I (2001), pp. 787-877
  * 32 S.V. Raghurama Rao, S. Deshpande, Computational Fluid Dynamics Journal of Japan Society of CFD 4 (1995), p.p. 415.
  * 33 V.V. Rusanov, Calculation of interaction of non steady shock waves with obstacles, NRC, Division of Mechanical Engineering (1962).
  * 34 K.M. Shyue, An Efficient Shock-Capturing Algorithm for Compressible Multicomponent Problems, Journal of Computational Physics 142 (1998), p.p. 208–242.
  * 35 K.M. Shyue, A Fluid-Mixture Type Algorithm for Compressible Multicomponent Flow with van der Waals Equation of State, Journal of Computational Physics 156 (1999), p.p. 43–88.
  * 36 K.M. Shyue, A Fluid-Mixture Type Algorithm for Compressible Multicomponent Flow with Mie Gruneisen Equation of State, Journal of Computational Physics 171 (2001), p.p. 678–707.
  * 37 K.M. Shyue, A fluid-mixture type algorithm for barotropic two fluid flow problems, Journal of Computational Physics 200 (2004), p.p. 718–748.
  * 38 R. Saurel, R. Abgrall, A simple method for compressible multifluid flows, SIAM Journal on Scientific Computing 21(3) (1999), p.p. 1115-1145.
  * 39 R. Saurel, M. Larini and J.C. Loraud, Exact and Approximate Riemann Solvers for Real Gases, J. Comput. Phys., 112 (1994), pp. 126-137.
  * 40 B. van Leer, Flux-vector splitting for the Euler Equations, in $8^{th}$ International Conference on Numerical Methods in Fluid Dynamics, Springer (1982), p.p. 507-512.
  * 41 B. van Leer, Upwind and high-resolution methods for compressible flow: From donor-cell to residual distribution schemes, Communications in Com- putational Physics, vol. 1, no. 2 (2006), pp. 192-206.
  * 42 N. Venkata Raghavendra, Discrete Velocity Boltzmann Schemes for Inviscid Compressible Flows, PhD Thesis, Indian Institute of Science, Bangalore, India, 2017.
  * 43 N. Venkata Raghavendra, S. V. Raghurama Rao, A Boltzmann scheme with physically relevant discrete velocities for Euler equations. arXiv:1612.07911v1 [physics.comp-ph] 23 Dec 2016.


#  XMM-Newton and Swift Observations of the Seyfert 1 AGN NGC 5940

Oluwashina K. Adegoke11affiliationmark:

###### Abstract

Probing the physics of the accretion flow around active galactic nuclei (AGN)
is crucial to understanding their emission mechanisms as well as being able to
constrain the geometrical and variability properties of the different regions
around them. The soft X-ray excess – usually observed below $\sim
2\,\mathrm{keV}$ in excess of the dominant X-ray powerlaw continuum – is one
prominent feature that is commonly seen in type 1 Seyfert AGN and therefore
readily provides a useful diagnostic of the accretion flow mechanism around
these systems. NGC 5940 is a Seyfert 1 AGN which reveals strong, prominent
soft X-ray excess below $\sim 2\,\mathrm{keV}$ as seen in both its XMM-Newton
and Swift observations. Model fit to the data revealed that this feature could
be equally well explained by the ionised partial covering, the thermal
Comptonisation and the blurred reflection models. Although the other models
cannot be decisively ruled out with the data at hand, the lack of significant
broad iron $K_{\alpha}$ as well as any significant emission/absorption line
features in the reflection grating spectrometer (RGS) data tend to favour the
thermal Comptonisation origin for the soft X-ray excess in NGC 5940.

††slugcomment: Not to appear in Nonlearned J., 45.00footnotetext: Department
of Mathematical & Physical Sciences, Afe Babalola University, P.M.B 5454, Ado-
Ekiti, Nigeria.

Keywords accretion, accretion discs – galaxies: active – galaxies: individual:
NGC 5940 – galaxies: nuclei – X-rays: galaxies.

##  1 Introduction

Active galactic nuclei (AGN) have been known to be luminous in all bands of
the electromagnetic spectrum. The optical/UV emission of AGN is believed to be
produced in an optically thick standard accretion disc (Shakura and Sunyaev,
1973; Novikov and Thorne, 1973) due to viscous heating as material in the disc
spiral down towards the central blackhole. The dominant X-ray powerlaw
emission – which is a significant fraction of an AGN’s bolometric luminosity –
is theorised to be produced from the Compton-upscattering of seed disc
optical/UV photons in a compact (size $\sim 10-20\,R_{g}$, where $R_{g}$ is
the gravitational radius), hot (having electron temperature $KT_{e}\sim
100\,\mathrm{keV}$), optically thin (with optical depth $\tau<1$) electron
plasma, called the corona (Haardt and Maraschi, 1993; Adegoke et al., 2019).
Other complex features imprinted on the X-ray spectra of AGN include the soft
X-ray excess emission below $\sim 2\,\mathrm{keV}$, iron $K_{\alpha}$ emission
complex at $\sim 6.4\,\mathrm{keV}$ and a Compton reflection “hump” which
peaks normally at $\sim 20\,\mathrm{keV}$.

The soft X-ray excess is a smooth emission component in the
$0.3-2.0\,\mathrm{keV}$ energy range that is commonly observed in Seyfert 1
AGN. It was first reported over three decades ago and was thought at the time
to be the Wien’s tail of the multi-temperature disc blackbody emission
observed in the UV band (see e.g., Singh et al., 1985; Arnaud et al., 1985;
Pounds et al., 1986; Walter and Fink, 1993). However, it was soon shown that
the predicted temperature ($\sim 0.1-0.2\,\mathrm{keV}$) is much higher than
expected for a standard accretion disc around a supermassive blackhole. It was
equally found that the temperature is fairly constant and not related to the
luminosity and mass of the central blackhole as against expectations that disc
temperatures should be correlated with blackhole masses (see e.g., Gierliński
and Done, 2004; Piconcelli et al., 2005; Miniutti et al., 2009). Gierliński
and Done (2004) proposed that the soft excess could be an artefact of smeared
absorption from partially ionised material, the model however requires
unphysically large smearing velocity ($\sim 0.5c$, where $c$ is the speed of
light) for disc wind which has been shown to be difficult to generate from
hydrodynamic simulations (Schurch and Done, 2007). Thus, the origin of the
soft excess is still being debated even today. Most recently, two competing
models have tried to explain the origin of the soft X-ray excess – the thermal
Comptonisation and the relativistically blurred reflection models. The thermal
Comptonisation model posits that the soft X-ray excess results from the
inverse-Compton scattering of disc optical/UV photons in an optically thick
($\tau\sim 10$), warm ($KT_{e}\sim 0.1-1\mathrm{keV}$) corona (e.g., Magdziarz
et al., 1998; Done et al., 2012; Adegoke et al., 2017), different from the hot
($\sim 100\,\mathrm{keV}$) corona responsible for the dominant X-ray powerlaw.
The relativistic reflection model on the other hand explains the soft X-ray
excess to be the result of blurred reflection from partially ionised disc
material in the vicinity of the blackhole due to its strong gravity and the
high velocities of the accreting material (e.g., Crummy et al., 2006; García
et al., 2014; Dauser et al., 2014). While both models tend to fit the data
equally well in many cases, testing these models for a large number of sources
and with data covering a wide range of the X-ray and other bands will be
crucial to narrowing down on the actual origin of the soft X-ray excess as
well as its possible dependence on the accretion flow properties of individual
systems (see e.g., Jin et al., 2012; Petrucci et al., 2018; Gliozzi and
Williams, 2020). For such a large sample, a systematic search for the presence
or absence of other features such as emission/absorption lines will help in
identifying the preferred models. It is equally plausible that both thermal
Comptonisation and X-ray reflection contribute simultaneously to the overall
broadband X-ray spectra in AGN such that while the Compton reflection hump is
due to relativistic reflection, thermal Comptonisation produces the soft X-ray
excess, as was recently reported for the Seyfert 1 AGN Ton S180 (Matzeu et
al., 2020).

NGC 5940 is a nearby Seyfert 1 AGN at a redshift of $0.034$ and luminosity
distance of $147\,\mathrm{MPc}$ (Robinson et al., 2019). Its mass has been
estimated to be $1.085^{+0.182}_{-0.167}\times 10^{7}\,M_{\odot}$ in the
reverberation mapping database111http://www.astro.gsu.edu/AGNmass/ of Bentz
and Katz (2015). This is based on the 2011 Lick Observatory $\rm{H}\beta$ and
$\rm{Fe}II$ lag measurements reported in Barth et al. (2013, 2015). Despite
its relative proximity, its spectral properties have not been probed in detail
until now, especially in the X-ray energy band.

This paper presents a detailed analysis of the XMM-Newton X-ray data of NGC
5940 – supplemented with data from the Swift X-ray Telescope. In Section 2, we
describe the observation and data reduction procedure. In Section 3, we detail
the data analysis approach used. We discuss the data analysis result in
Section 4, and in Section 5, we draw conclusions and summarise our results.

##  2 Observations and Data Reduction

NGC 5940 was observed by XMM-Newton (Jansen et al., 2001) on the 23rd of
February, 2012 for a duration of $\sim 33\,\mathrm{ks}$ (ObsID: 0670040601).
All cameras onboard the telescope were used during the observation. The
European Photon Imaging Cameras (EPIC) (Strüder et al., 2001; Turner et al.,
2001) were operated in the full frame mode. While EPIC-pn was operated with
thin filter, both MOS1 and MOS2 were operated with medium filters. The optical
monitor (OM) (Mason et al., 2001) was operated in the image mode using the
UVM2 and UVW1 filters.

The data was processed using the Science Analysis System (sas v.17.0.0) with
updated Current Calibration Files (CCF). The sas task EVSELECT was employed to
extract the event list. The data were subsequently screened and filtered for
intervals of high particle background after which a good time interval (GTI)
file was created. To check for possible pile-up in the data, a filtered event
list was extracted from the raw event list – from a radius of 40 arcseconds
centred on the source coordinates – using the GTI file after which the
EPATPLOT command was applied. The check revealed the presence of significant
pile-up in the pn data. To correct for the pile-up, X-ray loading was first
corrected for, – since X-ray loading normally accompanies piled-up
observations – after which the EVSELECT command was applied on the new event
list using the GTI file to create a cleaned event file. Source photons were
extracted from a region of radius 40 arcseconds centred on the source position
while background was extracted from a circular source-free region of radius 80
arcseconds on the same CCD. The EVSELECT task was used to generate the X-ray
source and background spectra. The SAS tasks ARFGEN and RMFGEN were
subsequently used to compute the ancillary and photon redistribution matrices.
Pile-up correction was implemented in the generation of the redistribution
matrix which was calculated from the frequency and spectrum of the incoming
photons. The resulting spectra were then grouped using the SPECGROUP command
to have a minimum of 25 counts per bin to facilitate use of the $\chi^{2}$
minimisation technique.

The reduction of the reflection grating spectrometer (RGS) (den Herder et al.,
2001) data followed standard procedure employing the SAS task RGSPROC. The
accuracy of the source coordinates was checked from the source list and found
to be in order. The data were subsequently checked against periods of high
particle background by plotting the lightcurve of the pure background events.
For improved signal-to-noise ratio, first order RGS1 and RGS2 spectra were
combined using the task RGSCOMBINE and then grouped to have a minimum of 30
counts per bin for use with XSPEC $\chi^{2}$ statistics.

Swift (Gehrels et al., 2004) has observed NGC 5940 a couple of times since
2008. Although none of the snapshots were simultaneous with the 2012 XMM-
Newton observation of the source, four observations carried out between 2010
and 2011 were combined into one spectra of total exposure $\sim
9.8\,\mathrm{ks}$ and used for spectral analysis. This is based on the
assumption that the spectral properties of the source had not changed
significantly between the time of the 2012 XMM-Newton and the earlier
2010/2011 Swift observations. The Swift X-ray Telescope (XRT; Burrows et al.,
2005) data reduction was carried out with the aid of the online Swift/XRT
product generator hosted at the University of
Leicester222https://www.swift.ac.uk/index.php (Evans et al., 2009).

##  3 Spectral Analysis

We used the EPIC-pn X-ray data for our spectral analysis due to its higher
signal-to-noise ratio compared to EPIC-MOS. We carried out the spectral
analysis on NGC 5940 using XSPEC v.12.11.1 (Arnaud, 1996). We used the
$\chi^{2}$ statistics and for a single parameter of interest, errors were
quoted at the 90% confidence level – corresponding to $\Delta\chi^{2}=2.706$ –
unless otherwise stated. In all the fitting, we have added contributions from
galactic absorption column of $N_{H}=3.6\times 10^{20}\,\mathrm{cm^{-2}}$
obtained from the recent HI21cm measurements (HI4PI Collaboration et al.,
2016), modelled with tbabs (Wilms et al., 2000).

Fig. 1 : Residual plot revealing strong soft X-ray excess below
$2\,\mathrm{keV}$ based on a powerlaw fit to the $2-10\,\mathrm{keV}$ EPIC-pn
data of NGC 5940 extrapolated to $0.3\,\mathrm{keV}$.

We began by fitting the $2-10\,\mathrm{keV}$ band with a baseline model
consisting of a simple powerlaw modified by galactic absorption (i.e.
tbabs*powerlaw). This provided a very good fit to the data with photon index
$\Gamma=1.78$ and $\chi^{2}/dof$ of 91.7/93, where $dof$ represents the
degrees of freedom. Inclusion of a gaussian line to model possible iron
$K_{\alpha}$ features did not improve the fit, giving $\chi^{2}/dof=88.3/90$.
When extrapolated down to $0.3\,\mathrm{keV}$, the data revealed strong soft
X-ray emission in the source below $\sim 2\,\mathrm{keV}$ in excess of the
powerlaw as shown in the residual plot of Fig. 1, providing an unacceptable
fit with $\chi^{2}/dof$ of 223.6/140. To test existing models on the origin of
this soft excess emission, we fitted the entire $0.3-10\,\mathrm{keV}$ band
with phenomenological and subsequently physically motivated models.

###  3.1 Phenomenological Model

To model the soft X-ray excess, we added the multicolour disc blackbody model
diskbb (see e.g., Mitsuda et al., 1984; Makishima et al., 1986) to the
baseline model (i.e. tbabs*(diskbb+powerlaw)). This provided an acceptable fit
to the data with $\chi^{2}/dof=159.7/138$. The inner disc temperature
$T_{in}\sim 0.18\,\mathrm{keV}$ is high for an accretion disc around a
$10^{7}\,M_{\odot}$ blackhole. The photon index $\Gamma$ is 1.89. Adding a
gaussian line at $6.4\,\mathrm{keV}$ to the model slightly improved the fit
with $\Delta\chi^{2}/\Delta dof=-12/-3$, giving $\chi^{2}/dof=147.4/135$.
However, in the subsequent models, addition of a gaussian line at $\sim
6.4\,\mathrm{keV}$ did not improve the fits in any way and was therefore
excluded.

Fig. 2 : Unfolded spectra of best fit for each of the models applied to the
$0.3-10\,\mathrm{keV}$ XMM-Newton EPIC-pn data of NGC 5940. The upper panels
in each plot show the model-data fit – the black solid lines represent the
composite model in each of the plots. The lower panels in each plot show the
residuals.

###  3.2 Physical Models

####  3.2.1 Ionised partial covering: zxipcf

The shape of the residual plot for the powerlaw fit to the
$0.3-10\,\mathrm{keV}$ data seem to show features that may result from complex
partial covering absorption along the line of sight (shown in Fig. 2 – top,
leftmost plot). In this scenario, the intrinsic powerlaw continuum is obscured
by an ionised absorber along the line of sight. The observed spectrum is then
a combination of the obscured and the direct powerlaw components (see e.g.,
Reeves et al., 2008; Gallo et al., 2011; Tripathi et al., 2019). We tested
this possibility by multiplying the absorbed powerlaw with the ionised partial
covering absorption model zxipcf (i.e. tbabs*zxipcf*powerlaw). This provided a
very good fit to the data with $\chi^{2}/dof=135.9/137$, covering fraction
$C_{f}=0.42$ and line of sight column density $N_{h}\sim 1.6\times
10^{23}\,\mathrm{cm^{-2}}$, the best-fit parameters are shown in Table 1.

Table 1 : The best-fit parameter values from fitting different models to the
$0.3-10\,\mathrm{keV}$ XMM-Newton EPIC-pn data of NGC 5940.

Model/Parameter & Best-fit values |
---|---
Model | tbabs*(diskbb+powerlaw)
$T_{\mathrm{in}}(\mathrm{keV})$ | $0.185\pm{0.02}$
Normalisation - diskbb | $51.28^{+24.41}_{-15.51}$
Photon index $\Gamma$ | $1.89\pm 0.03$
Normalisation ($10^{-3}$) - powerlaw | $1.35\pm 0.05$
$\chi^{2}/dof$ | 159.7/138
Model | tbabs*(zxipcf*powerlaw)
Column density $N_{h}$ ($10^{22}$) - zxipcf | $15.74^{+24.49}_{-9.96}$
Ionisation parameter log($\xi$) | $-0.20^{+1.41}_{-P}$
Covering fraction $C_{Frac}$ | $0.42^{+0.15}_{-0.09}$
Photon index $\Gamma$ | $2.09\pm 0.02$
Normalisation ($10^{-3}$) - powerlaw | $2.62^{+0.96}_{-0.34}$
$\chi^{2}/dof$ | 135.9/137
Model | tbabs*(optxagnf)
Mass ($\times 10^{7}\,M_{\odot}$) |  $1.1$(f)
Luminosity distance ($\mathrm{MPc}$) |  $147$(f)
log($\frac{L}{L_{Edd}}$) | $-1.26^{+0.60}_{-0.06}$
spin a* |  $0.998$(f)
Corona radius $R_{cor}$ | $18.12^{+P}_{-11.14}$
Warm corona $kT_{SE}$ ($\mathrm{keV}$) | $0.40^{+0.14}_{-0.09}$
Optical depth $\tau$ | $12.92^{+2.55}_{-2.05}$
Photon index $\Gamma$ | $1.71^{+0.09}_{-0.12}$
Powerlaw fraction $f_{pl}$ | $0.77^{+0.08}_{-0.21}$
$\chi^{2}/dof$ | 135.6/136
Model | tbabs*(relxillcp)
spin a* |  $0.998$(f)
Inclination |  $30\deg$(f)
Photon index $\Gamma$ | $1.92\pm 0.02$
Ionisation parameter log($\xi$) | $3.00^{+0.10}_{-0.15}$
Iron abundance $A_{Fe}$ |  $1.0$(f)
Hot corona $kT_{e}$ ($\mathrm{keV}$) |  $100$(f)
Reflection frac | $0.55^{+0.16}_{-0.13}$
Normalisation ($10^{-5}$) | $1.89^{+0.11}_{-0.14}$
$\chi^{2}/dof$ | 164/138
Model | tbabs*(relxill+nthcomp)
spin a* | $0.991^{+P}_{-0.418}$
Inclination |  $30\deg$(f)
Photon index $\Gamma$\- relxill | $2.12^{+0.01}_{-0.04}$
Ionisation parameter log($\xi$) | $0.02^{+0.30}_{-P}$
Iron abundance $A_{Fe}$ |  $1.0$(f)
High energy cutoff $E_{cut}$ ($\mathrm{keV}$) |  $300$(f)
Reflection frac $R_{frac}$ |  $-1$(f)
Normalisation ($10^{-5}$) - relxill | $5.65^{+0.92}_{-1.57}$
Photon index $\Gamma$ \- nthcomp |  $2.12$(t)
Hot corona $kT_{e}$ ($\mathrm{keV}$) |  $100$(f)
Disk blackbody $kT_{bb}$ ($\mathrm{keV}$) |  $0.1$(f)
Normalisation ($10^{-3}$) - nthcomp | $1.49\pm 0.02$
$\chi^{2}/dof$ | 134.1/138

Note: “f” implies that the value is frozen, “t” implies that the parameter
value is tied to another parameter and “P” implies that the error value is
pegged at the hard limit.

####  3.2.2 Warm corona model: optxagnf

The warm coronal model (Done et al., 2012) posits that the gravitational
energy released at each annulus of the disc is radiated as thermal blackbody
emission down to the radius of the corona $R_{cor}$. The gravitational
potential energy within this radius is assumed to no longer be completely
thermalised, giving rise to a warm ($kT_{SE}\sim 0.2\,\mathrm{keV}$),
optically thick ($\tau>1$) plasma responsible for the soft X-ray excess
emission and a hot ($kT_{e}\sim 100\,\mathrm{keV}$), optically thin ($\tau<1$)
corona responsible for the dominant X-ray powerlaw continuum above $\sim
2\,\mathrm{keV}$. Important parameters of the model include the blackhole mass
$M_{BH}$, the dimensionless spin a*, the electron temperature $kT_{SE}$ for
the soft Comptonisation component and its optical depth $\tau$, the photon
index of the powerlaw X-ray continuum $\Gamma$, the Eddington ratio in log
scale log($L/L_{Edd}$) and the luminosity distance $D_{L}$ to the source. In
fitting with this model, we fixed the blackhole mass at $\sim 1.1\times
10^{7}\,M_{\odot}$ (Robinson et al., 2019), the disc outer radius at its
default value and $D_{L}$ at $147\,\mathrm{MPc}$. Following Crummy et al.
(2005), we fixed and froze the spin parameter to the maximum allowed value of
a*$=0.998$ as it was found to be unconstrained when freed. This model resulted
in an excellent fit with $\chi^{2}/dof=135.6/136$. The best-fit values of the
Eddington ratio, $R_{cor}$, $kT_{SE}$, optical depth and the powerlaw photon
index are $0.05$, $18.08\,R_{g}$, $0.40\,\mathrm{keV}$, $12.92$ and $1.71$
respectively. All the fit parameters are shown in Table 1.

####  3.2.3 Blurred reflection: relxill

To a good extent, reflection models have been successful in modelling some of
the prominent features imprinted on the X-ray spectrum of AGN including the
soft X-ray excess, the hard X-ray excess as well as line profiles like the
iron $K_{\alpha}$ emission feature seen at $\sim 6.4\,\mathrm{keV}$ (see e.g.,
Fabian et al., 2004; Walton et al., 2013; Jiang et al., 2019). To test whether
blurred relativistic reflection in the immediate vicinity of the AGN can
explain the origin of the prominent soft X-ray excess emission, we used the
relxill model. This model combines the capabilities of the reflection code
xillver (García et al., 2014) with the relativistic ray-tracing code relline
(Dauser et al., 2014). It models the irradiation of the accretion disc by a
broken powerlaw emissivity of the form $\epsilon\propto r^{-q_{1}}$ between
$r_{in}$ and $r_{br}$ and $\epsilon\propto r^{-q_{2}}$ between $r_{br}$ and
$r_{out}$; where $q_{1}$ and $q_{2}$ are the inner and outer emissivity
indices, $r$ is the radius of the accretion disc while $r_{in}$, $r_{br}$ and
$r_{out}$ are the inner, break and outer radii of the accretion disc
respectively. More information about the model and its different flavours are
available on its documentation webpage333http://www.sternwarte.uni-
erlangen.de/~dauser/research/relxill/. We started by applying the relxillcp
flavour of the model which is more physically consistent as the reflection
fraction in this case is computed using the primary continuum implemented with
the nthcomp Comptonisation model (Zdziarski et al., 1996; Życki et al., 1999).

During fitting, the indices $q_{1}$ and $q_{2}$ were kept at their default
values of $3$. Also, $r_{br}$, inclination, $r_{in}$ and $r_{out}$ were left
frozen at their default values. Iron abundance relative to solar abundance
$A_{Fe}$ and the hot corona temperature were fixed at $1$ and
$100\,\mathrm{keV}$ respectively. The spin parameter a* was fixed at its
maximum value of 0.998. The model provided an acceptable fit to the data with
$\chi^{2}/dof=164/138$ and a steeper photon index $\Gamma=1.92$. Allowing the
inner emissivity index $q_{1}$, break radius $r_{br}$, spin, inclination and
the iron abundance $A_{Fe}$ parameters to be free improved the fit
significantly, with $\Delta\chi^{2}/\Delta dof=-20/-4$ giving
$\chi^{2}/dof=144/133$. This configuration additionally gave $r_{br}\sim
4\,R_{g}$, $q_{1}\sim 9.8$, a*$=0.899$, $\Gamma=1.96$, inclination $\sim
35\deg$ and $A_{Fe}\sim 4$. As evident, the inner emissivity index is rather
extreme in this case which may indicate strong inner disc reflection due to
strong general relativistic effects in the immediate vicinity of the
blackhole. The modelled iron abundance is also considerably super-solar.

We afterwards replaced the relxillcp with the relxill+nthcomp model where the
illuminating powerlaw is modelled with nthcomp and the blurred reflection only
with relxill. For this, the reflection fraction parameter $R_{frac}$ was set
to $-1$. Most parameters of the relxill and relxillcp models are the same with
the exception of the coronal electron temperature $kT_{e}$ in relxillcp which
is replaced by the high energy cutoff $E_{cut}$ in relxill. The photon index
of nthcomp was tied to that of the relxill model. This provided an excellent
fit to the data with $\Delta\chi^{2}=-10$ compared to the previous model for
four additional degrees of freedom, giving $\chi^{2}/dof=134.1/137$ and a
slightly steeper photon index value $\Gamma=2.12$. Table 1 contains all the
fitted parameter values of the model.

###  3.3 RGS spectral analysis

We analysed the RGS spectra of NGC 5940 to check for possible prominent
emission and absorption line features in the $0.3-2\,\mathrm{keV}$ energy
range associated with the soft X-ray excess. This is because the RGS has a
much higher energy resolution compared to both EPIC-pn and MOS and as such
these reflection features should be evident in the RGS spectra of the source
if blurred relativistic reflection is predominantly what produces the soft
X-ray excess. We modelled the combined first order RGS spectra with a simple
powerlaw modified by galactic absorption. This provided a very good fit with
$\chi^{2}/dof=216.9/221$ and $\Gamma=2.12$. Figure 3 shows the powerlaw
modelled RGS (1+2) stacked spectra of NGC 5940, with the residual shown in the
lower panel of the plot. This spectral analysis result suggested that the soft
X-ray spectra in NGC 5940 is largely featureless, consistent with a simple
powerlaw. Following Matzeu et al. (2020), we systematically scanned blindly
the $6.5-29\,\mathrm{\textup{\AA}}$ wavelength band of the RGS spectra for
absorption and emission line features. To do this, we included a Gaussian line
profile in the absorbed powerlaw model, allowing for both positive and
negative normalisations. The search did not result in the detection of any
prominent features that is significant.

Fig. 3 : The upper panel shows a powerlaw fit (in black) to the first order
RGS1+2 spectrum (in red) of NGC 5940. The lower panel shows the residual plot.

###  3.4 Swift XRT spectral analysis

We combined Swift XRT spectra from four observations (ObsIDs: 00037386002,
00037386003, 00037386004, 00037386005) of NGC 5940 carried out between 2010
and 2011 for improved statistics. We grouped the data to have a minimum of 5
counts per bin for use with XSPEC $\chi^{2}$ statistic. We started by fitting
the $1.5-7.0\,\mathrm{keV}$ spectra with a simple powerlaw. When extrapolated
down to $0.3\,\mathrm{keV}$, the presence of strong soft X-ray excess was
confirmed to be evident in the spectra as shown in Fig. 4. However, both the
powerlaw photon index $\Gamma$ and tbabs column density $N_{H}$ could not be
constrained due to the poorer data quality when compared to XMM-Newton.

Fig. 4 : The upper panel shows a powerlaw fit (in black) to the
$1.5-7.0\,\mathrm{keV}$ Swift XRT spectrum (in red) of NGC 5940 extrapolated
down to $0.3\,\mathrm{keV}$. The lower panel shows the residual plot.

##  4 Discussion

We carried out spectral analysis on the $2012$ XMM-Newton observation of the
Seyfert 1 AGN NGC 5940. The data revealed the presence of strong soft X-ray
emission in the source below $2\,\mathrm{keV}$ in excess of the dominant X-ray
powerlaw. Here, we investigated the possible origin of this emission
component.

We started by fitting a phenomenological model to the data, namely; the
multicolour disc blackbody, to model the soft X-ray excess with the addition
of the powerlaw model for the X-ray emission above $\sim 2\,\mathrm{keV}$. As
is the case with most fits using the blackbody model, the inferred temperature
of $\sim 0.18\,\mathrm{keV}$ is too high for an accretion disc around a
blackhole mass of $\sim 10^{7}\,M_{\odot}$. Standard accretion disc theory
predicts that the inner disc temperature of such a system is given by the
relation

$T_{eff}\sim 6.3\times
10^{5}\left(\dfrac{\dot{M}}{\dot{M}_{Edd}}\right)^{1/4}\left(\dfrac{M}{10^{8}M_{\odot}}\right)^{-1/4}\left(\dfrac{R}{R_{s}}\right)^{-3/4},$

where $\dot{M}/\dot{M}_{Edd}$ is the Eddington-scaled accretion rate, $M$ is
the blackhole mass, $R$ and $R_{s}$ are respectively the disc annular radius
and the Schwarzschild radius. Thus, for a Schwarzschild blackhole of $\sim
10^{7}\,M_{\odot}$, at the inner edge $T_{eff}$ should be of order
$~{}30\,\mathrm{eV}$ (e.g., Bhattacharyya et al., 2014; Adegoke et al., 2017).
We subsequently modelled the lower X-ray energy features using the ionised
partial covering model – multiplied by the powerlaw model. This model posits
that the lower energy features may include contributions from partial covering
of the X-ray source due to partially ionised material along the line of sight
(e.g., Miller et al., 2006; Reeves et al., 2008). This model provided a very
good fit to the data and has been employed to study the soft X-ray excess and
the spectral variability in PG 0844+349 (Gallo et al., 2011) and Zw 229.015
(Tripathi et al., 2019) among others with similarly good fits. However, a
compact absorber would imply the presence of a strong iron fluorescence line
which is not observed in the data. It is plausible that the geometry of the
absorber may be responsible for the non-detection of the line. Longer and
possibly joint multiwavelength observations of the source will throw more
light on this.

We then applied two more physically motivated models – the warm Comptonisation
and the blurred reflection models – that have been applied to model
successfully the soft excess in several Seyfert 1 AGN e.g., 1H 0707-495
(Fabian et al., 2009), Mrk 509 (Mehdipour et al., 2011), ESO 113-G010 (Cackett
et al., 2013), II Zw 177 (Pal et al., 2016), Ton S180 (Matzeu et al., 2020).
Both the thermal Comptonisation and the blurred reflection models provided
acceptable fit to the data implying they both explained the origin of the soft
excess.

The warm Comptonisation model (optxagnf) considers the soft X-ray excess to
originate due to repeated inverse-Compton scattering of seed disc photons in a
warm, optically thick corona. The fitted optical depth $\tau$ and temperature
$kT_{SE}$ of the warm corona were respectively $\sim 13$ and
$0.4\,\mathrm{keV}$ with a moderately high accretion rate value of about $5\%$
of Eddington. This inferred value of accretion rate supports earlier analysis
carried out on the source (Castelló-Mor et al., 2016). The model favoured a
maximally spinning blackhole with corona radius, photon index and power
fraction of the powerlaw having values $\sim 18\,R_{g}$, $\sim 1.7$ and $77\%$
respectively. The modelled $0.3-10\,\mathrm{keV}$ flux and luminosity came out
to be $8.2\times 10^{-12}\,\mathrm{ergs\,cm^{-2}\,s^{-1}}$ and $2.2\times
10^{43}\,\mathrm{ergs\,s^{-1}}$ respectively.

The blurred reflection model (relxill) on the other hand considers that strong
X-ray reflection in the inner region of the accretion disc, blurred and
smoothed due to strong general relativistic effects in this region is
responsible for the soft X-ray excess emission. Fit with the relxillcp flavour
of the model gave reflection fraction value of $\sim 0.5$ implying moderately
strong X-ray reflection in the disc innermost region. The model equally
favoured a maximally spinning blackhole having photon index $\sim 1.9$. The
relxill+nthcomp version of the model provided a very good fit to the data with
photon index value of $\sim 2$. The value of the ionisation parameter here is
much smaller at log$\epsilon=0.02$. An additional evidence for blurred inner
disc reflection origin for the soft X-ray excess is the expected presence of
broad iron $K_{\alpha}$ line at $\sim 6.4\,\mathrm{keV}$ as this line is also
believed to be produced due to strong X-ray reflection in the disc (Tanaka et
al., 1995). The line is however not detectable in this observation – to any
level of significance. Moreso, the RGS spectra with its higher resolution
should normally reveal some of the prominent line features that blend to
produce the soft excess, but the RGS data of this source appear to be largely
featureless, fitted by a simple powerlaw. Longer focused observation in the
X-rays and optical/UV bands may be required to probe this further.

##  5 Conclusion & Summary

We studied the spectral properties of the nearby Seyfert 1 AGN NGC 5940 using
predominantly its 2012 XMM-Newton observation supplemented with data from
Swift. This source has been spared of scrutiny over the past years despite its
proximity and availability of usable data. The data revealed strong soft X-ray
excess emission below $\sim 2\,\mathrm{keV}$ as is commonly seen in many
unobscured AGN. Using state-of-the-art spectral models, we went ahead to probe
the origin of the soft X-ray excess emission.

The phenomenological model applied provided acceptable fits to the XMM-Newton
EPIC-pn data of the source. From the disc blackbody plus powerlaw model, the
model-predicted temperature is much higher than expected for accretion discs
around supermassive blackholes. The model for partial covering from partially
ionised material tend to provide a very good fit to the data for the origin of
the soft excess. However, longer duration, high resolution multiwavelength
observations will be required to probe further the possibility of partial
covering in the source. Such observations will better probe additional
emission/absorption and other signatures that should accompany an ionised
partial covering scenario

Two additional physically motivated models were also applied – The thermal
Comptonisation and the blurred reflection models. We found that both models
provided acceptable fits and thus can explain the origin of the soft X-ray
excess emission well although with suggestive preference for the thermal
Comptonisation model. This is because other features that should normally
accompany strong reflection – like the broad iron $K_{\alpha}$ line as well as
emission line features in the RGS spectra – were not observed in our analysis.

Timing analysis may provide useful insight into understanding the importance
of thermal Componisation and/or blurred reflection in the accretion flow
physics of the source. For example, if cross-correlation analysis reveals that
optical/UV disc photons lead the X-rays in their variability, then thermal
Comptonisation could play a dominant role in the origin of the soft excess
(Adegoke et al., 2019). If on the contrary, X-rays lead, then reprocessing and
therefore blurred reflection may play a more dominant role (Krolik et al.,
1991; Cackett et al., 2007). It is equally possible that both thermal
Comptonisation and relativistically blurred reflection play a role in the
overall spectral energy distribution of the source such that while warm
Comptonisation explains the origin of the soft X-ray excess, reflection may
explain the line features as well as the hard X-ray excess. Simultaneous
observations of the source at higher X-ray energies with instruments like the
NuSTAR space telescope and at lower optical/UV energies will be crucial for a
more detailed probe of the spectral behaviour of this peculiar source.

Acknowledgements This work made use of data supplied by the XMM-Newton Science
Archive and the UK Swift Science Data Centre at the University of Leicester.
The anonymous referee is appreciated for useful comments that improved the
quality of the manuscript.

## Additional Statements

## Funding

None

## Data Availability

Details of all the data analysed in this study are included in the article.

## Declarations

## Conflict of interest

The author has no competing interests to declare regarding the content of this
article.

Informed consent None

## References

  * Adegoke et al. (2017) Adegoke, O., Rakshit, S., Mukhopadhyay, B.: Mon. Not. R. Astron. Soc. 466, 3951 (2017). https://doi.org/10.1093 /mnras/stw3320. 1612.06817
  * Adegoke et al. (2019) Adegoke, O., Dewangan, G.C., Pawar, P., Pal, M.: Astrophys. J. Lett. 870(2), 13 (2019). https://doi.org/10.3847 /2041-8213/aaf8ab. 1812.08056
  * Arnaud (1996) Arnaud, K.A.: In: Jacoby, G.H., Barnes, J. (eds.) Astronomical Data Analysis Software and Systems V. Astronomical Society of the Pacific Conference Series, vol. 101, p. 17 (1996)
  * Arnaud et al. (1985) Arnaud, K.A., Branduardi-Raymont, G., Culhane, J.L., Fabian, A.C., Hazard, C., McGlynn, T.A., Shafer, R.A., Tennant, A.F., Ward, M.J.: Mon. Not. R. Astron. Soc. 217, 105 (1985). https://doi.org/10.1093/mnras/217.1.105
  * Barth et al. (2013) Barth, A.J., Pancoast, A., Bennert, V.N., Brewer, B.J., Canalizo, G., Filippenko, A.V., Gates, E.L., Greene, J.E., Li, W., Malkan, M.A., Sand, D.J., Stern, D., Treu, T., Woo, J.-H., Assef, R.J., Bae, H.-J., Buehler, T., Cenko, S.B., Clubb, K.I., Cooper, M.C., Diamond-Stanic, A.M., Hönig, S.F., Joner, M.D., Laney, C.D., Lazarova, M.S., Nierenberg, A.M., Silverman, J.M., Tollerud, E.J., Walsh, J.L.: Astrophys. J. 769(2), 128 (2013). https://doi.org/10.1088/0004-637X/769/2/128. 1304.4643
  * Barth et al. (2015) Barth, A.J., Bennert, V.N., Canalizo, G., Filippenko, A.V., Gates, E.L., Greene, J.E., Li, W., Malkan, M.A., Pancoast, A., Sand, D.J., Stern, D., Treu, T., Woo, J.-H., Assef, R.J., Bae, H.-J., Brewer, B.J., Cenko, S.B., Clubb, K.I., Cooper, M.C., Diamond-Stanic, A.M., Hiner, K.D., Hönig, S.F., Hsiao, E., Kandrashoff, M.T., Lazarova, M.S., Nierenberg, A.M., Rex, J., Silverman, J.M., Tollerud, E.J., Walsh, J.L.: Astrophys. J. Suppl. Ser. 217(2), 26 (2015). https://doi.org/10.1088/0067-0049/217/2/26. 1503.01146
  * Bentz and Katz (2015) Bentz, M.C., Katz, S.: Publ. Astron. Soc. Pac. 127(947), 67 (2015). https://doi.org/10.1086/679601. 1411.2596
  * Bhattacharyya et al. (2014) Bhattacharyya, S., Bhatt, H., Bhatt, N., Singh, K.K.: Mon. Not. R. Astron. Soc. 440(1), 106 (2014). https://doi.org /10.1093/mnras/stu239. 1301.1150
  * Burrows et al. (2005) Burrows, D.N., Hill, J.E., Nousek, J.A., Kennea, J.A., Wells, A., Osborne, J.P., Abbey, A.F., Beardmore, A., Mukerjee, K., Short, A.D.T., Chincarini, G., Campana, S., Citterio, O., Moretti, A., Pagani, C., Tagliaferri, G., Giommi, P., Capalbi, M., Tamburelli, F., Angelini, L., Cusumano, G., Bräuninger, H.W., Burkert, W., Hartner, G.D.: Space Sci. Rev. 120(3-4), 165 (2005). https://doi.org/10.1007/s11214-005-5097-2. astro-ph/0508071
  * Cackett et al. (2007) Cackett, E.M., Horne, K., Winkler, H.: Mon. Not. R. Astron. Soc. 380, 669 (2007). https://doi.org/10.1111 /j.1365-2966.2007.12098.x. 0706.1464
  * Cackett et al. (2013) Cackett, E.M., Fabian, A.C., Zogbhi, A., Kara, E., Rey- nolds, C., Uttley, P.: Astrophys. J. Lett. 764(1), 9 (2013). https://doi.org/10.1088/2041-8205/764/1/L9. 1210.7874
  * Castelló-Mor et al. (2016) Castelló-Mor, N., Netzer, H., Kaspi, S.: Mon. Not. R. Astron. Soc. 458(2), 1839 (2016). https://doi.org/10.1093 /mnras/stw445. 1601.07177
  * Crummy et al. (2005) Crummy, J., Fabian, A.C., Brandt, W.N., Boller, T.: Mon. Not. R. Astron. Soc. 361(4), 1197 (2005). https://doi.org /10.1111/j.1365-2966.2005.09258.x. astro-ph/0506119
  * Crummy et al. (2006) Crummy, J., Fabian, A.C., Gallo, L., Ross, R.R.: Mon. Not. R. Astron. Soc. 365(4), 1067 (2006). https://doi.org /10.1111/j.1365-2966.2005.09844.x. astro-ph/0511457
  * Dauser et al. (2014) Dauser, T., Garcia, J., Parker, M.L., Fabian, A.C., Wilms, J.: Mon. Not. R. Astron. Soc. 444, 100 (2014). https://doi .org/10.1093/mnrasl/slu125. 1408.2347
  * den Herder et al. (2001) den Herder, J.W., Brinkman, A.C., Kahn, S.M., Bran- duardi-Raymont, G., Thomsen, K., Aarts, H., Audard, M., Bixler, J.V., den Boggende, A.J., Cottam, J., Decker, T., Dubbeldam, L., Erd, C., Goulooze, H., Güdel, M., Guttridge, P., Hailey, C.J., Janabi, K.A., Kaastra, J.S., de Korte, P.A.J., van Leeuwen, B.J., Mauche, C., McCalden, A.J., Mewe, R., Naber, A., Paerels, F.B., Peterson, J.R., Rasmussen, A.P., Rees, K., Sakelliou, I., Sako, M., Spodek, J., Stern, M., Tamura, T., Tandy, J., de Vries, C.P., Welch, S., Zehnder, A.: Astron. Astrophys. 365, 7 (2001). https://doi.org/10.1051/0004-6361:20000058
  * Done et al. (2012) Done, C., Davis, S.W., Jin, C., Blaes, O., Ward, M.: Mon. Not. R. Astron. Soc. 420(3), 1848 (2012). https://doi.org/ 10.1111/j.1365-2966.2011.19779.x. 1107.5429
  * Evans et al. (2009) Evans, P.A., Beardmore, A.P., Page, K.L., Osborne, J.P., O’Brien, P.T., Willingale, R., Starling, R.L.C., Burrows, D.N., Godet, O., Vetere, L., Racusin, J., Goad, M.R., Wiersema, K., Angelini, L., Capalbi, M., Chincarini, G., Gehrels, N., Kennea, J.A., Margutti, R., Morris, D.C., Mountford, C.J., Pagani, C., Perri, M., Romano, P., Tanvir, N.: Mon. Not. R. Astron. Soc. 397(3), 1177 (2009). https://doi.org/10.1111/j.1365-2966.2009.14913.x. 0812.3662
  * Fabian et al. (2004) Fabian, A.C., Miniutti, G., Gallo, L., Boller, T., Tanaka, Y., Vaughan, S., Ross, R.R.: Mon. Not. R. Astron. Soc. 353(4), 1071 (2004). https://doi.org/10.1111/j.1365-2966.2004.08036.x. astro-ph/0405160
  * Fabian et al. (2009) Fabian, A.C., Zoghbi, A., Ross, R.R., Uttley, P., Gallo, L.C., Brandt, W.N., Blustin, A.J., Boller, T., Caballero-Garcia, M.D., Larsson, J., Miller, J.M., Miniutti, G., Ponti, G., Reis, R.C., Reynolds, C.S., Tanaka, Y., Young, A.J.: Nature 459(7246), 540 (2009). https://doi.org/10.1038/na- ture08007
  * Gallo et al. (2011) Gallo, L.C., Grupe, D., Schartel, N., Komossa, S., Miniutti, G., Fabian, A.C., Santos-Lleo, M.: Mon. Not. R. Astron. Soc. 412(1), 161 (2011). https://doi.org/10.1111/j.1365-2966.2010.17894.x. 1010.4453
  * García et al. (2014) García, J., Dauser, T., Lohfink, A., Kallman, T.R., Steiner, J.F., McClintock, J.E., Brenneman, L., Wilms, J., Eikmann, W., Reynolds, C.S., Tombesi, F.: Astrophys. J. 782(2), 76 (2014). https://doi.org/10.1088/0004-637X/782/2/76. 1312.3231
  * Gehrels et al. (2004) Gehrels, N., Chincarini, G., Giommi, P., Mason, K.O., Nousek, J.A., Wells, A.A., White, N.E., Barthelmy, S.D., Burrows, D.N., Cominsky, L.R., Hurley, K.C., Marshall, F.E., Mészáros, P., Roming, P.W.A., Angelini, L., Barbier, L.M., Belloni, T., Campana, S., Caraveo, P.A., Chester, M.M., Citterio, O., Cline, T.L., Cropper, M.S., Cummings, J.R., Dean, A.J., Feigelson, E.D., Fenimore, E.E., Frail, D.A., Fruchter, A.S., Garmire, G.P., Gendreau, K., Ghisellini, G., Greiner, J., Hill, J.E., Hunsberger, S.D., Krimm, H.A., Kulkarni, S.R., Kumar, P., Lebrun, F., Lloyd-Ronning, N.M., Markwardt, C.B., Mattson, B.J., Mushotzky, R.F., Norris, J.P., Osborne, J., Paczynski, B., Palmer, D.M., Park, H.-S., Parsons, A.M., Paul, J., Rees, M.J., Reynolds, C.S., Rhoads, J.E., Sasseen, T.P., Schaefer, B.E., Short, A.T., Smale, A.P., Smith, I.A., Stella, L., Tagliaferri, G., Takahashi, T., Tashiro, M., Townsley, L.K., Tueller, J., Turner, M.J.L., Vietri, M., Voges, W., Ward, M.J., Willingale, R., Zerbi, F.M., Zhang, W.W.: Astrophys. J. 611(2), 1005 (2004). https://doi.org/10.1086/422091. astro-ph/0405233
  * Gierliński and Done (2004) Gierliński, M., Done, C.: Mon. Not. R. Astron. Soc. 349(1), 7 (2004). https://doi.org/10.1111/j.1365-2966.2004.07687 .x. astro-ph/0312271
  * Gliozzi and Williams (2020) Gliozzi, M., Williams, J.K.: Mon. Not. R. Astron. Soc. 491(1), 532 (2020). https://doi.org/10.1093/mnras/stz30 05. 1910.12115
  * Haardt and Maraschi (1993) Haardt, F., Maraschi, L.: Astrophys. J. 413, 507 (1993). https://doi.org/10.1086/173020
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., Keller, R., Kerp, J., Lenz, D., Winkel, B., Bailin, J., Calabretta, M.R., Dedes, L., Ford, H.A., Gibson, B.K., Haud, U., Janowiecki, S., Kalberla, P.M.W., Lockman, F.J., McClure-Griffiths, N.M., Murphy, T., Nakanishi, H., Pisano, D.J., Staveley-Smith, L.: Astron. Astrophys. 594, 116 (2016). https://doi.org/10.1051/0004-6361/201629178. 1610.06175
  * Jansen et al. (2001) Jansen, F., Lumb, D., Altieri, B., Clavel, J., Ehle, M., Erd, C., Gabriel, C., Guainazzi, M., Gondoin, P., Much, R., Munoz, R., Santos, M., Schartel, N., Texier, D., Vacanti, G.: Astron. Astrophys. 365, 1 (2001). https://doi.org/10.1051/0004-6361:20000036
  * Jiang et al. (2019) Jiang, J., Fabian, A.C., Dauser, T., Gallo, L., García, J.A., Kara, E., Parker, M.L., Tomsick, J.A., Walton, D.J., Reynolds, C.S.: Mon. Not. R. Astron. Soc. 489(3), 3436 (2019). https://doi.org/10.1093/mnras/stz2326. 1908.0727 2
  * Jin et al. (2012) Jin, C., Ward, M., Done, C.: Mon. Not. R. Astron. Soc. 425(2), 907 (2012). https://doi.org/10.1111/j.1365-2966.2012.21272.x. 1205.1846
  * Krolik et al. (1991) Krolik, J.H., Horne, K., Kallman, T.R., Malkan, M.A., Edelson, R.A., Kriss, G.A.: Astrophys. J. 371, 541 (1991). https://doi.org/10.1086/169918
  * Magdziarz et al. (1998) Magdziarz, P., Blaes, O.M., Zdziarski, A.A., Johnson, W.N., Smith, D.A.: Mon. Not. R. Astron. Soc. 301(1), 179 (1998). https://doi.org/10.1046/j.1365-8711.1998.02 015.x
  * Makishima et al. (1986) Makishima, K., Maejima, Y., Mitsuda, K., Bradt, H.V., Remillard, R.A., Tuohy, I.R., Hoshi, R., Nakagawa, M.: Astrophys. J. 308, 635 (1986). https://doi.org/10.1086/16 4534
  * Mason et al. (2001) Mason, K.O., Breeveld, A., Much, R., Carter, M., Cordova, F.A., Cropper, M.S., Fordham, J., Huckle, H., Ho, C., Kawakami, H., Kennea, J., Kennedy, T., Mittaz, J., Pandel, D., Priedhorsky, W.C., Sasseen, T., Shirey, R., Smith, P., Vreux, J.-M.: Astron. Astrophys. 365, 36 (2001). https://doi.org/10.1051/0004-6361:20000044. astro-ph/0011216
  * Matzeu et al. (2020) Matzeu, G.A., Nardini, E., Parker, M.L., Reeves, J.N., Braito, V., Porquet, D., Middei, R., Kammoun, E., Lusso, E., Alston, W.N., Giustini, M., Lobban, A.P., Joyce, A.M., Igo, Z., Ebrero, J., Ballo, L., Santos-Lleó, M., Schartel, N.: Mon. Not. R. Astron. Soc. 497(2), 2352 (2020). https://doi.org/10.1093/mnras/staa2076. 2007.065 75
  * Mehdipour et al. (2011) Mehdipour, M., Branduardi-Raymont, G., Kaastra, J.S., Petrucci, P.O., Kriss, G.A., Ponti, G., Blustin, A.J., Paltani, S., Cappi, M., Detmers, R.G., Steenbrugge, K.C.: Astron. Astrophys. 534, 39 (2011). https://doi.org/10.10 51/0004-6361/201116875. 1107.0659
  * Miller et al. (2006) Miller, L., Turner, T.J., Reeves, J.N., George, I.M., Porquet, D., Nandra, K., Dovciak, M.: Astron. Astrophys. 453(1), 13 (2006). https://doi.org/10.1051/0004-6361:20065276. astro-ph/0605130
  * Miniutti et al. (2009) Miniutti, G., Ponti, G., Greene, J.E., Ho, L.C., Fabian, A.C., Iwasawa, K.: Mon. Not. R. Astron. Soc. 394(1), 443 (2009). https://doi.org/10.1111/j.1365-2966.2008.14334.x. 0812.1652
  * Mitsuda et al. (1984) Mitsuda, K., Inoue, H., Koyama, K., Makishima, K., Matsuoka, M., Ogawara, Y., Shibazaki, N., Suzuki, K., Tanaka, Y., Hirano, T.: Publ. Astron. Soc. Jpn. 36, 741 (1984)
  * Novikov and Thorne (1973) Novikov, I.D., Thorne, K.S.: In: Black Holes (Les Astres Occlus), p. 343 (1973)
  * Pal et al. (2016) Pal, M., Dewangan, G.C., Misra, R., Pawar, P.K.: Mon. Not. R. Astron. Soc. 457, 875 (2016). https://doi.org/10. 1093/mnras/stw009. 1601.00503
  * Petrucci et al. (2018) Petrucci, P.-O., Ursini, F., De Rosa, A., Bianchi, S., Cappi, M., Matt, G., Dadina, M., Malzac, J.: Astron. Astrophys. 611, 59 (2018). https://doi.org/10.1051/0 004-6361/201731580. 1710.04940
  * Piconcelli et al. (2005) Piconcelli, E., Jimenez-Bailón, E., Guainazzi, M., Schartel, N., Rodríguez-Pascual, P.M., Santos-Lleó, M.: Astron. Astrophys. 432(1), 15 (2005). https://doi.org/10.1051 /0004-6361:20041621. astro-ph/0411051
  * Pounds et al. (1986) Pounds, K.A., Warwick, R.S., Culhane, J.L., de Korte, P.A.J.: Mon. Not. R. Astron. Soc. 218, 685 (1986). https://doi.org/10.1093/mnras/218.4.685
  * Reeves et al. (2008) Reeves, J., Done, C., Pounds, K., Terashima, Y., Hayashida, K., Anabuki, N., Uchino, M., Turner, M.: Mon. Not. R. Astron. Soc. 385(1), 108 (2008). https://doi.org/10.1111 /j.1745-3933.2008.00443.x. 0801.1587
  * Robinson et al. (2019) Robinson, J.H., Bentz, M.C., Johnson, M.C., Courtois, H.M., Ou-Yang, B.: Astrophys. J. 880(2), 68 (2019). https://doi.org/10.3847/1538-4357/ab29f9. 1906.07062
  * Schurch and Done (2007) Schurch, N.J., Done, C.: Mon. Not. R. Astron. Soc. 381(4), 1413 (2007). https://doi.org/10.1111/j.1365-2966.2007.12 336.x. 0706.1885
  * Shakura and Sunyaev (1973) Shakura, N.I., Sunyaev, R.A.: Astron. Astrophys. 500, 33 (1973)
  * Singh et al. (1985) Singh, K.P., Garmire, G.P., Nousek, J.: Astrophys. J. 297, 633 (1985). https://doi.org/10.1086/163560
  * Strüder et al. (2001) Strüder, L., Briel, U., Dennerl, K., Hartmann, R., Kendzio- rra, E., Meidinger, N., Pfeffermann, E., Reppin, C., Aschenbach, B., Bornemann, W., Bräuninger, H., Burkert, W., Elender, M., Freyberg, M., Haberl, F., Hartner, G., Heuschmann, F., Hippmann, H., Kastelic, E., Kemmer, S., Kettenring, G., Kink, W., Krause, N., Müller, S., Oppitz, A., Pietsch, W., Popp, M., Predehl, P., Read, A., Stephan, K.H., Stötter, D., Trümper, J., Holl, P., Kemmer, J., Soltau, H., Stötter, R., Weber, U., Weichert, U., von Zanthier, C., Carathanassis, D., Lutz, G., Richter, R.H., Solc, P., Böttcher, H., Kuster, M., Staubert, R., Abbey, A., Holland, A., Turner, M., Balasini, M., Bignami, G.F., La Palombara, N., Villa, G., Buttler, W., Gianini, F., Lainé, R., Lumb, D., Dhez, P.: Astron. Astrophys. 365, 18 (2001). https://doi.org/10.1051/0004-6361:20000066
  * Tanaka et al. (1995) Tanaka, Y., Nandra, K., Fabian, A.C., Inoue, H., Otani, C., Dotani, T., Hayashida, K., Iwasawa, K., Kii, T., Kunieda, H., Makino, F., Matsuoka, M.: Nature 375(6533), 659 (1995). https://doi.org/10.1038/375659a0
  * Tripathi et al. (2019) Tripathi, S., Waddell, S.G.H., Gallo, L.C., Welsh, W.F., Chiang, C.-Y.: Mon. Not. R. Astron. Soc. 488(4), 4831 (2019). https://doi.org/10.1093/mnras/stz1988. 1907.07048
  * Turner et al. (2001) Turner, M.J.L., Abbey, A., Arnaud, M., Balasini, M., Barbera, M., Belsole, E., Bennie, P.J., Bernard, J.P., Bignami, G.F., Boer, M., Briel, U., Butler, I., Cara, C., Chabaud, C., Cole, R., Collura, A., Conte, M., Cros, A., Denby, M., Dhez, P., Di Coco, G., Dowson, J., Ferrando, P., Ghizzardi, S., Gianotti, F., Goodall, C.V., Gretton, L., Griffiths, R.G., Hainaut, O., Hochedez, J.F., Holland, A.D., Jourdain, E., Kendziorra, E., Lagostina, A., Laine, R., La Palombara, N., Lortholary, M., Lumb, D., Marty, P., Molendi, S., Pigot, C., Poindron, E., Pounds, K.A., Reeves, J.N., Reppin, C., Rothenflug, R., Salvetat, P., Sauvageot, J.L., Schmitt, D., Sembay, S., Short, A.D.T., Spragg, J., Stephen, J., Strüder, L., Tiengo, A., Trifoglio, M., Trümper, J., Vercellone, S., Vigroux, L., Villa, G., Ward, M.J., Whitehead, S., Zonca, E.: Astron. Astrophys. 365, 27 (2001). https://doi.org/10.1051/0004-6361:20000087. astro-ph/0011498
  * Walter and Fink (1993) Walter, R., Fink, H.H.: Astron. Astrophys. 274, 105 (1993)
  * Walton et al. (2013) Walton, D.J., Nardini, E., Fabian, A.C., Gallo, L.C., Reis, R.C.: Mon. Not. R. Astron. Soc. 428(4), 2901 (2013). https://doi.org/10.1093/mnras/sts227. 1210.4593
  * Wilms et al. (2000) Wilms, J., Allen, A., McCray, R.: Astrophys. J. 542(2), 914 (2000). https://doi.org/10.1086/317016. astro-ph/0008425
  * Zdziarski et al. (1996) Zdziarski, A.A., Johnson, W.N., Magdziarz, P.: Mon. Not. R. Astron. Soc. 283(1), 193 (1996). https://doi.org/10.10 93/mnras/283.1.193. astro-ph/9607015
  * Życki et al. (1999) Życki, P.T., Done, C., Smith, D.A.: Mon. Not. R. Astron. Soc. 309(3), 561 (1999). https://doi.org/10.1046/j.1365-8711.1999.02885.x. astro-ph/9904304