# A High Order Cartesian Grid, Finite Volume Method for Elliptic Interface
Problems

Will Thacher  Hans Johansen  Daniel Martin

###### Abstract

We present a higher-order finite volume method for solving elliptic PDEs with
jump conditions on interfaces embedded in a 2D Cartesian grid. Second, fourth,
and sixth order accuracy is demonstrated on a variety of tests including
problems with high-contrast and spatially varying coefficients, large
discontinuities in the source term, and complex interface geometries. We
include a generalized truncation error analysis based on cell-centered Taylor
series expansions, which then define stencils in terms of local discrete
solution data and geometric information. In the process, we develop a simple
method based on Green’s theorem for computing exact geometric moments directly
from an implicit function definition of the embedded interface. This approach
produces stencils with a simple bilinear representation, where spatially-
varying coefficients and jump conditions can be easily included and finite
volume conservation can be enforced.

###### keywords:

Elliptic Interface Problem , High Order , Embedded Boundary , Cut Cell ,
Finite Volume , Jump conditions , Discontinuous coefficients , Variable
coefficients

††journal: Journal of Computational Physics

[inst1]organization= Graduate Student Researcher, Applied Science and
Technology Group, University of California Berkeley, city=Berkeley,state=CA,
postcode=94720, country=United States

[inst2]organization= Applied Numerical Algorithms Group, Lawrence Berkeley
National Laboratory, city=Berkeley,state=CA, postcode=94720, country=United
States

##  1 Introduction

Elliptic PDEs with discontinuities in the source term, coefficients and
solution form an important class of equations in computational science and
engineering. These equations arise from mathematical models of multi-material
systems, multi-phase flows, crystal growth, and many other physical processes
[1]. Solving such equations numerically is not straightforward because the
accuracy of the scheme is typically based on smoothness assumptions that do
not in general apply at the interface.

Numerous schemes have been proposed to solve this problem based on finite
difference, finite volume, and finite element formulations. These methods can
roughly be classified into those that treat the interface explicitly by
creating elements that conform to the shape of the interface, or those that
represent the interface implicitly by “embedding” it onto a non-conforming
mesh (see [2] for a thorough review and further references). In the finite
element realm, methods such as [3] body-fit the mesh to the interface whereas
methods such as [4] use a fixed mesh and modify basis functions where the
interface crosses elements. A widely used and influential method in the finite
difference category is the Immersed Interface Method (IIM) [5]. The IIM uses
standard Cartesian grid finite difference stencils away from the interface and
modifies stencils near the interface using one-sided Taylor series expansions
that incorporate jump conditions. The Ghost Fluid Method [6] extrapolates the
solution across the interface to nearby grid points by incorporating jump
conditions so that standard stencils can still be used at all grid points.
These various finite difference methods are closely related to various schemes
for imposing boundary conditions; jump conditions can be thought of as a sort
of implicit boundary condition that depends on the solution itself.

This paper is concerned with the third category: finite volume schemes. These
methods are conservative (in the sense that the divergence theorem is applied
to a control volume), and have well-studied stability properties [7]. The
embedded boundary (EB) method of [8] combines the implicit and explicit
interface representations: the interface is embedded onto a Cartesian grid,
forming cut-cell volumes of arbitrary shapes where it intersects rectangular
cells. The elliptic equation is then discretized in flux divergence form using
techniques developed in [7] and [9] with appropriate modifications made at the
interface to enforce jump conditions. The method developed in [8] is second-
order accurate $L^{\infty}$ and $L^{1}$ norm, but is difficult to extend to
higher order accuracy. For low-order methods, cell averages can be treated as
point values to second-order accuracy, so finite difference type schemes can
be employed to create stencils. This is not the case for higher order finite
volume methods; integration must be performed over arbitrarily-shaped “cut
cells.” Techniques such as choosing a midpoint or centroid as a quadrature
rule for surface integrals of the flux is not sufficient for higher-order
accuracy.

Many of these difficulties are being addressed by recent developments in
higher-order finite volume and EB methods, which are summarized in [10]. One
example is the use of weighted least-squares interpolation for stencil
construction in complex geometries ([11], [12]) as well as the derivation of
high-order stencils on Cartesian grids [13]. Given the close relationship
between boundary conditions and jump conditions, we propose extending the
methodology of [12] to the elliptic interface problem. The primary
contribution of this research is a finite volume method for the variable
coefficient 2D elliptic interface problem that is 1) high order accurate and
2) conservative. In the process, we have also created approaches for 1) an
efficient technique for generating exact geometric information from an
implicit function and 2) a method for building high-order finite volume
stencils for variable coefficient elliptic operators on arbitrary cut-cell
meshes.

The outline of the paper is as follows: In Section 2, we define mesh and
geometric quantities and give a general truncation error analysis which allows
us to design stencils of arbitrarily high order. In Section 3, we describe in
detail our method for constructing stencils. In Section 4, we present results
that validate the approach using a series of model problems that test
different aspects of the scheme.

##  2 Discretization

$V_{+,\mathbf{i}}$$V_{-,\mathbf{i}}$$A_{B,\mathbf{i}}$$\hat{\mathbf{n}}$$A_{+,\mathbf{i}-\frac{1}{2}\mathbf{e}_{y}}$$A_{+,\mathbf{i}+\frac{1}{2}\mathbf{e}_{y}}$$A_{-,\mathbf{i}+\frac{1}{2}\mathbf{e}_{y}}$$A_{+,\mathbf{i}-\frac{1}{2}\mathbf{e}_{x}}$$A_{+,\mathbf{i}+\frac{1}{2}\mathbf{e}_{x}}$$A_{-,\mathbf{i}+\frac{1}{2}\mathbf{e}_{x}}$
Figure 1: Cut cell geometric quantities that make up the finite volume
notation.

Let $\Omega$ be a physical domain that is divided into the subdomains
$\Omega^{+}$ and $\Omega^{-}$ by an interface $\Gamma$. We consider the
variable coefficient elliptic interface problem for $u(x)$:

$\displaystyle\alpha u-\nabla\cdot\left(\beta\nabla u\right)$
$\displaystyle=f\hbox{ on }\Omega$ (1) $\displaystyle\left[u\right]$
$\displaystyle=w\hbox{ on }\Gamma$ (2)
$\displaystyle\left[\beta\partial_{\mathbf{n}}u\right]$ $\displaystyle=v\hbox{
on }\Gamma\,.$ (3)

Here, $\left[\cdot\right]$ denotes a jump in some quantity at the interface:
$\left[u\right]=u^{+}(\mathbf{x})-u^{-}(\mathbf{x})$ at some point
$\mathbf{x}\in\Gamma$, and the term $\partial_{\mathbf{n}}u\equiv\nabla
u\cdot\hat{\mathbf{n}}$ represents a flux at this boundary with unit normal
$\hat{\mathbf{n}}(\mathbf{x})$. Finally, coefficients
$\alpha^{\pm}(\mathbf{x}),\beta^{\pm}(\mathbf{x})$ and the source term
$f^{\pm}(\mathbf{x})$ vary in space and may be discontinuous across $\Gamma$.

The domain $\Omega$ is discretized into a Cartesian mesh of square control
volumes (or “cells”) $V_{p,\mathbf{i}}$, $\mathbf{i}\in\mathbb{Z}^{2}$, that
have centroids $\mathbf{x}_{p,\mathbf{i}}$ and side lengths of scale $h$, the
grid spacing (see Figure 1). We indicate $p\in\left\\{+,-\right\\}$ to specify
a subdomain of $\Omega_{p}$; $p$ can often be thought of as the phase or
material type of a physical quantity. We assume that each cell
$V_{p,\mathbf{i}}$ may have up to four grid-aligned faces, which we label
$A_{p,\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{d}}$, where $\mathbf{e}_{d}$ is the
unit vector in direction $d$.

Any cell that is intersected by $\Gamma$, the “embedded boundary” (EB), is
called a “cut” cell. We make the following assumptions to simplify the
geometric considerations. First, a cut cell consists of only two control
volumes $V_{+,\mathbf{i}}$ and $V_{-,\mathbf{i}}$ divided by a portion of the
EB, denoted by $A_{B,\mathbf{i}}$. The unit normal vector $\hat{\mathbf{n}}$
on $\Gamma$ points from $\Omega^{+}$ to $\Omega^{-}$. So, along with
$A_{p,\mathbf{i}\pm\frac{1}{2}\mathbf{e}_{d}}$ as the grid-aligned faces of
each portion of the cut cell, each cut cell must have a total of at least 3,
and at most 5, faces.

Because we are using a finite volume formulation, we should define additional
geometric quantities that will be useful throughout this paper: a geometric
“moment” is an integral of a centered monomial over some specified region. We
define four moments corresponding to four components of the geometry:

$\displaystyle m_{{p,\mathbf{i}}}^{\mathbf{q}}$
$\displaystyle=\int_{V_{p,\mathbf{i}}}\mathbf{(x-\bar{x})^{q}}\ dV$ (4)
$\displaystyle m_{{p,\mathbf{i}\pm\frac{1}{2}e_{d}}}^{\mathbf{q}}$
$\displaystyle=\int_{A_{p,\mathbf{i}\pm\frac{1}{2}e_{d}}}\mathbf{(x-\bar{x})^{q}}\
dA$ (5) $\displaystyle m_{{B,\mathbf{i}}}^{\mathbf{q}}$
$\displaystyle=\int_{A_{B,\mathbf{i}}}\mathbf{(x-\bar{x})^{q}}\ dA$ (6)
$\displaystyle m_{{B,\mathbf{i},d}}^{\mathbf{q}}$
$\displaystyle=\int_{A_{B,\mathbf{i},d}}\mathbf{(x-\bar{x})^{q}}\ \hat{n}_{d}\
dA\ ,$ (7)

where $\mathbf{q}=\left[q_{x}\ q_{y}\right]$ is a vector of non-negative
integers, and we use the multi-index notation
$\mathbf{(x-\bar{x})^{q}}=(x-\bar{x})^{q_{x}}(y-\bar{y})^{q_{y}}$. The multi-
indices have sum of at most $P$, $|\mathbf{q}|\leq P$, and are ordered
lexicographically: $\left\\{00,10,20,...P0,01,11,...0P\right\\}$. This allows
us to refer to $\mathbf{v}\left[\mathbf{q}\right]$ as the $\mathbf{q}^{th}$
entry of a vector $\mathbf{v}$. Thus by definition, the volume of cell
$\mathbf{i}$ is $\left|V_{p,\mathbf{i}}\right|=m_{{p,\mathbf{i}}}^{00}$, and
the centroid $\bar{\mathbf{x}}_{p,\mathbf{i}}$ of $V_{p,\mathbf{i}}$ is
$\frac{1}{|V_{p,\mathbf{i}}|}\left[m_{{p,\mathbf{i}}}^{10},m_{{p,\mathbf{i}}}^{01}\right]$.
Similarly, $m_{{B,\mathbf{i}}}^{00}$ is the area of the EB, and
$m_{{B,\mathbf{i},x}}^{00}$ is its $x$ normal component-weighted area, or $x$
direction _cross-section_. For ease of notation, throughout this paper we
ignore $\bar{\mathbf{x}}$, although in practice it is the cell-center of each
full Cartesian cell. Finally, our method for calculating these moments in 2D
is exact to within roundoff errors, and is detailed in the appendix.

Two types of variables are stored on the mesh: cell-averaged quantities
$\left\langle
u\right\rangle_{p,\mathbf{i}}=\frac{1}{|V_{p,\mathbf{i}}|}\int_{V_{p,\mathbf{i}}}u\
dV$, and centroid-centered quantities
$u_{p,\mathbf{i}}=u(\mathbf{x}_{p,\mathbf{i}})$. The coefficients $\alpha$ and
$\beta$ are given as point values at the centroids of cells. The right-hand-
side function is provided as cell-averaged values, $\langle f\rangle$, of
sufficient accuracy, and we solve for cell-averaged values $\langle u\rangle$.

Within this context, our finite volume scheme solves the discrete system:

$\displaystyle\left\langle\alpha
u\right\rangle_{p,\mathbf{i}}-\left\langle\nabla\cdot\beta\nabla
u\right\rangle_{p,\mathbf{i}}$ $\displaystyle=\left\langle
f\right\rangle_{p,\mathbf{i}}$ (8) $\displaystyle\left[u\right]_{\mathbf{i}}$
$\displaystyle=w_{\mathbf{i}}$ (9)
$\displaystyle\left[\beta\partial_{\mathbf{n}}u\right]_{\mathbf{i}}$
$\displaystyle=v_{\mathbf{i}}\ .$ (10)

for $\left\langle u\right\rangle_{p,\mathbf{i}}$ in each volume
$V_{p,\mathbf{i}}$ in the mesh, subject to problem-specific boundary
conditions. The system of equations we solve will have one degree of freedom
in full cells and two degrees of freedom in cut cells.
$\left[\cdot\right]_{\mathbf{i}}$ denotes the integral of the jump of a
quantity across the EB in cut cell $\mathbf{i}$:

$\displaystyle\left[u\right]_{\mathbf{i}}=\int_{A_{B,\mathbf{i}}}u^{+}-u^{-}dA\
.$ (11)

The objective of the following section is to provide a general truncation
error analysis that will allow us to discretize (8) – (10) to high order
accuracy.

###  2.1 Error Analysis

In the finite volume or finite difference context, a _stencil_ approximates
some functional $G(u)$ by a linear combination of local information about the
function $u$. This functional is typically point values or integrals over some
region of derivative of $u$. The function information, or data, can include
boundary conditions, jump conditions, point values of the function, or
averages of the function over some nearby region. Let $\mathbf{d}$ denote this
vector of local function data, and $\mathbf{s}$ be the vector of stencil
values corresponding to each of the data points in $\mathbf{d}$. The
truncation error $\tau$ of this stencil is defined as:

$\displaystyle\tau=\mathbf{s}^{T}\mathbf{d}-G(u)\ .$ (12)

For the present problem, $G$ will be an integral of $u$, or some combination
of its partial derivatives, over a one or two dimensional region. If we can
approximate $u$ using a truncated Taylor series, $G(u)$ can be approximated
with a linear combination of Taylor series coefficients of $u$ up to the
desired order of accuracy.

Throughout this section we drop the subscripts $p$ and $\mathbf{i}$ except
where necessary for clarity. We can express $u$ locally as a Taylor series
expansion and remainder term:

$\displaystyle u(\mathbf{x})=\sum_{|\mathbf{q}|\leq
P}\frac{1}{\mathbf{q!}}u^{(\mathbf{q})}(\mathbf{\bar{x}})\mathbf{x}^{\mathbf{q}}+O\left(h^{P+1}\right),$
(13)

where, again using multi-index notation,
$\mathbf{x}^{\mathbf{q}}=x^{q_{x}}y^{q_{y}}$, $\mathbf{q!}=q_{x}!q_{y}!$, and
$u^{(\mathbf{q})}=\frac{\partial^{q_{x}}\partial^{q_{y}}}{\partial
x^{q_{x}}\partial y^{q_{y}}}u$. The Taylor polynomial is then just

$\displaystyle u(\mathbf{x})$ $\displaystyle=\sum_{|\mathbf{q}|\leq
P}c^{\mathbf{q}}_{u}\mathbf{x}^{\mathbf{q}}+O\left(h^{P+1}\right)\hbox{, where
}$ (14) $\displaystyle c^{\mathbf{q}}_{u}$
$\displaystyle=\frac{1}{\mathbf{q!}}u^{(\mathbf{q})}(\mathbf{\bar{x}})\,.$

Integrating the flux divergence term in (8) over the discrete volume $V$ and
applying Gauss’ theorem we obtain:

$\displaystyle\int_{V}\nabla\cdot\beta\nabla u\ dV$
$\displaystyle=\left[\sum_{\pm,\ d}(\pm
1)\int_{A_{\pm\frac{1}{2}e_{d}}}\beta\frac{\partial u}{\partial x_{d}}\
dA\right]+\int_{A_{B}}\beta\nabla u\cdot\hat{\mathbf{n}}\ dA\ .$ (15)

Expressing $\beta$ and $u$ as Taylor expansions, for a surface integral over
any (EB or grid-aligned) face $A$ we have:

$\displaystyle\int_{A}\beta\nabla u\cdot\hat{\mathbf{n}}\ dA$
$\displaystyle=\int_{A}\left(\sum_{|\mathbf{r}|\leq
P}c_{\beta}^{\mathbf{r}}\mathbf{x}^{\mathbf{r}}\right)\left(\sum_{|\mathbf{q}|\leq
P}c_{u}^{\mathbf{q}}\left[\frac{\partial\mathbf{x}^{\mathbf{q}}}{\partial
x}\hat{n}_{x}+\frac{\partial\mathbf{x}^{\mathbf{q}}}{\partial
y}\hat{n}_{y}\right]\right)+O(h^{P})\ dA$ (16)
$\displaystyle=\int_{A}\left(\sum_{|\mathbf{r}|\leq
P}c_{\beta}^{\mathbf{r}}\mathbf{x}^{\mathbf{r}}\right)\left(\sum_{|\mathbf{q}|\leq
P}c_{u}^{\mathbf{q}}\left[q_{x}\mathbf{x}^{\mathbf{q}-\mathbf{e}_{x}}\hat{n}_{x}+q_{y}\mathbf{x}^{\mathbf{q}-\mathbf{e}_{y}}\hat{n}_{y}\right]\right)+O(h^{P})\
dA$ (17) $\displaystyle=\sum_{|\mathbf{q}|\leq
P}\left[\sum_{|\mathbf{r}|\leq(P-|\mathbf{q}|)}c_{\beta}^{\mathbf{r}}\left(q_{x}m_{A,x}^{\mathbf{q}+\mathbf{r}-\mathbf{e}_{x}}+q_{y}m_{A,y}^{\mathbf{q}+\mathbf{r}-\mathbf{e}_{y}}\right)\right]c_{u}^{\mathbf{q}}+O(h^{P+1})\
.$ (18)

In order to compute $\left\langle\nabla\cdot\beta\nabla u\right\rangle$, we
apply (18) to each surface integral in (15) and divide by the cell volume.
Since $\left|V\right|$ is an $O(h^{2})$ quantity, the error term for the cell-
averaged flux divergence term is $O(h^{P-1})$.

Similarly, for the linear term in (8), we have:

$\displaystyle\left\langle\alpha u\right\rangle$
$\displaystyle=\frac{1}{\left|V\right|}\int_{V}\sum_{|\mathbf{r}|\leq
P}c_{\alpha}^{\mathbf{r}}\mathbf{x}^{\mathbf{r}}\sum_{|\mathbf{q}|\leq
P}c_{u}^{\mathbf{q}}\mathbf{x}^{\mathbf{q}}+O(h^{P+1})dV$ (19)
$\displaystyle=\frac{1}{\left|V\right|}\sum_{|\mathbf{q}|\leq
P}\left[\sum_{|\mathbf{r}|\leq(P-|\mathbf{q}|-2)}c_{\alpha}^{\mathbf{r}}m^{\mathbf{r}+\mathbf{q}}\right]c_{u}^{\mathbf{q}}+O(h^{P+1})\
.$ (20)

We see that our approximation to the functionals $G$ of interest can be
written in the general form:

$\displaystyle G(u)=\sum_{\left|\mathbf{q}\right|\leq
P}g^{\mathbf{q}}c_{u}^{\mathbf{q}}+O(h^{R})=\mathbf{g}^{T}\mathbf{c}_{u}+O(h^{R})\
,$ (21)

where $\mathbf{c}_{u}$ is the vector of approximate Taylor series coefficients
for $u$, and $\mathbf{g}$ is the vector of terms that result from operating on
these coefficients, shown in brackets in (18) and (20). The order $R$ error
term for the functional approximation may be different than the order $P+1$
error term for the Taylor series of $u$ as a result of differentiation and
integration.

To calculate the truncation error $\tau$ in (12), we must fill in the vector
$\mathbf{d}$ with function data. Suppose, for example, that the function data
are cell-averaged values of $u$. Then we can write:

$\displaystyle\left\langle u\right\rangle_{\mathbf{j}}$
$\displaystyle=\frac{1}{|V_{\mathbf{j}}|}\int_{V_{\mathbf{j}}}\sum_{|\mathbf{q}|\leq
P}c^{\mathbf{q}}_{u}\mathbf{x}^{\mathbf{q}}+O\left(h^{P+1}\right)dV$ (22)
$\displaystyle=\sum_{|\mathbf{q}|\leq
P}\frac{m_{\mathbf{j}}^{\mathbf{q}}}{|V_{\mathbf{j}}|}c^{\mathbf{q}}_{u}+O\left(h^{P+1}\right)=\mathbf{m}_{\mathbf{j}}^{T}\mathbf{c}_{u}+O\left(h^{P+1}\right)\
,$ (23)

where $\mathbf{m}_{\mathbf{j}}$ is the vector of cell-averaged volume moments
for any cell $\mathbf{j}$. If we have $n$ such cell-averaged values, then we
can write:

$\displaystyle\mathbf{d}=\begin{bmatrix}\mathbf{m}_{\mathbf{j}_{1}}^{T}\\\
\mathbf{m}_{\mathbf{j}_{2}}^{T}\\\ ...\\\
\mathbf{m}_{\mathbf{j}_{n}}^{T}\end{bmatrix}\mathbf{c}_{u}+O(h^{P+1})=\mathbf{M}\mathbf{c}_{u}+O(h^{P+1})\
,$ (24)

where $\mathbf{M}$ is the “moment matrix” whose rows are the vectors
$\mathbf{m}_{\mathbf{j}_{i}}^{T}$. Inserting this into (12) we have:

$\displaystyle\tau=\mathbf{s}^{T}\mathbf{M}\mathbf{c}_{u}+O(h^{P+1})-\left(\mathbf{g}^{T}\mathbf{c}_{u}+O(h^{R})\right)\
.$ (25)

To ensure that $\tau$ is $O(h^{\min(R,P+1)})$, we must have:

$\displaystyle\mathbf{s}^{T}\mathbf{M}\mathbf{c}_{u}$
$\displaystyle=\mathbf{g}^{T}\mathbf{c}_{u}\,.$ (26)

This must hold for any $u$ with a Taylor series represented by arbitrary
$\mathbf{c}_{u}$, so that

$\displaystyle\mathbf{M}^{T}\mathbf{s}$ $\displaystyle=\mathbf{g}\,.$ (27)

Note that this linear system is a general form for the _method of undetermined
coefficients_ , as explained in [14] (Chapter 1.2): stencil weights are chosen
so that the sum of Taylor series terms of the stencil exactly matches the
Taylor series terms of the functional up to some order. In this case, we are
working with cell averages rather than pointwise evaluations, but the
principle is the same. If the neighborhood of local function data is chosen
such that this linear system is exactly determined or undetermined, the least
norm solution to this system is given by the pseudoinverse of $\mathbf{M}$:

$\displaystyle\mathbf{s}=(\mathbf{M}^{T})^{\dagger}\mathbf{g}\,.$ (28)

If this system is undetermined, there are infinitely many stencils that will
have the same order of truncation error. However, there is another profitable
way to view this stencil construction process that justifies using the least-
norm stencil. Suppose we want to interpolate the data stored in the vector
$\mathbf{d}$ with a degree $P$ polynomial. If the matrix $\mathbf{M}$ is full
rank, we can determine the coefficients of this polynomial by solving, in a
least-squares sense, the linear system in (24):

$\displaystyle\mathbf{c}_{u}$ $\displaystyle=\mathbf{M}^{\dagger}\mathbf{d}\
.$ (29)

Inserting this into (21), we have:

$\displaystyle G(u)$
$\displaystyle=\mathbf{g}^{T}\mathbf{M}^{\dagger}\mathbf{d}+O\left(h^{\min(R,P+1)}\right)\
.$ (30)

Finally, inserting this into the truncation error expression (12):

$\displaystyle\mathbf{s}^{T}\mathbf{d}$
$\displaystyle=\mathbf{g}^{T}\mathbf{M}^{\dagger}\mathbf{d}\implies$ (31)
$\displaystyle\mathbf{s}$ $\displaystyle=(\mathbf{M}^{T})^{\dagger}\mathbf{g}\
,$ (32)

where the implication follows from the fact this must hold for any
$\mathbf{d}$.

Thus we can view our stencils as originating from either the undetermined
system $\mathbf{M}^{T}\mathbf{s}=\mathbf{g}$, in which $\mathbf{s}$ is chosen
to cancel lower order Taylor series terms, or the overdetermined system
$\mathbf{M}\mathbf{c}_{u}=\mathbf{d}$, in which we fit an interpolating
polynomial to local data. This error analysis is quite general; we have made
no mention of the shape of the volumes in the mesh, only that we know their
geometric moments to sufficient accuracy. The interface jump conditions and
boundary conditions are considered to be pieces of function data that can be
used to build stencils near an interface or boundary. Taking the
overdetermined perspective, by doing so we will constrain our interpolating
polynomials, and therefore the numerical solution, to match boundary and jump
conditions.

##  3 Stencil Construction

In the previous section we showed that if we approximate the solution $u$ with
a polynomial whose coefficients are mapped from local function data by
$\mathbf{M}^{\dagger}$, then our stencil will take the simple form (28). In
this section we will describe in detail our method for computing the moment
matrix $\mathbf{M}$ and the vector $\mathbf{g}$ of terms that result from
operating on the Taylor polynomials approximating $u$. From (18) and (20), we
see that to achieve a truncation error of order $P-1$ for the two terms in
(8), we need to “calculate” Taylor coefficients $c_{\phi}^{\mathbf{q}}$ for
$\phi\in\left\\{u,\beta,\alpha\right\\}$ up to order $P$. In this section we
again drop the subscript $p$ except where necessary for clarity.

Both the linear term and the flux divergence term can be expressed as linear
combinations of the Taylor series terms $c_{u}^{\mathbf{q}}$, as seen in (20)
and (18). For the linear term, let $\mathbf{g}_{\alpha,\mathbf{i}}$ be the
vector whose $\mathbf{q}^{th}$ entry is

$\displaystyle\mathbf{g}_{\alpha,\mathbf{i}}\left[\mathbf{q}\right]=\sum_{|\mathbf{r}|\leq(P-|\mathbf{q}|-2)}c_{\alpha}^{\mathbf{r}}m_{\mathbf{i}}^{\mathbf{r}+\mathbf{q}}\
.$ (33)

Applying (32), we have that a stencil for the linear term is given by:

$\displaystyle\mathbf{s}_{\alpha,\mathbf{i}}=(\mathbf{M}_{u,\mathbf{i}}^{T})^{\dagger}\mathbf{g}_{\alpha,\mathbf{i}}$
(34)

where the subscript $\alpha,\mathbf{i}$ indicates a stencil for the linear
term, $\left\langle\alpha u\right\rangle$, over cell $\mathbf{i}$. For the
integral of $\beta\nabla u\cdot\hat{\mathbf{n}}$ over a face $A$, we let
$\mathbf{g}_{\beta,A}$ be the vector whose $\mathbf{q}^{th}$ entry is

$\displaystyle\mathbf{g}_{\beta,A}\left[\mathbf{q}\right]=\sum_{|\mathbf{r}|\leq(P-|\mathbf{q}|)}c_{\beta}^{\mathbf{r}}\left(q_{x}m_{A,x}^{\mathbf{q}+\mathbf{r}-\mathbf{e}_{x}}+q_{y}m_{A,y}^{\mathbf{q}+\mathbf{r}-\mathbf{e}_{y}}\right)\
.$ (35)

A stencil for the flux surface integral (16) is then given by:

$\displaystyle\mathbf{s}_{\beta,A}=(\mathbf{M}_{u,\mathbf{i}}^{T})^{\dagger}\mathbf{g}_{\beta,A}\
.$ (36)

Each entry in $\mathbf{g}_{\alpha,\mathbf{i}}$ involves the Taylor
coefficients of the $\alpha$ field, and likewise each $\mathbf{g}_{\beta,A}$
entry of involves the Taylor coefficients of the $\beta$ field. We must
generate these coefficients for each volume $V$, which can be done with the
same moment matrix formulation that we use to compute the Taylor coefficients
$\mathbf{c}_{u}$.

Recall that $\alpha$ is given as point values at the centroids of cells. Let
$\mathcal{N}_{p,\mathbf{i}}$ be some neighborhood of cells
$\mathbf{j}_{k}\in\left\\{\mathbf{j}_{1},\dots,\mathbf{j}_{n}\right\\}$ in
phase $p$ around volume $V_{p,\mathbf{i}}$. A cell $V_{p,\mathbf{j}}$ in this
neighborhood has centroid $\mathbf{x}_{p,\mathbf{j}}$, with
$\alpha_{p,\mathbf{j}}=\alpha(\mathbf{x}_{p,\mathbf{j}})$. Our function data
will consist of $\alpha_{p,\mathbf{j}}$ at each point $\mathbf{j}_{k}$ in the
neighborhood, which we compile into the vector
$\mathbf{d}_{\alpha,p,\mathbf{i}}$. Let $\mathbf{M}_{\alpha,p,\mathbf{i}}$ be
the variable coefficient interpolation matrix whose rows consist of monomials
in the Taylor expansion of $\alpha$ evaluated at $\mathbf{x}_{p,\mathbf{j}}$:

$\displaystyle\mathbf{M}_{\alpha,p,\mathbf{i}}\left[k,\mathbf{q}\right]=\mathbf{x}_{p,\mathbf{j}_{k}}^{\mathbf{q}}\
.$ (37)

Letting $\mathbf{c}_{\alpha,p}$ be the vector of Taylor series coefficients,
we have

$\displaystyle\mathbf{c}_{\alpha,p}=\mathbf{M}_{\alpha,p,\mathbf{i}}^{\dagger}\mathbf{d}_{\alpha,p,\mathbf{i}}\
.$ (38)

We can do the same for $\beta$ to calculate $\mathbf{c}_{\beta,p}$.

Therefore we have that:

$\displaystyle\mathbf{g}_{\alpha,\mathbf{i}}$
$\displaystyle=\mathbf{G}_{\alpha,\mathbf{i}}\mathbf{M}_{\alpha,\mathbf{i}}^{\dagger}\mathbf{d}_{\alpha,\mathbf{i}}\
,$ (39)

where $\mathbf{G}_{\alpha,\mathbf{i}}$ is the matrix whose
$\mathbf{q},\mathbf{r}$ entry is $m_{\mathbf{i}}^{\mathbf{r}+\mathbf{q}}$.
Combining this with (34), we can finally write:

$\displaystyle\mathbf{s}_{\alpha,\mathbf{i}}^{T}$
$\displaystyle=\mathbf{d}_{\alpha,\mathbf{i}}^{T}(\mathbf{M}_{\alpha,\mathbf{i}}^{\dagger})^{T}\mathbf{G}_{\alpha,\mathbf{i}}^{T}\mathbf{M}_{u,\mathbf{i}}^{\dagger}\
.$ (40)

This stencil is bilinear in $\alpha$ and $u$, which is appropriate for the
bilinear functional we are trying to approximate. We can obtain a similar
formulation for the flux stencil $\mathbf{s}_{\beta,A}^{T}$ for each face $A$
of the cell. Thus constructing stencils is just a matter of constructing the
moment matrices for $u,\beta$ and $\alpha$. We now describe how to construct
moment matrices at, near, and away from the interface.

Figure 2: For the $P=2$ scheme, the regular cell footprint is a standard five-
point Laplacian, and if any point in the footprint is a cut cell, it is then
“irregular.” Cut cells are shown with dark shading, irregular cells with light
shading, and the remaining white cells are “regular.”

###  3.1 Moment Matrices

We partition our cells into three subsets: cut cells $\Omega_{C}$, irregular
cells $\Omega_{I}$, and regular cells $\Omega_{R}$. Cut cells are intersected
by the EB. Irregular cells are not intersected by the EB, but at least one
cell in the stencil footprint for a regular cell is intersected by the EB. See
Figure 2. The method for constructing the moment matrices is different for
each of these three types of cells.

####  3.1.1 Regular Cells

x (a)

$x^{6}$$x^{5}y$$x^{4}y^{2}$$x^{3}y^{3}$$x^{2}y^{4}$$xy^{5}$$y^{6}$$x^{5}$$x^{4}y$$x^{3}y^{2}$$x^{2}y^{3}$$xy^{4}$$y^{5}$$x^{4}$$x^{3}y$$x^{2}y^{2}$$xy^{3}$$y^{4}$$x^{3}$$x^{2}y$$xy^{2}$$y^{3}$$x^{2}$$xy$$y^{2}$$x$$y$$1$
(b)

Figure 3: Figure (a) shows increasing footprints for the moment matrix for
$P=2,4,6$ schemes with increasingly lighter shades. The footprints “support”
the corresponding (same shading) even monomial terms in Pascal’s triangle,
(b). As explained in section 3.1.1, odd moments of order $P$ are not needed to
achieve order $P$ truncation error. In addition, we require at most order
$P-1$ monomial terms for the coefficient $\beta$. The cell face
$A_{\mathbf{i}-\mathbf{e}_{x}}$ is marked with a red x, and the the dashed,
dotted, and solid blue outlines (for $P=2,4,6$, respectively), show the
footprint for the modified moment matrix
$\mathbf{M}_{\beta,A_{\mathbf{i}-\mathbf{e}_{x}}}$ in (42). The corresponding
dashed, dotted, and solid blue boxed monomials in (b) show the order $P-1$
terms that are supported by these footprints.

The vast majority of cells will be regular and will all have the same bilinear
stencil, meaning we only have to solve for this stencil once. Since regular
cells are squares, the integral of any monomial error term with odd degree
over a regular cell is $0$. This means that for the the flux divergence term,
in regular cells we can achieve a truncation error of order $P$ using an order
$P$ polynomial. We can achieve an order $P$ cell averaged linear term with an
order $P-2$ polynomial.

Furthermore, we do not need to calculate all of the Taylor polynomial
coefficients up to a given order: if $P$ is even, then the highest order
Taylor series coefficients that we need are those $c^{\mathbf{q}}$ such that
$|\mathbf{q}|=P$ and $q_{x},q_{y}$ are both even. This follows from the fact
that our operator does not involve any mixed derivatives: if $|\mathbf{q}|=P$
for $P$ even, then if $q_{x}$ is odd, $q_{y}$ must be as well, so taking
derivatives in only one of the dimensions will leave at least one of these
powers odd. Although we will get even cross terms from the flux divergence,
the even cross terms that arise from combining derivatives of odd order $P$
moments with lower order moments are on the order $P$ and are therefore not
needed. These cancellations are typical for centered finite differences, but
here they arrive through symmetries in moments and polynomial coefficients. As
shown in Figure 3, these coefficients can be supported with a stencil
footprint consisting of cells whose centroids are a Manhattan distance of
$\frac{P}{2}h$ from the center cell. The columns of the moment matrix
$\mathbf{M}_{u,\mathbf{i}}$ correspond to all monomials with either
$|\mathbf{q}|<P$ or $|\mathbf{q}|=P$ and $q_{x},q_{y}$ are both even. Each row
is simply the cell-averaged moments $\mathbf{m}_{\mathbf{j}}^{T}$ for each
cell $\mathbf{j}$ in the stencil, and the resulting matrix is square.

For the matrices $\mathbf{M}_{\alpha,\mathbf{i}}$ and
$\mathbf{M}_{\beta,\mathbf{i}}$ we use the same footprint as
$\mathbf{M}_{u,\mathbf{i}}$. Construction of $\mathbf{s}_{\alpha,\mathbf{i}}$
is then straightforward. The flux divergence term is slightly more
complicated. For each of the four faces of the cell
$A_{\mathbf{i}\pm\mathbf{e}_{d}}$, if we can use the same flux stencil on each
cell’s face, then the stencil will guarantee conservation, that is that its
contribution to one cell will be the negative of its contribution to its
neighbor sharing the same face. For example, consider the face
$A_{\mathbf{i}-\mathbf{e}_{x}}$, the left hand vertical face of the cell. For
the unit normal, we have $\hat{n}_{x}=-1$ and $\hat{n}_{y}=0$, so (18) reduces
to the simpler form:

$\int_{A_{\mathbf{i}-\mathbf{e}_{x}}}\beta\nabla u\cdot\hat{\mathbf{n}}\
dA=\sum_{|\mathbf{q}|\leq
P}\left(\sum_{|\mathbf{r}|\leq(P-|\mathbf{q}|)}c_{\beta}^{\mathbf{r}}\left(-\int_{-\frac{h}{2}}^{\frac{h}{2}}q_{x}\left[-\frac{h}{2},y\right]^{\mathbf{q}+\mathbf{r}-\mathbf{e}_{x}}dy\right)\right)c_{u}^{\mathbf{q}}+O(h^{P+1})\
,$ (41)

meaning the integral in parentheses is only non-zero if $q_{x}>0$ and
$q_{y}+r_{y}$ is even. In addition, for accuracy requirements we only need
moments with $\left|\mathbf{r}+\mathbf{q}-1\right|<P$. This means we do not
need to calculate any coefficients $\mathbf{c}_{\beta}^{\mathbf{r}}$ such that
$|\mathbf{r}|=P$; in other words, an order $P-1$ Taylor approximation to
$\beta$ will suffice. Formally, we can multiply the matrix
$\mathbf{M}_{\beta,\mathbf{i}}$ on the left and right by matrices that zero
out the proper columns and rows, giving us the modified moment matrix:

$\displaystyle\mathbf{M}_{\beta,A_{\mathbf{i}-\mathbf{e}_{x}}}=\mathbf{P}_{L}\mathbf{M}_{\beta,\mathbf{i}}\mathbf{P}_{R}\
,$ (42)

where $\mathbf{P}_{R}$ eliminates all columns corresponding to moments with
order greater than $P-1$, as well as columns corresponding to moments of order
$P-1$ such that $q_{y}$ is odd. $\mathbf{P}_{L}$ eliminates rows that are not
necessary to support these moments. We can likewise adjust
$\mathbf{d}_{\beta,\mathbf{i}}$ and $\mathbf{G}_{\beta,\mathbf{i}}$ to account
for these modifications. The contributions from $c_{u}^{\mathbf{q}}$ are
already symmetric about the face because they all represent centered
differences. This process results in a flux stencil footprint consisting of
all cells whose centroids are a Manhattan distance of $\frac{P-1}{2}h$ from
the centroid of the face $A_{\mathbf{i}-\mathbf{e}_{x}}$ (see Figure 3).
However, these simplifications to achieve the minimal stencil footprint rely
on symmetry arguments, so they do not apply generally to irregular and cut
cells.

$V_{+,\mathbf{i}}$$V_{-,\mathbf{i}}$$\mathcal{N}_{+,\mathbf{i}}$$\mathcal{N}_{-,\mathbf{i}}$
(a)

$V_{+,\mathbf{j}}$$V_{-,\mathbf{k}}$$\mathcal{N}_{+,\mathbf{j}}$$\mathcal{N}_{-,\mathbf{k}}$
(b)

Figure 4: Figure (a) shows neighborhoods $\mathcal{N}_{\pm,\mathbf{i}}$ used
to construct interpolation matrices for order $P=2$, around cut cell
$\mathbf{i}$, which contains two volumes $V_{\pm,\mathbf{i}}$ bordering the
interface. Figure (b) shows neighborhoods surrounding full cells $\mathbf{j}$
and $\mathbf{k}$ that don’t contain the interface, but are “irregular,”
meaning the regular stencil for order $P=2$ would be inconsistent ($O$th-
order).

####  3.1.2 Irregular Cells

Let cell $\mathbf{i}$ be an irregular cell in phase $p$, meaning its regular
cell footprint contains at least one cell which is intersected by the EB. This
invalidates the truncation error analysis for the regular cell stencil, so we
adopt a more general method for construction of moment matrices for irregular
cells. Let $\mathcal{N}_{p,\mathbf{i}}$ be a neighborhood of cells in phase
$p$ around cell $\mathbf{i}$ (see Figure 4(b) ). Our data vector
$\mathbf{d}_{u,p,\mathbf{i}}$ will consist of cell averaged values
$\left\langle u\right\rangle_{\mathbf{j}}$ for each cell
$\mathbf{j}\in\mathcal{N}_{p,\mathbf{i}}$. Each row of the corresponding
moment matrix $\mathbf{M}_{u,p,\mathbf{i}}$ will simply be
$\mathbf{m}_{\mathbf{j}}^{T}$, the row of cell averaged volume moments for
each cell $\mathbf{j}$. For irregular cells, the moment matrices
$\mathbf{M}_{\alpha,p,\mathbf{i}}$ and $\mathbf{M}_{\beta,p,\mathbf{i}}$ are
identical. We form one of each per cell and use a different
$\mathbf{G}_{\beta,p,\mathbf{i}\pm\mathbf{e}_{d}}$ for each face to create a
stencil for the flux integral along a face. The columns of these moment
matrices consist of all moments up to order $P$.

####  3.1.3 Cut Cells

In cut cells, the moment matrix is additionally used to enforce jump
conditions. When interpreted as an overdetermined system, our interpolating
polynomials are being constrained to satisfy interface matching conditions.
Alternatively as an underdetermined system, we use jump conditions as data to
cancel truncation error terms. It is necessary to enforce these jump
conditions so that our discrete operator is not degenerate.

Let $\mathcal{N}_{p,\mathbf{i}}$ be two neighborhoods of cells in their
respective phases around the cut cell $\mathbf{i}$. See Figure 4 (a). As with
irregular cells, we can form two moment matrices $\mathbf{M}_{u,p,\mathbf{i}}$
and data vectors $\mathbf{d}_{u,p,\mathbf{i}}$ of cell averaged values. For
each cut cell $\mathbf{j}$ in
$\mathcal{N}_{+,\mathbf{i}}\cup\mathcal{N}_{-,\mathbf{i}}$, we want to enforce
the two jump conditions (9) and (10). Expressing these in terms of moments and
Taylor coefficients, we have

$\displaystyle\int_{A_{B,\mathbf{j}}}u^{+}-u^{-}\ dA$
$\displaystyle=\sum_{|\mathbf{q}|\leq
P}(c_{u,+}^{\mathbf{q}}-c_{u,-}^{\mathbf{q}})m_{B,\mathbf{j}}^{\mathbf{q}}+O(h^{P+2})\
,$ (43)

and for the jump in the flux:

$\int_{A_{B,\mathbf{j}}}(\beta^{+}\nabla u^{+}-\beta^{-}\nabla
u^{-})\cdot\hat{\mathbf{n}}\ dA=\sum_{|\mathbf{r}|,|\mathbf{q}|\leq
P}\left(c_{\beta,+}^{\mathbf{r}}c_{u,+}^{\mathbf{q}}-c_{\beta,-}^{\mathbf{r}}c_{u,-}^{\mathbf{q}}\right)\left[q_{x}m_{B,\mathbf{j},x}^{\mathbf{q}+\mathbf{r}-\mathbf{e}_{x}}+q_{y}m_{B,\mathbf{j},y}^{\mathbf{q}+\mathbf{r}-\mathbf{e}_{y}}\right]+O(h^{P+1})\
.$ (44)

These expressions are both linear in the coefficients $c_{u,p}^{\mathbf{q}}$,
and they couple both sets of coefficients by interpolating $u^{+}$ and $u^{-}$
simultaneously. The resulting matrices $\mathbf{M}_{J,+,\mathbf{i}}$ and
$\mathbf{M}_{J,-,\mathbf{i}}$ have rows with the jump condition expressions
and columns corresponding to $\mathbf{c}_{u,+}$ and $\mathbf{c}_{u,-}$,
respectively. The vector $\mathbf{d}_{J,\mathbf{i}}$ consists of the given
jump condition data. Finally we form the moment matrix
$\mathbf{M}_{\mathbf{i}}$ which is defined as:

$\displaystyle\mathbf{M}_{\mathbf{i}}=\begin{bmatrix}\mathbf{M}_{u,+,\mathbf{i}}&0\\\
0&\mathbf{M}_{u,-,\mathbf{i}}\\\
\mathbf{M}_{J,+,\mathbf{i}}&\mathbf{M}_{J,-,\mathbf{i}}\end{bmatrix}\ ,$ (45)

and we solve for both sets of coefficients simultaneously:

$\displaystyle\begin{bmatrix}\mathbf{c}_{u,+}\\\
\mathbf{c}_{u,-}\end{bmatrix}$
$\displaystyle=\mathbf{M}_{\mathbf{i}}^{\dagger}\begin{bmatrix}\mathbf{d}_{u,+,\mathbf{i}}\\\
\mathbf{d}_{u,-,\mathbf{i}}\\\
\mathbf{d}_{J,\mathbf{i}}\end{bmatrix}=\begin{bmatrix}(\mathbf{M}_{\mathbf{i}}^{\dagger})_{V}&(\mathbf{M}_{\mathbf{i}}^{\dagger})_{J}\end{bmatrix}\begin{bmatrix}\mathbf{d}_{u,+,\mathbf{i}}\\\
\mathbf{d}_{u,-,\mathbf{i}}\\\
\mathbf{d}_{J,\mathbf{i}}\end{bmatrix}=(\mathbf{M}_{\mathbf{i}}^{\dagger})_{V}\begin{bmatrix}\mathbf{d}_{u,+,\mathbf{i}}\\\
\mathbf{d}_{u,-,\mathbf{i}}\end{bmatrix}+(\mathbf{M}_{\mathbf{i}}^{\dagger})_{J}\mathbf{d}_{J,\mathbf{i}}\
.$ (46)

If values in $\mathbf{d}_{J,\mathbf{i}}$ are nonzero, then forming a stencil
using these coefficients will result in adding a scalar to the right hand side
$\left\langle f\right\rangle_{p,\mathbf{i}}$. For example, for the linear term
we would have:

$\displaystyle\mathbf{s}_{\alpha,p,\mathbf{i}}^{T}$
$\displaystyle=\left[\mathbf{d}_{\alpha,p,\mathbf{i}}^{T}(\mathbf{M}_{\alpha,p,\mathbf{i}}^{\dagger})^{T}\mathbf{G}_{\alpha,p,\mathbf{i}}^{T}(\mathbf{M}_{\mathbf{i}}^{\dagger})_{V}\right]+\left[\mathbf{d}_{\alpha,p,\mathbf{i}}^{T}(\mathbf{M}_{\alpha,p,\mathbf{i}}^{\dagger})^{T}\mathbf{G}_{\alpha,p,\mathbf{i}}^{T}(\mathbf{M}_{\mathbf{i}}^{\dagger})_{J}\mathbf{d}_{J,\mathbf{i}}\right]\
,$ (47)

where the second term in brackets is a scalar. To form the moment matrices
$\mathbf{M}_{\alpha,p,\mathbf{i}}$ and $\mathbf{M}_{\beta,p,\mathbf{i}}$, we
can just use the same neighborhoods as the $\mathbf{M}_{u,p,\mathbf{i}}$; we
do not need to couple these systems because there are no external constraints
on the jumps in coefficients.

####  3.1.4 Conservation

For each cut cell and irregular cell, we have shown how to obtain a stencil
for $\int_{A}\beta\nabla u\cdot\hat{\mathbf{n}}\ dA$ on each face $A$ of any
cell. However, in order to have a conservative method, we must have only one
flux stencil per non-EB face. A simple solution to this problem is to average
the flux stencils between pairs of neighboring irregular and cut cells.
Although this results in larger stencils, it has the advantage of coupling a
layer of irregular cells to the interface jump conditions. This is because the
irregular cell stencils that border cut cells will share stencil information
with cut cells, which incorporate interface jump conditions. We reiterate that
this is not a significant issue because the density of the linear system is
dominated by the size of the regular cell stencil. For regular cells we do not
have to average with neighbors because our flux stencils were created
individually for each face and are symmetric about that face. At a cell face
which is shared between an irregular and regular cell, we use the regular cell
flux stencil.

###  3.2 Neighborhood Selection and Weighting

In general, neighborhoods need to be chosen so that resulting moment matrices
are overdetermined. So as to not perform some sort of search based on local
geometry, we opt to make the neighborhood sufficiently large to accommodate a
reasonably smooth geometry. For any irregular or cut cell in phase $p$, we let
$\mathcal{N}_{p,\mathbf{i}}$ be those cells in phase $p$ that lie in the the
square of cells with side length $2P+1$ surrounding cell $i$. See Figure 4. As
is done in [12] and [15], we employ a weighted least-squares approach to force
stencil weights to decay with distance faster than the growth of the highest
polynomial term. To each piece of data in $\mathbf{d}$ we assign a weight that
is inversely related to its distance from the centroid of volume
$V_{p,\mathbf{i}}$. If $\delta_{\mathbf{j}}$ is this distance, then the
corresponding row of the moment matrix and $\mathbf{d}$ are multiplied by
$w_{\mathbf{j}}$, where

$\displaystyle w_{\mathbf{j}}=\frac{1}{(1+\delta_{\mathbf{j}})^{P+1}}\ ,$ (48)

where $P$ is the order of the scheme. This forms a diagonal weight matrix
$\mathbf{W}$, which gives us the weighted least-squares solution:

$\displaystyle\mathbf{c}_{u}=(\mathbf{W}\mathbf{M})^{\dagger}\mathbf{W}\mathbf{d}\
.$ (49)

This weighting does not affect the truncation error: with some matrix algebra
we can see that this weighted least-squares solution is equivalent to a change
of basis in the undetermined formulation:

$\displaystyle\mathbf{s}$
$\displaystyle=\mathbf{W}\left(\mathbf{M}^{T}\mathbf{W}\right)^{\dagger}\mathbf{g}\implies$
(50) $\displaystyle\mathbf{M}^{T}\mathbf{W}\mathbf{W}^{-1}\mathbf{s}$
$\displaystyle=\mathbf{g}\ .$ (51)

Since $\mathbf{M}^{T}\mathbf{s}=\mathbf{g}$ still holds, the truncation error
is unaffected by the weighting. This has proved to be an effective tool for
controlling the spectrum and conditioning of the discrete operator [12]. We
compute the pseudoinverse using the SVD algorithm in _LAPACK_ [16].

###  3.3 Solver and Software Implementation

We assemble the stencils to form the linear system

$\displaystyle\mathbf{L}\mathbf{u}=\mathbf{f}+\mathbf{r}\ ,$ (52)

where $\mathbf{r}$ represents the contribution to the right hand side from the
jump conditions. Following [12], we precondition this system by left
multiplying with the diagonal matrix whose $\left(i,i\right)$ entry is
$\frac{|V_{\mathbf{i}}|}{h^{2}}$; i.e. we multiply each row by the volume
fraction of that cell. This simple preconditioner eliminates the volume
scaling associated with very small volume fractions.

We solve the linear system using Krylov subspace methods and preconditioners
provided by the _PETSc_ library [17], [18]. Since the linear system is non-
symmetric, we use BiCG-Stab or GMRES. We have experimented with the _PETSc_
algebraic multigrid and block Jacobi preconditioners. One of these options is
typically sufficient, but if they fail, we use the direct solver _SuperLU_
[19]. In future research we will develop a geometric multigrid preconditioner
similar to that in [8] or [20]. Our method is well-suited for geometric
multigrid because the Taylor series formulation makes interpolation
straightforward. However, we emphasize that an efficient solver is not the
focus of this particular paper. The algorithm is implemented using the
_Chombo_ software library [21], which allows for large-scale parallelization
of the algorithm. Visualizations are created using _VisIt_ [22].

##  4 Numerical Tests

We validate our method with a series of numerical tests. The goals of this
section are to:

  1. 1.

Validate the truncation error analysis and measure the solution error,

  2. 2.

Demonstrate consistent accuracy for problems with large coefficient and
solution jumps, and

  3. 3.

Demonstrate convergence on non-trivial interface geometries.

Although our scheme is designed to be have arbitrary order of accuracy, we
have evaluated it for just $P\in\left\\{2,4,6\right\\}$. We measure the error
$\left\langle e\right\rangle_{\mathbf{i}}$ as discrete cell averages, and
evaluate it using discrete $L^{p}$ norms:

$\displaystyle\left\|e\right\|_{1}$
$\displaystyle=\int_{\Omega}|e|dV=\sum_{p,\mathbf{i}}\left|\left\langle
e\right\rangle_{p,\mathbf{i}}\right||V_{p,\mathbf{i}}|$ (53)
$\displaystyle\left\|e\right\|_{\infty}$
$\displaystyle=\max_{p,\mathbf{i}}\left\langle e\right\rangle_{p,\mathbf{i}}\
.$ (54)

(a)

(b)

(c)

Figure 5: Results of the ellipse boundary tests for $P=4$, $n=512$. Plots of
the (a) exact solution, (b) absolute value of truncation error (log scale),
and (c) absolute value of solution error (log scale). Note that the truncation
error is concentrated at the interface, and its influence on the solution
error is much smoother after inverting the elliptic operator.

3264128256512$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$10^{2}$$n$ErrorTruncation
Error, $L^{\infty}$ Norm$P=6$ $O(h^{5})$$P=4$ $O(h^{3})$$P=2$ $O(h)$ (a)

3264128256512$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$10^{2}$$n$ErrorTruncation
Error, $L^{1}$ Norm$P=6$ $O(h^{6})$$P=4$ $O(h^{4})$$P=2$ $O(h^{2})$ (b)

3264128256512$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$n$ErrorSolution
Error, $L^{\infty}$ Norm$P=6$ $O(h^{6})$$P=4$ $O(h^{4})$$P=2$ $O(h^{2})$ (c)

3264128256512$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$n$ErrorSolution
Error, $L^{1}$ Norm$P=6$ $O(h^{6})$$P=4$ $O(h^{4})$$P=2$ $O(h^{2})$ (d)

Figure 6: Truncation error $L^{\infty}$ norm (a), $L^{1}$ norm (b), and
solution error norms (c) and (d), respectively, for the ellipse geometry. We
observe the expected rate of convergence for both quantities in both norms.
The truncation error is dominated by behavior at the interface, which because
it is codomension-1 smaller, achieves one order higher in $L^{1}$ norm (see
Figure 5).

###  4.1 Truncation and Solution Error Validation

Our physical domain for all tests is $\Omega=[-1,1]^{2}$, intersected by some
interface $\Gamma$. Our first interface $\Gamma$ is an ellipse with major axis
of length $\frac{\pi}{4}$ and minor axis of length $\frac{\pi}{8}$. The
superscript $+$ refers to quantities enclosed by the interface and the
superscript $-$ refers to quantities on the exterior of the interface. We test
our discretization using the method of manufactured solutions, such that
$\alpha^{\pm},\beta^{\pm},u^{\pm}$ are all constructed as linear combinations
of periodic functions $p_{k_{x},k_{y}}$ on the square:

$\displaystyle p_{k_{x},k_{y}}=\cos^{2}\left(\pi
k_{x}x\right)\sin^{2}\left(\pi k_{y}y\right)\ ,$ (55)

and let $k_{x},k_{y}\in\left\\{-2,-1,1,2\right\\}$, creating 20 total basis
functions. We randomly generate different sets of coefficients
$c_{k_{x},k_{y}}\in[-1,1]$ for each of the six functions
$\alpha^{\pm},\beta^{\pm},u^{\pm}$. An exact $f^{\pm}$ is formed by applying
the exact differential operator to $u^{\pm}$. The variable coefficient fields
are offset by a positive constant so that they are nonnegative everywhere.
Figure 5(a) shows the solution for one particular example.

The global truncation error $\mathbf{e}$ and solution error $\mathbf{t}$ are
defined as

$\displaystyle\mathbf{e}$ $\displaystyle=\mathbf{u}^{e}-\mathbf{u}$ (56)
$\displaystyle\mathbf{t}$
$\displaystyle=\left(\mathbf{L}\mathbf{u}^{e}-\mathbf{r}\right)-\mathbf{L}^{e}\mathbf{u}^{e}\
,$ (57)

where $\mathbf{u}^{e}$ is the exact solution and $\mathbf{L}^{e}$ is the exact
operator, so that the solution error satisfies the equation

$\displaystyle\mathbf{L}\mathbf{e}=\mathbf{t}\ .$ (58)

Note that (58) implies homogeneous jump conditions on the error; there is no
contribution to the right-hand-side of this system from jump conditions.

Our truncation error analysis predicts an order $P-1$ truncation error in cut
and irregular cells and an order $P$ truncation error in regular cells. In
Figure 5(b) we see that truncation error (for $P=4$) is almost entirely
concentrated in cut and irregular cells. This is reflected in Figure 6(a), as
the max norm of the truncation error converges at order $P-1$, the expected
rate for cells at and near the interface. However, the number of cut and
irregular cells is of order $h^{-1}$ because it is a codimension one smaller
region, while the number of regular cells is scales like $h^{-2}$. Therefore
for the $L^{1}$ norm we have:

$\displaystyle\left\|t\right\|_{1}$
$\displaystyle=\sum_{\mathbf{i}\in\Omega_{C}\cup\Omega_{I}}\left|\left\langle
t\right\rangle_{p,\mathbf{i}}\right||V_{p,\mathbf{i}}|+\sum_{\mathbf{i}\in\Omega_{R}}\left|\left\langle
t\right\rangle_{p,\mathbf{i}}\right||V_{p,\mathbf{i}}|=\sum_{\mathbf{i}\in\Omega_{C}\cup\Omega_{I}}O(h^{P+1})+\sum_{\mathbf{i}\in\Omega_{R}}O(h^{P+2})$
(59) $\displaystyle=O(h^{-1})O(h^{P+1})+O(h^{-2})O(h^{P+2})=O(h^{P})\ .$ (60)

This is confirmed in 6 (b); we see clean order $P$ convergence for the $L^{1}$
norm of the truncation error. Given that the truncation error at the interface
is orders of magnitude greater than truncation error elsewhere, the $L^{1}$
norm of the truncation error is also dominated by behavior at the interface.

Although we do not have an analytical bound on
$\left\|\mathbf{L}^{-1}\right\|$, based on the analysis in [7] and the results
in [12], [8], [23] and others, we expect the solution error to converge at
order $P$ in both norms. This behavior is shown in Figures 5(b) and (c); the
solution error is roughly of the same order of magnitude everywhere in the
domain. We plan on further analysis to explore the combined effects of the
homogeneous jump conditions imposed on the error equation and the regularity
of the elliptic operator, but the empirical results demonstrate the desired
convergence rates.

1e-41e-21e01e21e4$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$\frac{\beta^{-}}{\beta^{+}}$ErrorSolution
Error, $L^{1}$ Norm$P=2$$P=4$$P=6$ (a)

1e-41e-21e01e21e4$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$10^{2}$$\frac{\beta^{-}}{\beta^{+}}$ErrorSolution
Error, $L^{1}$ Norm$P=2$$P=4$$P=6$ (b)

Figure 7: Solution error for discontinuous diffusion coefficient tests. Figure
(a) is for the cosine geometry, whereas (b) is the annulus geometry (we only
plot $L^{1}$ error because the $L^{\infty}$ norms behave nearly identically).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8: Plots of the absolute value of truncation error (left column) and
absolute value of solution error (right column) for $P=2$, $n=512$, on a log
scale. Figures show: (a) and (b) cosine geometry with
$\frac{\beta^{-}}{\beta^{+}}=10^{-4}$; (c) and (d) annulus geometry with
$\frac{\beta^{-}}{\beta^{+}}=10^{4}$; and (e) and (f) annulus geometry with
$\frac{\beta^{-}}{\beta^{+}}=10^{-4}$. Note that the larger, rough trunction
error near the interface becomes smoother, rapidly decaying error in the
solution (except (f), see text for discussion).

###  4.2 Discontinuous Diffusion Coefficient

Next we test the robustness of our scheme on problems with large jumps in the
diffusion coefficient, as is in common in the literature (see [8], [23], [11],
[5]). We set the linear term coefficient $\alpha^{\pm}=0$, and let the
diffusion coefficients be constant, varying the ratio
$\frac{\beta^{-}}{\beta^{+}}$ from $10^{4}$ to $10^{-4}$. Specifically, we fix
$\beta^{+}=1$ and vary $\beta^{-}$ from $1$ to $10^{4}$, and vice versa. The
manufactured solution is the same as in the truncation and solution error
tests. We are interested in studying the relationship between the solution
error and the ratio of diffusion coefficients, so we fix the grid spacing at
$h=256^{-1}$. However, we introduce two geometries that expose different error
characteristics.

The first interface geometry is the zero level set of the function

$\displaystyle\psi(x,y)=\frac{1}{4}\cos(\pi y)+x+\frac{\pi}{100}\ ,$ (61)

which is simply a cosine in the xy plane. We impose periodic boundary
conditions in the $y$ direction and Dirichlet boundary conditions in the $x$
direction. The $+$ region is the to the right of the interface for this
geometry. We impose Dirichlet boundary conditions by filling layers of ghost
cells with exact solution values. In this case both phases are tied to
boundary conditions, so we do not expect any significant difference between
large and small value of $\frac{\beta^{-}}{\beta^{+}}$. For the second test,
we use the “annulus” interface geometry given in section 3.1 of [23]. This
interface shape can be seen in Figures 8(c) and (e). In the former case, the
$+$ (interior) phase has a relatively large coefficient ratio ($10^{4}$),
while in the latter it is the inverse ($10^{-4}$). Given that in these cases
the $-$ phase has domain boundary conditions while the $+$ phase does not, we
expect to see different error behavior as $\frac{\beta^{-}}{\beta^{+}}$
varies.

For the cosine geometry, we see that accuracy is mostly unaffected by changes
in the diffusion coefficient ratio (Figure 7(a)). In Figure 8(a) and (b), the
truncation error again is larger at the interface but the solution error is
smooth, similar to the previous tests. We observed no impact of the
conditioning of the discretization matrix on the solution accuracy as the
ratio of coefficients varies over 8 orders of magnitude.

For the second test, instead we see that the error is about 4 orders of
magnitude higher when we have the diffusion coefficient on the interior of the
domain is much larger than the diffusion coefficient on the exterior of the
domain ($\beta^{-}/\beta^{+}=10^{-4}$). This result is consistent with results
reported in Figure 2 of [8], as well as Figure 17b of [11]. Through potential
theoretic arguments, we believe this is a result of solution errors in the
interior region not being “tied down” to any domain boundary condition, as in
the cosine test. Because of the gradient jump condition (3), any gradient
_errors_ in the interior are multiplied by $\beta^{+}/\beta^{-}=10^{4}$ in
their contribution to the exterior domain gradients at the interface, forcing
the interior solution there to “drift” in proportion. However, with this
scaling the convergence rates are still retained, but with an error constant
reflecting this ratio in diffusion coefficients.

###  4.3 Discontinuous Solution

We perform a similar test with a solution that has large jumps at the
interface. $\beta^{\pm}$ and $u^{\pm}$ are the same as in the solution and
truncation error test, and we again set $\alpha^{\pm}=0$. The $u^{\pm}$ fields
are multiplied by scaling factors $s^{\pm}$ to create large jumps in the
solution, and we use the same two geometries as the discontinuous coefficients
test. We observe $L^{1}$ and $L^{\infty}$ errors that are proportional to the
larger of the two scaling coefficients (see Figure 9). This scaling does not
magnify the error because it appears as a large discontinuity in the source
term, as well as in the jump conditions, which both contribute only to the
right hand side of the linear system. This test highlights the importance of
the having two separate degrees of freedom in each cut cell, from which we are
able to accurately reproduce a solution and gradients which jump by up to 4
orders of magnitude across the interface.

1e-41e-21e01e21e4$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$10^{2}$$10^{4}$$\frac{s^{-}}{s^{+}}$ErrorSolution
Error, $L^{\infty}$ Norm$P=2$$P=4$$P=6$ (a)

1e-41e-21e01e21e4$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$10^{2}$$10^{4}$$\frac{s^{-}}{s^{+}}$ErrorSolution
Error, $L^{\infty}$ Norm$P=2$$P=4$$P=6$ (b)

Figure 9: Comparison of solution maximum errors for two geometries (a) cosine,
and (b) annulus, and over a range of solution scaling ratios, fixing grid size
at $n=128$ cells and varying the order of accuracy $P$ (note that the $L^{1}$
error is nearly identical).

3264128256$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$n$ErrorSolution Error,
$L^{\infty}$ Norm$\frac{\beta^{-}}{\beta^{+}}=10^{-3}$
$\frac{\beta^{-}}{\beta^{+}}=10^{3}$ $\frac{\alpha^{-}}{\alpha^{+}}=10^{-3}$
$\frac{\alpha^{-}}{\alpha^{+}}=10^{3}$ $\frac{f^{-}}{f^{+}}=10^{-3}$
$\frac{f^{-}}{f^{+}}=10^{3}$$O(h^{4})$$O(h^{4})$ (a)

3264128256$10^{-10}$$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$$n$ErrorSolution
Error, $L^{1}$ Norm (b)

Figure 10: $L^{\infty}$ and $L^{1}$ norm of errors for Richardson convergence
test for a variety of scaling factor ratios in linear term coefficient,
diffusion coefficient, and source term. We fix $P=4$ and observe fourth order
convergence in the $L^{1}$ norm and nearly fourth order convergence in the max
norm.

###  4.4 Imposing Homogeneous Jump Conditions

Lastly, we test the ability of the method to impose jump conditions as a
constraint. We let $\alpha^{\pm},\beta^{\pm}$ be the same as in the truncation
and solution error test, and use the manufactured solution $u^{\pm}$ from that
test as the source term $f^{\pm}$. We impose homogeneous (zero) jump
conditions and focus on the annulus geometry. For this test we fix $P=4$ and
test the scheme with a variety of coefficient and source term scalings. The
error is measured by using the $n=512$ numerical solution as the exact
solution, with the results in Figure 10. We observe roughly fourth order
convergence for all tests, although there is more variation in convergence
compared to the manufactured solution tests.

##  5 Conclusion

We have developed a finite volume method for the variable coefficient elliptic
interface problem and demonstrated up to sixth order accurate on a variety of
test problems. In developing this method we gave a general truncation error
analysis that justifies the use of stencils based on least-squares
interpolation. Our stencils are derived from cell-centered Taylor polynomials
which are implicitly defined in terms of local values of the solution and
interface jumps or boundary conditions, where appropriate. To enforce
conservation, we choose a single flux on each face which is an average of the
flux calculated in neighboring cells from the respective Taylor polynomials.
In cells away from the interface, we take advantage of standard finite volume
symmetries to build stencils with a minimal footprint.

Future research will involve 1) building an efficient geometric multigrid
solver along the lines of [20], 2) extending our method to three dimensions,
and 3) incorporating adaptive mesh refinement. The prior method presented in
[12] accomplishes these for Poisson’s equation in more complex geometries,
including boundaries with kinks, using “smoothed” constructive solid geometry
capabilities of the Chombo software library [21]. This would enable this
method to be used in discretizations for large scale science applications.
Further exploration is also needed of the theory of undetermined stencil
systems (using the moment matrix transpose, $\mathbf{M}^{T}$). This paper has
shown that building a stencil of a given order and truncation error still
allows infinitely many valid stencils; this fact could be exploited to promote
sparsity or alter the conditioning or stability of the operator, and we are
drafting a paper with analysis that may provide specific algorithmic guidance.

## Funding:

This work is supported by the U.S. Department of Energy, Office of Science,
Advanced Scientific Computing Research, Base Math Program, under contract
number DE-AC02-05CH11231.

##  Appendix A Geometry Generation Algorithm

We specify the interface as a zero level set of an implicit function
$\psi(\mathbf{x})$. Therefore, for sufficiently smooth $\phi$ compared to the
grid resolution, we assume we can identify cut cells by evaluating $\psi$ at
the four corners of each cell. If any of these values have different signs,
the cell is tagged as a cut cell. If the interface intersects one face of the
cell multiple times, or there are more than two faces intersected by the
interface, we consider the geometry to be under-resolved and could refine the
mesh or adjust the boundary without inducing significant errors. Given these
assumptions, when the interface intersects a cell it creates a region which is
bounded on one side by the interface and on two or three sides by the edges of
a square cell (see Figure 1). Volume moments (4) are defined as integrals over
this region. Area moments, defined in (6) and (7), are integrals over the
portion of the EB that intersects the cut cell. We compute these integrals by
approximating the interface with piecewise line segments, and then apply a
formula for the integral of monomials along line segments. The vertices of the
line segments are roots of $\psi$, which we find using a simple root finder
such as the secant method. By refining this interface iteratively into $2^{n}$
line segments, we can calculate a convergent sequence $m_{n}$ of moment
approximations that stops when $\left|m_{n+1}-m_{n}\right|$ reaches machine
precision. The convergence of this sequence is accelerated using Richardson
extrapolation, which in this case is often referred to as Romberg integration.

A formula for the integral of $x^{p}y^{q}$ over an arbitrary polygon can be
derived from Green’s theorem:

$\displaystyle\int_{P}\frac{\partial f(x,y)}{\partial
x}dV=\int_{C}f(x,y)\hat{n}_{x}dA=\int_{C}f(x,y)dy\ ,$ (62)

where $C$ is the boundary of the polygon $P$. Let
$f=\frac{x^{p+1}}{p+1}y^{q}$, giving us:

$\displaystyle\int_{P}x^{p}y^{q}dV=\int_{C}\frac{x^{p+1}}{p+1}y^{q}dy\ .$ (63)

We parameterize each edge segment by:

$\displaystyle x(t)$ $\displaystyle=(x_{k+1}-x_{k})t+x_{k}=\Delta xt+x_{k}$
(64) $\displaystyle y(t)$ $\displaystyle=(y_{k+1}-y_{k})t+y_{k}=\Delta
yt+y_{k}\ ,$ (65)

where $t$ goes from $0$ to $1$ and $(x_{k},y_{k})$ are the ordered vertices of
the polygon. The formula for the right hand side of (63) along a single line
segment $C_{k}(t)$ is obtained by a binomial expansion:

$\displaystyle\int_{C_{k}}\frac{x^{p+1}}{p+1}y^{q}dy$
$\displaystyle=\frac{1}{p+1}\int_{0}^{1}(\Delta xt+x_{k})^{p+1}(\Delta
yt+y_{k})^{q}\Delta ydt$ (66)
$\displaystyle=\sum_{i=0}^{p+1}\sum_{j=0}^{q}\frac{{p+1\choose i}{q\choose
j}}{(p+1)(p+2+q-i-j)}(x_{k}^{i}y_{k}^{j})(\Delta x)^{p+1-i}(\Delta y)^{q+1-j}\
.$ (67)

We follow a similar procedure for area integrals:

$\displaystyle\int_{C}x^{p}y^{q}$
$\displaystyle=\int_{0}^{1}x(t)^{p}y(t)^{q}dC$ (68)
$\displaystyle=\sum_{i=0}^{p}\sum_{j=0}^{q}\frac{{p\choose i}{q\choose
j}}{(p+1+q-i-j)}(x_{k}^{i}y_{k}^{j})(\Delta x)^{p-i}(\Delta y)^{q-j}\Delta C\
,$ (69)

and to calculate area integrals times unit normals we multiply equation (69)
by $n_{x}=\dfrac{\Delta y}{\Delta C}$ or $n_{y}=\dfrac{-\Delta x}{\Delta C}$,
where $\Delta C=\sqrt{\Delta x^{2}+\Delta y^{2}}$. (The tangent vector is
rotated 90 degrees clockwise). The integrals of interest are obtained by
adding the integrals along all line segments of the polygon in the case of
volume moments, or just along the interface in the case of area moments.

## References

  * Li [2003] Z. Li,  An overview of the immersed interface method and its applications,  Taiwanese Journal of Mathematics 7 (2003) 1–49.
  * Gibou et al. [2013] F. Gibou, C. Min, R. Fedkiw,  High resolution sharp computational methods for elliptic and parabolic problems in complex geometries,  J. Sci. Comput. 54 (2013) 369–413.
  * Babuska [1970] I. Babuska,  The finite element method for elliptic equations with discontinuous coefficients,  Computing 5 (1970) 207–213.
  * Li [1998] Z. Li,  The immersed interface method using a finite element formulation,  Applied Numerical Mathematics 27 (1998) 253–267.
  * Leveque and Li [1994] R. J. Leveque, Z. Li,  The immersed interface method for elliptic equations with discontinuous coefficients and singular sources,  SIAM Journal on Numerical Analysis 31 (1994) 1019–1044.
  * Fedkiw et al. [1999] R. P. Fedkiw, T. Aslam, B. Merriman, S. Osher,  A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method),  Journal of Computational Physics 152 (1999) 457–492.
  * Johansen and Colella [1998] H. Johansen, P. Colella,  A cartesian grid embedded boundary method for poisson’s equation on irregular domains,  Journal of Computational Physics 147 (1998) 60–85.
  * Crockett et al. [2011] R. Crockett, P. Colella, D. Graves,  A cartesian grid embedded boundary method for solving the poisson and heat equations with discontinuous coefficients in three dimensions,  Journal of Computational Physics 230 (2011) 2451–2469.
  * SCH [2006] A cartesian grid embedded boundary method for the heat equation and poisson’s equation in three dimensions,  Journal of Computational Physics 211 (2006) 531–550.
  * Colella [2016] P. Colella,  High-order finite-volume methods on locally-structured grids,  Discrete and Continuous Dynamical Systems 36 (2016) 4247–4270.
  * Chen and Strain [2008] T. Chen, J. Strain,  Piecewise-polynomial discretization and krylov-accelerated multigrid for elliptic interface problems,  Journal of Computational Physics 227 (2008) 7503–7542.
  * Devendran et al. [2017] D. Devendran, D. Graves, H. Johansen, T. Ligocki,  A fourth-order Cartesian grid embedded boundary method for Poisson’s equation,  Communications in Applied Mathematics and Computational Science 12 (2017) 51 – 79.
  * Zhang et al. [2012] Q. Zhang, H. Johansen, P. Colella,  A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation,  SIAM Journal on Scientific Computing 34 (2012) B179–B201.
  * LeVeque [2007] R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, Society for Industrial and Applied Mathematics, 2007.
  * Overton-Katz et al. [2022] N. Overton-Katz, X. Gao, S. Guzik, O. Antepara, D. T. Graves, H. Johansen,  A fourth-order embedded boundary finite volume method for the unsteady stokes equations with complex geometries,  arXiv (2022).
  * Barker et al. [2001] V. A. Barker, L. S. Blackford, J. Dongarra, J. D. Croz, S. Hammarling, M. Marinova, J. Waśniewski, P. Yalamov, LAPACK95 Users’ Guide, Society for Industrial and Applied Mathematics, 2001.
  * Balay et al. [2022] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, J. Zhang, PETSc/TAO Users Manual, Technical Report ANL-21/39 - Revision 3.17, Argonne National Laboratory, 2022\.
  * Balay et al. [1997] S. Balay, W. D. Gropp, L. C. McInnes, B. F. Smith,  Efficient management of parallelism in object oriented numerical software libraries,  in: E. Arge, A. M. Bruaset, H. P. Langtangen (Eds.), Modern Software Tools in Scientific Computing, Birkhäuser Press, 1997, pp. 163–202.
  * Li and Demmel [2003] X. S. Li, J. W. Demmel,  SuperLU_DIST: A scalable distributed-memory sparse direct solver for unsymmetric linear systems,  ACM Trans. Mathematical Software 29 (2003) 110–140.
  * Devendran et al. [2014] D. Devendran, D. T. Graves, H. Johansen,  A hybrid multigrid algorithm for poisson’s equation using an adaptive, fourth order treatment of cut cells,  Technical Report LBNL-1004329, LBNL (2014).
  * Adams et al. [2021] M. Adams, P. Colella, D. Graves, J. Johnson, H. Johansen, N. Keen, T. Ligocki, D. Martin, P. McCorquodale, D. Modiano, P. Schwartz, T. Sternberg, B. V. Straalen, Chombo software package for AMR applications: design document, Technical Report, April 2021.
  * Childs et al. [2012] H. Childs, E. Brugger, B. Whitlock, J. Meredith, S. Ahern, D. Pugmire, K. Biagas, M. C. Miller, C. Harrison, G. H. Weber, H. Krishnan, T. Fogal, A. Sanderson, C. Garth, E. W. Bethel, D. Camp, O. Rubel, M. Durant, J. M. Favre, P. Navratil, VisIt: An End-User Tool For Visualizing and Analyzing Very Large Data, 2012.
  * Bochkov and Gibou [2020] D. Bochkov, F. Gibou,  Solving elliptic interface problems with jump conditions on cartesian grids,  Journal of Computational Physics 407 (2020) 109269.


,

# Tunneling Analysis of Kerr-Newman Black Hole-Like Solution in Rastall Theory

Riasat Ali<EMAIL_ADDRESS>Department of Mathematics, GC University
Faisalabad Layyah Campus, Layyah-31200, Pakistan  Kazuharu Bamba
<EMAIL_ADDRESS>Faculty of Symbiotic Systems Science, Fukushima
University, Fukushima 960-1296, Japan  Syed Asif Ali Shah
<EMAIL_ADDRESS>Department of Mathematics and Statistics, The
University of Lahore 1-Km Raiwind Road, Sultan Town Lahore 54000, Pakistan
Muhammad Junaid Saleem<EMAIL_ADDRESS>NCBA and E Multan Campus
Multan-66000, Pakistan

###### Abstract

Hamilton-Jacobi ansatz is used to analyze boson charged particles tunneling in
Rastall gravity Kerr- Newman black hole (BH) surrounded by perfect fluid
matter through the horizon. The geometrical BH parameters are observed, how
these are affecting the radiation utilizing the Lagrangian gravity equation.
In this article, the corrected Hawking temperature is computed by considering
the effects of quantum gravity. In our analysis, the Rastall gravity BH
solution surround by perfect fluid matter effects on Hawking radiation is
analyzed graphically. Moreover, the stability and instability of Rastall
gravity BH are investigated. The influence of quantum gravity and rotation
parameters on BH radiation are also observed by graphical representation.

Rastall Gravity Black Hole; Lagrangian Gravity Equation; Hawking Radiation

##  I Introduction

Hawking 1  observed that a black hole (BH) radiates all kinds of particles and
the radiated spectrum of particles are absolutely thermal and the spectrum of
emission probability for different types of particles has been studied. The
phenomenon of tunneling is established on the physical action of particles
theory which gives BHs radiation. The law of energy conservation and momentum
conservation has been adapted to this phenomenon. The particular temperature
value at which vector particles are emitted from the BH horizon is called
boson radiation. The Hawking radiation for scalar and fermion particles from
different spacetimes have been investigated 01 ; 2 . The tunneling probability
of charged scalars and fermions outgoing from the event horizon from these BHs
have been studied. The temperature and entropy in the cosmological horizon of
Schwarzschild de-Sitter spacetime have been investigated in 3 . The charged
particles radiated through tunneling from 4-dimensional as well as
5-dimensional BHs and Vaidya BH have been investigated 5 ; 6 , however, the
Hawking temperatures and tunneling rates which depend on the properties of
space-time have been observed.

The noncommutative acoustic BHs of quantum-corrected finite entropy have been
discussed 7  and the corrections in terms of electric charge are conserved.
The energy and mass of outgoing fermions and Kerr and Kerr-Newman BHs have
been investigated 8 ; 9 ; 10  and the tunneling particles would be stopped at
some particular range and it is also concluded that the tunneling rate of
radiated particles and Hawking temperature for both BHs which is depends on
different parameters. The Hawking radiation spectrum under the influence of
quantum gravity has been analyzed from different types of BH in x1 ; x2 . The
energy spectrum and the Hawking temperature contribute holomorphic and anti-
holomorphic functions, GUP and conformal gravity effects on tunneling
radiation from various BHs have been observed 11 ; 12 ; 13 ; 14 , which are
consistent with original results. The BHs tunneling radiation with the
assumption of GUP 15 , unitary theory 16 , loop quantum theory17  and
coordinate system 18  have been examined and found that which depends on GUP,
unitary theory, string and loop quantum theory but independent on coordinate
system.

The tunneling particles from the BH, BTZ like BH and black ring have been
examined 19 ; 20 ; 21  and it contains the hypothesis; BHs and black rings
have both unstable and stable properties. Many authors computed R1 ; R2 ; R3 ;
R4 ; R5 ; R6 ; R7 ; R8 ; R9 ; R10 ; R11 ; R11! ; R12 ; R13 ; R14 ; R15 ; R16 ;
R17 ; R18 ; R19 ; R20 ; R21 ; R22 ; R23  the tunneling radiation of boson and
fermion particles to get the Hawking temperature for different wormholes and
BHs. The researcher computed 24 ; 25  the quantum gravity effect on tunneling
radiation and it is a phenomenon for different types of BHs.

The Newmananis algorithmic rule to analyze the Hayward BH solution with
rotation parameter has been computed in N1 .The temperature for regular
Hayward BH metric with rotation parameter has been derived with the help of
semi-classical phenomenon. An extension for temperature of 3 and 4 dimensional
BH metric as well as higher dimensional BHs metric under the influences of
quantum gravity have been analyzed in N2 ; N3 . The Rastall rotating BH with
surrounding quintessence and dust parameter has been discussed in N4  and
showed that the Rastall theory metric solution is different from the standard
metric solution.

The charged bosons particles tunneling with spin equals $1$ that is $Z$ and
$W_{\pm}$ play a significant role for the standard model of electro-weak
interaction so that the radiation of charged boson particles should be
importance in the observation of tunneling radiation. Many different
approaches have been analyzed for the investigation of the tunneling spectrum
through the horizons of BHs. Visser 26  analyzed that the Rastall gravity is
utterly equivalent to standard Einstein general relativity which is obtained
the artificially splitting the physical conserved stress-energy tensor into
two non-conserved pieces.

This paper aims to study the quantum gravity effects on BH thermodynamics.
This paper is organized as follows. In the Sec. (II) deals the introduction of
Rastall gravity Kerr-Newman BH surrounded by perfect fluid matter. We analyze
the tunneling radiation of BH by assuming the Lagrangian equation in the
scenario of GUP and calculated the temperature. The graphical representation
of temperature in terms of rotation and quantum gravity are discussed in Sec.
(III). The last section deals with our results.

##  II Kerr-Newman Black Hole in Rastall Gravity Theory

By applying Newman-Janis algorithmic rule to study the rotating Rastall BH
space-time, which analyzes by mass ($M$), rotation parameter($a$), surrounding
fluid structure parameter ($\sigma$), state parameter of surrounding fluid $u$
and Rastall coupling parameter $\alpha$. The Rastall BH metric in the Boyer-
Lindquist coordinates from 27 ; J1  can be expressed as

$ds^{2}=-Adt^{2}+Bdr^{2}+Cd\theta^{2}+Dd\phi^{2}+2Fdtd\phi,$ (1)

where $A$, $B$, $C$, $D$ and $E$ are given by the following equations:

$\displaystyle A$ $\displaystyle=$ $\displaystyle 1-\frac{2Mr+\sigma
r^{\gamma}}{\Sigma},~{}~{}~{}B=\frac{\Sigma}{\lambda_{r}},~{}~{}~{}C=\Sigma,$
$\displaystyle D$ $\displaystyle=$ $\displaystyle
sin^{2}\theta\left(r^{2}+a^{2}+\frac{(2Mr+\sigma
r^{\gamma})a^{2}sin^{2}\theta}{\Sigma}\right),$ $\displaystyle F$
$\displaystyle=$ $\displaystyle\frac{-asin^{2}\theta\left(2Mr+\sigma
r^{\gamma}\right)}{\Sigma}.$

Here, these $\lambda_{r}$, $\Sigma$ and $\gamma$ parameters are defined as

$\displaystyle\lambda_{r}$ $\displaystyle=$ $\displaystyle
r^{2}-2Mr+a^{2}-\sigma r^{\gamma},$ $\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle a^{2}cos^{2}\theta+a^{2},$ $\displaystyle\gamma$
$\displaystyle=$ $\displaystyle{\frac{1-3u}{1-3\alpha(1+u)}}$

It is important to mention that, if $a=0$ and $-1<u<\frac{-1}{3}$, the metric
(1) constitutes the Kerr BH surrounded by quintessence J2  and if $\sigma$ and
$a$ both are zero, we get Schwarzschild BH solution in Ref. J1 .

In order to study the tunneling radiation of charged particles via the BH
horizon, we shall analyze electromagnetic effects of the Lagrangian gravity
equation. The electromagnetic-charged fields, which describe the particle’s
motion in the Lagrangian charged equation 29

$\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}\Psi^{\nu\mu})+\frac{i}{h}eA_{\mu}\Psi^{\nu\mu}+\frac{m^{2}}{h^{2}}\Psi^{\nu}+\hbar^{2}\beta\partial_{0}\partial_{0}\partial_{0}(g^{00}\Psi^{0\nu})-\hbar^{2}\beta\partial_{i}\partial_{i}\partial_{i}(g^{ii}\Psi^{i\nu})=0,$
(2)

where $g$, $\Psi^{\mu\nu}$ and $m$ are the matrix determinant of coefficients,
particle mass and anti-symmetric tensor, respectively.

$\Psi_{\nu\mu}=(1-\hbar^{2}\beta\partial^{2}_{\nu})\partial{\nu}\Psi_{\mu}-(1-\hbar^{2}\beta\partial^{2}_{\mu})\partial{\mu}\Psi_{\nu}+(1-\hbar^{2}\beta\partial^{2}_{\nu})\frac{i}{h}eA_{\nu}\Psi_{\mu}-(1-\hbar^{2}\beta\partial^{2}_{\mu})\frac{i}{h}eA_{\mu}\Psi_{\nu}.$

Here, $\nabla_{\mu}$, $e$ and $A_{\mu}$ are the covariant derivative, particle
charged and BH potential, respectively. So the tunneling of positive and
negative particles is similar but opposite sign $(W_{+}=-W_{-})$. For
explanation, here we shall study the positive particles tunneling and this
tunneling can be changed to negative particles tunneling. The values of
$\Psi^{\mu}$ and $\Psi^{\nu\mu}$ are given by

$\displaystyle\Psi^{\nu}$ $\displaystyle=$
$\displaystyle\Psi_{\mu}g^{\mu\nu},$ $\displaystyle\Psi^{0}$ $\displaystyle=$
$\displaystyle\frac{D\Psi_{0}-F\Psi_{3}}{AD+F^{2}},~{}~{}~{}\Psi^{1}=\frac{\Psi_{1}}{B},~{}~{}~{}\Psi^{2}=\frac{\Psi_{2}}{C},~{}~{}~{}\Psi^{3}=\frac{A\Psi_{3}-F\Psi_{0}}{AD+F^{2}},$
$\displaystyle\Psi^{pq}$ $\displaystyle=$
$\displaystyle\Psi_{\mu\nu}g^{p\mu}g^{q\nu}$ $\displaystyle\Psi^{01}$
$\displaystyle=$
$\displaystyle\frac{D\Psi_{01}-F\Psi_{31}}{B(AD+F^{2})},~{}~{}~{}\Psi^{02}=\frac{D\Psi_{02}-F\Psi_{32}}{C(AD+F^{2})},~{}~{}~{}\Psi^{03}=\frac{-\Psi_{03}}{AD+F^{2}},~{}~{}~{}\Psi^{12}=\frac{\Psi_{12}}{BC},$
$\displaystyle\Psi^{13}$ $\displaystyle=$
$\displaystyle\frac{A\Psi_{13}-F\Psi_{10}}{B(AD+F^{2})},~{}~{}~{}\Psi^{23}=\frac{A\Psi_{23}-F\Psi_{20}}{C(AD+F^{2})}.$

Applying the WKB approximation,

$\psi_{\nu}=k_{\nu}\exp[\frac{i}{\hbar}N_{0}(t,r,\theta,\phi)+\Sigma\hbar^{n}N_{n}(t,r,\theta,\phi)],$
(3)

The term $\hbar$ is assumed only for the lowest ($1^{st}$) order and the
higher (higher order$>1$) dictate contributions are neglected in the
Lagrangian gravity Eq. (2), we get set of field equations. The variables
separation technique in 29 , we can take

$N_{0}=-E_{0}t+W(r)+V(\theta)+J\phi,$ (4)

where $E_{0}=E-\omega\Omega$, $\omega$ and $E$ represent particle’s angular
momentum and energy respectively. The set of field equations, we can get a
matrix field equation

$K(k_{0},k_{1},k_{2},k_{3})^{T}=0,$

which implies the ”$K$” is label as $4\times 4$ matrix, whose elements are
given as follows:

$\displaystyle K_{00}$ $\displaystyle=$
$\displaystyle\frac{D}{B}\dot{W}^{2}W_{1}-\frac{D}{C}J^{2}J_{1}-\dot{V}^{2}V_{1}-eA_{3}[2\dot{V}+2\beta\dot{V}^{3}+eA_{3}+\beta
eA_{3}\dot{V}^{2}]-m^{2}D,$ $\displaystyle K_{01}$ $\displaystyle=$
$\displaystyle-\frac{D}{B}\dot{W}E_{0}[1+\beta E_{0}^{2}-eA_{0}E_{0}-\beta
eA_{0}E_{0}]-\frac{F}{B}\dot{W}[\dot{V}+\beta\dot{V}^{3}+eA_{3}+\beta
eA_{3}\dot{V}^{2}],$ $\displaystyle K_{02}$ $\displaystyle=$
$\displaystyle-\frac{D}{C}J[E_{0}-\beta E_{0}^{{}^{3}}+eA_{0}+\beta
eA_{0}E_{0}]-\frac{F}{C}J[\dot{V}+\beta\dot{V}^{3}+eA_{3}+\beta\dot{V}^{2}],$
$\displaystyle K_{03}$ $\displaystyle=$
$\displaystyle\frac{F}{B}\dot{W}^{2}W_{1}+\frac{F}{C}J^{2}J_{1}-E_{0}E_{1}\dot{V}+eA_{0}\dot{V}E_{1}-eA_{3}E_{0}E_{1}+e^{2}A_{0}A_{3}E_{1}+m^{2}F,$
$\displaystyle K_{10}$ $\displaystyle=$
$\displaystyle-D\dot{W}E_{0}W_{1}-F\dot{V}\dot{W}W_{1}+DeA_{0}\dot{W}W_{1}-FeA_{3}\dot{W}W_{1},$
$\displaystyle K_{11}$ $\displaystyle=$ $\displaystyle-
DE_{0}^{2}E_{1}-DeA_{0}E_{0}E_{1}-F\dot{V}E_{0}V_{1}-FeA_{3}E_{0}V_{1}-\frac{F_{0}}{C}J^{2}J_{1}-A\dot{V}^{2}V_{1}$
$\displaystyle-
A_{0}eA_{3}\dot{V}V_{1}-F\dot{V}E_{0}E_{1}+FeA_{0}\dot{V}[2+\beta
E_{0}^{2}+\beta\dot{V}^{2}]+DeA_{0}E_{0}E_{1}-De^{2}A^{2}_{0}E_{1}$
$\displaystyle+Fe^{2}A_{0}A_{3}V_{1}-AeA_{3}\dot{V}V_{1}-FeA_{3}E_{0}E_{1}+Fe^{2}A_{0}A_{3}E_{1}-m^{2}F_{0}$
$\displaystyle K_{12}$ $\displaystyle=$
$\displaystyle\frac{F_{0}}{C}[\dot{W}J+\beta\dot{W}^{3}J],$ $\displaystyle
K_{13}$ $\displaystyle=$ $\displaystyle
F\dot{W}E_{0}W_{1}+A\dot{W}\dot{V}W_{1}+AeA_{3}[E_{0}^{2}+\beta\dot{W}^{3}],$
$\displaystyle K_{20}$ $\displaystyle=$ $\displaystyle-
DE_{0}J_{1}-FJ\dot{V}J_{1}+DeA_{0}JJ_{1}-FeA_{3}JJ_{1},$ $\displaystyle
K_{21}$ $\displaystyle=$ $\displaystyle\frac{F_{0}}{B}\dot{W}JJ_{1},$
$\displaystyle K_{22}$ $\displaystyle=$ $\displaystyle-
DE_{0}^{2}E_{1}+DeA_{0}E_{0}E_{1}-F\dot{V}E_{0}V_{1}-FeA_{3}E_{0}J_{1}-\frac{F_{0}}{B}[\dot{W}^{2}-\beta\dot{W}^{4}]$
$\displaystyle-\dot{V}^{2}AV_{1}-AeA_{3}\dot{V}V_{1}-F\dot{V}E_{0}E_{1}+FeA_{0}\dot{V}[1+E_{1}+\beta\dot{V}^{2}]$
$\displaystyle+DeA_{0}E_{0}E_{1}-De^{2}A^{2}_{0}E_{1}+Fe^{2}A_{0}A_{3}V_{1}-AeA_{3}\dot{V}V_{1}-FeA_{3}E_{0}E_{1}$
$\displaystyle+Fe^{2}A_{0}A_{3}E_{1}-Ae^{2}A^{2}_{3}V_{1}-m^{2}F_{0},$
$\displaystyle K_{23}$ $\displaystyle=$ $\displaystyle
FJE_{0}J_{1}+AJ\dot{V}J_{1}+AeA_{3}JJ_{1}-FeA_{0}JJ_{1},$ $\displaystyle
K_{30}$ $\displaystyle=$
$\displaystyle-\dot{V}E_{0}V_{1}-eA_{3}E_{0}V_{1}+\frac{F}{B}\dot{W}^{2}W_{1}+\frac{F}{C}J^{2}J_{1}+eA_{0}\dot{V}V_{1}+e^{2}A_{0}A_{3}V_{1}+m^{2}F,$
$\displaystyle K_{31}$ $\displaystyle=$
$\displaystyle\frac{A}{B}\dot{V}\dot{W}V_{1}+\frac{A}{B}eA_{3}\dot{W}V_{1}+\frac{F}{B}\dot{W}[E_{0}+\beta
E_{0}^{3}-eA_{0}-eA_{0}\beta E_{0}^{2}],$ $\displaystyle K_{32}$
$\displaystyle=$
$\displaystyle\frac{A}{C}J\dot{V}V_{1}+\frac{A}{C}eA_{3}JV_{1}+\frac{F}{C}J[E_{0}+\beta
E_{0}^{3}+eA_{0}+eA_{0}\beta E_{0}^{2}],$ $\displaystyle K_{33}$
$\displaystyle=$ $\displaystyle-
E_{0}^{2}E_{1}+eA_{0}E_{0}E_{1}-\frac{A}{B}\dot{W}^{2}W_{1}-\frac{A}{C}J^{2}J_{1}+eA_{0}$
$\displaystyle[E_{0}+\beta E_{0}^{3}-eA_{0}-\beta eA_{0}]-m^{2}A,$

where $\dot{W}=\frac{\partial W}{\partial r}$, $\dot{V}=\frac{\partial
V}{\partial\theta}$, $E_{1}=1+\beta E_{0}^{2}$, $J_{1}=1+\beta J^{2}$,
$W_{1}=1+\beta\dot{W}^{2}$, $V_{1}=1+\beta\dot{V}^{2}$ and $F_{0}=AD+F^{2}$
and the imaginary parts of $imW_{-}$ and $imW_{+}$ yield

$imW_{+}=-imW_{-}=\pm\int\sqrt{\frac{(E_{0}-eA_{0})^{2}+X_{1}[1+\beta\frac{X_{2}}{X_{1}}]}{-AB^{-1}}}dr$
(5)

Here, $W_{-}$ and $W_{+}$ represent the solution of absorbing and radiating
charged boson particles action respectively and $X_{1}=m^{2}A+AC^{-1}J^{2}$
and
$X_{2}=E_{0}^{4}-eA_{0}E_{0}^{3}+\frac{A}{B}\dot{W}^{4}+\frac{A}{C}J^{4}-eA_{0}E_{0}^{3}+e^{2}A^{2}_{0}$.
In the series of Taylor’s, we are expanding the functions $A(r)$ and $B(r)$
near the horizon, we get

$A(r_{+})\approx\acute{A}(r)(r-r_{+}),~{}~{}B(r_{+})\approx\acute{B}(r)(r-r_{+})$
(6)

Since applying Eqs. (5) and (6) and integrate the around the pole, we get

$imW_{+}=+\pi\iota\frac{(E_{0}-eA_{0})^{2}}{2\kappa(r_{+})}(1+\Xi\beta),$ (7)

where $\Xi$ is arbitrary parameter and $\kappa(r_{+})$ is a surface gravity of
BH. The surface gravity of BH defines as;

$\kappa(r_{+})=\frac{2M(M-r_{+})}{(a^{2}+a^{2}cos^{2}\theta)^{2}}$ (8)

Here, the tunneling probability $\Gamma$ forbidden trajectories of the
classically of the $s$-waves coming from inside to outside from the BH
horizon. Applying the WKB approximation, by the terms of classical action
$S_{0}$ of charged boson particles tunneling across the BH horizon as
trajectories up to leading order is,

$\displaystyle\Gamma$ $\displaystyle=$
$\displaystyle\frac{\Gamma_{emission}}{\Gamma_{absorption}}=\frac{\exp(-2imW_{+}-2imV)}{\exp(-2imW_{-}-2imV)}=e^{-4imW_{+}}$
(9) $\displaystyle=$ $\displaystyle\exp\left[-2\pi\frac{(E-\omega\Omega-
eA_{0})}{\kappa(r_{+})}(1+\Xi\beta)\right].$

Thus, for computing the Hawking temperature and we expand the action in terms
of particles energy $E$. Hence, the Hawking temperature is calculated at
linear order given by

$T\cong\frac{M(r_{+}-M)}{\pi(a^{2}+a^{2}cos^{2}\theta)^{2}}[1-\Xi\beta].$ (10)

The charged boson massive particles case is similar as for massless case for
Rastall gravity Kerr-Newman BH surround by perfect fluid matter from the above
Hawking temperature. Moreover, the temperature for spin-up particles are
similar $(-W_{+}=W_{-})$ as for spin-down case with the change $r_{+}$ into
$r_{-}$. We observe that both spin-down and spin-up charged boson particles
are radiate with like rate the temperature is assumed for this case.

##  III Graphical Temperature Analysis

In this section, we analyze graphical behavior of thermodynamical quantities
in Rastall gravity Kerr-Newman BH surrounded by perfect fluid matter. For this
purpose, we take a parameter $(\Xi=1)$ and observe the effects of the
thermodynamical quantities (quantum gravity and rotation parameter) on BH
radiation. The BH remains unstable if temperature increases or decreases. From
Fig. 1, the Hawking temperature is decreasing due to the more quantum gravity.
The stability of the BH depends on quantum gravity and also shows that BH
remains stable unless the quantum gravity parameter is zero $(\beta=0)$. Also,
temperature decreases in the approximation range $(0<r_{+}<5)$.
If the quantum gravity is assumed to be non-zero, then Hawking temperature
attains the similar value in figure 2. The Hawking temperature is decrease to
rotation parameter in the values 0.1-0.3. The gravity parameter is fixed and
varying the rotation parameter then we concluded that rotation parameter is
inversely proportional to Hawking temperature and stability of BH. The
stability of BH remains constant in very small range, i.e., approximation
$0<r_{+}<0.5$ and then temperature decreases with a wide variation. The
rotation parameter is inversely proportional to the temperature of BH by
keeping the quantum gravity constant in Fig. 2 for 2D. The temperature is
behaving constant when $r_{+}$ is very small, but at a certain value of
rotation parameter, it attain a constantly decrease value in Fig. 2 for 2D.
Hence our BH is stable when $r_{+}$ is very small range otherwise is unstable.

We concluded that BH radiation depends on BH mass, rotation parameter and
gravity parameter but keeping constant $cos^{2}\theta=1$.

Figure 1: $T_{H}$ versus $r_{+}$ for $M>0,~{}a=10,~{}~{}cos^{2}\theta=1$ and
$\Xi=1.$

Figure 2: $T_{H}$ versus $r_{+}$ for $M>0,~{}~{}\beta>0,~{}~{}cos^{2}\theta=1$
and $\Xi=1.$

##  IV Conclusions

In this paper, we have observed radiation spectrum by considering Hawking
temperature for Rastall gravity Kerr-Newman BH is surrounded by perfect fluid
matter. For this purpose, we have utilized the WKB approximation to the
Hamilton-Jacobi ansatz for massive bosons charged spin-1 particles. The
calculation yields, the Hawking temperature reliable with BH universality. In
our investigation, we have changed the Lagrangian field equation in curved
space time by effect of incorporate GUP quantum gravity. We have computed the
tunneling probability at event horizon as well as the equating Hawking
temperature. We have observed temperature in detail the effect of quantum
gravity and rotation parameters in graphically.

The paper have studied graphical investigation of the result and the BH will
emit away all types of particles as the black-body radiation and the BH lose
completely its information. In adjustment to break this problem, the after-
effect rotation and gravity to be modified. We assume into account the law of
conservation of momentum and energy. In our observation, the back reaction of
the radiated boson particles of the self-gravitational interaction and BH
geometry are reasonably ignored, the calculated Hawking temperatures depend on
quantum gravity and the leading terms of metric. The modified form of
temperature depends on BH mass, rotation parameters and quantum gravity
parameter. When the gravity parameter $(\beta=0)$ influences are ignored, we
have obtained the temperature of Rastall gravity Kerr-Newman BH. So, the
quantum gravity effect minimum then stability of BH are increased.

The quantum gravity and rotation are present in the BH, which shows that the
instability of BH. The BH is completely stable when rotation and quantum
gravity do not appear.
Acknowledgments
The work of KB was supported in part by the JSPS KAKENHI Grant Number
JP21K03547.

## Appendix

In the WKB approximation the term $\hbar$ is considered just for the $1^{st}$
order and the higher dictate contributions are neglected in the Lagrangian Eq.
(2), we get following field equations;

$\displaystyle\frac{D}{B}[k_{1}(\partial_{0}N_{0})(\partial_{1}N_{0})+\beta
k_{1}(\partial_{0}N_{0})^{3}(\partial_{1}N_{0})+k_{0}(\partial_{1}N_{0})^{2}-\beta
k_{0}(\partial_{1}N_{0})^{4}+eA_{0}k_{1}(\partial_{0}N_{0})^{2}(\partial_{1}N_{0})+\beta
eA_{0}k_{1}(\partial_{0}N_{0})^{2}(\partial_{1}N_{0})]$
$\displaystyle-\frac{F}{B}[k_{1}(\partial_{1}N_{0})(\partial_{3}N_{0})+\beta
k_{1}(\partial_{1}N_{0})(\partial_{3}N_{0})^{3}-k_{3}(\partial_{1}N_{0})^{2}-\beta
k_{3}(\partial_{1}N_{0})^{4}+eA_{3}k_{1}(\partial_{1}N_{0})+\beta
eA_{3}k_{1}(\partial_{1}N_{0})(\partial_{3}N_{0})^{2}]$
$\displaystyle+\frac{D}{C}[k_{2}(\partial_{0}N_{0})(\partial_{2}N_{0})+\beta
k_{2}(\partial_{0}N_{0})^{3}(\partial_{2}N_{0})-k_{o}(\partial_{2}N_{0})^{2}-\beta
k_{0}(\partial_{2}N_{0})^{4}+eA_{0}k_{2}(\partial_{2}N_{0})+\beta
eA_{0}k_{2}(\partial_{0}N_{0})^{2}(\partial_{2}N_{0})]$
$\displaystyle-\frac{F}{C}[k_{2}(\partial_{2}N_{0})(\partial_{3}N_{0})+\beta
k_{2}(\partial_{2}N_{0})(\partial_{3}N_{0})^{3}-k_{3}(\partial_{2}N_{0})^{2}-\beta
k_{3}(\partial_{2}N_{0})^{4}+eA_{3}k_{2}(\partial_{2}N_{0})+\beta
eA_{3}k_{2}(\partial_{2}N_{0})(\partial_{3}N_{0})^{2}]$
$\displaystyle+[k_{3}(\partial_{0}N_{0})+\beta
k_{3}(\partial_{0}N_{0})^{3}(\partial_{3}N_{0})-k_{0}(\partial_{3}N_{0})^{2}-\beta
k_{0}(\partial_{3}N_{0})^{4}+eA_{0}k_{3}(\partial_{3}N_{0})+\beta
eA_{0}k_{3}(\partial_{0}N_{0})^{2}(\partial_{3}N_{0})-eA_{3}k_{0}(\partial_{3}N_{0})$
$\displaystyle-\beta
eA_{3}k_{0}(\partial_{3}N_{0})^{3}]+eA_{3}[k_{3}(\partial_{0}N_{0})+\beta
k_{3}(\partial_{0}N_{0})^{3}-k_{0}(\partial_{3}N_{0})-\beta
k_{0}(\partial_{3}N_{0})^{3}+eA_{0}k_{3}+\beta
eA_{0}k_{3}(\partial_{0}N_{0})^{2}-eA_{3}k_{0}$ $\displaystyle-\beta
eA_{3}k_{0}(\partial_{3}N_{0})^{2}]=0,$ (11) $\displaystyle
D[k_{0}(\partial_{0}N_{0})(\partial_{1}N_{0})+\beta
k_{0}(\partial_{0}N_{0})(\partial_{1}N_{0})^{3}-k_{1}(\partial_{0}N_{0})^{2}-\beta
k_{1}(\partial_{0}N_{0})^{4}-eA_{0}k_{1}(\partial_{0}N_{0})-\beta
eA_{0}k_{1}(\partial_{0}N_{0})^{3}]-F[k_{3}(\partial_{0}N_{0})$
$\displaystyle(\partial_{1}N_{0})+\beta
k_{3}(\partial_{0}N_{0})(\partial_{1}N_{0})^{3}-k_{1}(\partial_{0}N_{0})(\partial_{3}N_{0})-\beta
k_{1}(\partial_{0}N_{0})(\partial_{3}N_{0})^{3}-eA_{3}k_{1}(\partial_{0}N_{0})-\beta
eA_{3}k_{1}(\partial_{0}N_{0})(\partial_{3}N_{0})^{2}]$
$\displaystyle+(AD-F^{2})C^{-1}[k_{2}(\partial_{1}N_{0})(\partial_{2}N_{0})+\beta
k_{2}(\partial_{1}N_{0})^{3}(\partial_{2}N_{0})-k_{1}(\partial_{2}N_{0})^{2}-\beta
k_{1}(\partial_{2}N_{0})^{4}]+A[k_{3}(\partial_{1}N_{0})(\partial_{3}N_{0})+\beta
k_{3}$
$\displaystyle(\partial_{1}N_{0})^{3}(\partial_{3}N_{0})-k_{1}(\partial_{3}N_{0})^{2}-\beta
k_{1}(\partial_{3}N_{0})^{4}-eA_{3}k_{1}(\partial_{3}N_{0})-\beta
eA_{3}k_{1}(\partial_{3}N_{0})^{3}]-F[k_{0}(\partial_{1}N_{0})(\partial_{3}N_{0})+\beta
k_{0}(\partial_{1}N_{0})^{3}$
$\displaystyle(\partial_{3}N_{0})-k_{1}(\partial_{0}N_{0})(\partial_{3}N_{0})-\beta
k_{1}(\partial_{0}N_{0})^{3}(\partial_{3}N_{0})-eA_{0}k_{1}(\partial_{3}N_{0})-\beta
eA_{0}k_{1}(\partial_{0}N_{0})^{2}(\partial_{3}N_{0})]+DeA_{0}[k_{0}(\partial_{1}N_{0})+\beta
k_{0}$ $\displaystyle(\partial_{1}N_{0})^{3}-k_{1}(\partial_{0}N_{0})-\beta
k_{1}(\partial_{0}N_{0})^{3}-eA_{0}k_{1}-\beta
eA_{0}k_{1}(\partial_{0}N_{0})^{2}]-FeA_{0}[k_{3}(\partial_{1}N_{0})+\beta
k_{3}(\partial_{1}N_{0})^{3}-k_{1}(\partial_{3}N_{0})-\beta k_{1}$
$\displaystyle(\partial_{3}N_{0})^{3}-eA_{3}k_{1}-\beta
eA_{3}k_{1}(\partial_{3}N_{0})^{2}]+eAA_{3}[k_{3}(\partial_{0}N_{0})^{2}+\beta
k_{3}(\partial_{1}N_{0})^{3}-k_{1}(\partial_{3}N_{0})-\beta
k_{1}(\partial_{3}N_{0})^{3}-eA_{3}k_{1}-\beta eA_{3}$ $\displaystyle
k_{1}(\partial_{3}N_{0})^{2}]-eFA_{3}[k_{0}(\partial_{1}N_{0})+\beta
k_{0}(\partial_{1}N_{0})^{3}-k_{1}(\partial_{0}N_{0})-\beta
k_{1}(\partial_{0}N_{0})^{3}-eA_{0}k_{1}-\beta
eA_{0}k_{1}(\partial_{0}N_{0})^{2}]-m^{2}k_{1}$ $\displaystyle(AD-F^{2})=0$
(12) $\displaystyle D[k_{0}(\partial_{0}N_{0})(\partial_{2}N_{0})+\beta
k_{0}(\partial_{0}N_{0})(\partial_{2}N_{0})^{3}-k_{2}(\partial_{0}N_{0})^{2}-\beta
k_{2}(\partial_{0}N_{0})^{4}-eA_{0}k_{2}(\partial_{0}N_{0})-\beta
eA_{0}k_{2}(\partial_{0}N_{0})^{3}]-F[k_{3}(\partial_{0}N_{0})$
$\displaystyle(\partial_{2}N_{0})+\beta
k_{3}(\partial_{0}N_{0})(\partial_{2}N_{0})^{3}-k_{2}(\partial_{0}N_{0})(\partial_{3}N_{0})-\beta
k_{2}(\partial_{0}N_{0})(\partial_{3}N_{0})^{3}-eA_{3}k_{2}(\partial_{0}N_{0})-\beta
eA_{3}k_{2}(\partial_{0}N_{0})(\partial_{3}N_{0})^{2}]$
$\displaystyle+(AD-F^{2})B^{-1}[k_{1}(\partial_{1}N_{0})(\partial_{2}N_{0})+\beta
k_{1}(\partial_{1}N_{0})(\partial_{2}N_{0})^{3}-k_{2}(\partial_{1}N_{0})^{2}-\beta
k_{2}(\partial_{1}N_{0})^{4}]+A[k_{3}(\partial_{2}N_{0})(\partial_{3}N_{0})+\beta
k_{3}$
$\displaystyle(\partial_{2}N_{0})^{3}(\partial_{3}N_{0})-k_{2}(\partial_{3}N_{0})^{2}-\beta
k_{3}(\partial_{3}N_{0})^{4}-eA_{3}k_{2}(\partial_{3}N_{0})-\beta
eA_{3}k_{2}(\partial_{3}N_{0})^{3}]-F[k_{0}(\partial_{2}N_{0})(\partial_{3}N_{0})+\beta
k_{0}(\partial_{2}N_{0})^{3}$
$\displaystyle(\partial_{3}N_{0})-eA_{0}k_{2}(\partial_{3}N_{0})-\beta
eA_{0}k_{2}(\partial_{0}N_{0})^{2}(\partial_{3}N_{0})-k_{2}(\partial_{0}N_{0})(\partial_{3}N_{0})-\beta
k_{2}(\partial_{0}N_{0})^{3}(\partial_{3}N_{0})]+DeA_{0}[k_{0}(\partial_{2}N_{0})+\beta
k_{0}$ $\displaystyle(\partial_{2}N_{0})^{3}-k_{2}(\partial_{0}N_{0})-\beta
k_{2}(\partial_{0}N_{0})^{3}-eA_{0}k_{2}-\beta
eA_{0}k_{2}(\partial_{0}N_{0})^{2}]-FeA_{0}[k_{3}(\partial_{2}N_{0})+\beta
k_{3}(\partial_{2}N_{0})^{3}-k_{2}(\partial_{3}N_{0})-\beta k_{2}$
$\displaystyle(\partial_{3}N_{0})^{3}-eA_{3}k_{2}-\beta
eA_{3}k_{2}(\partial_{3}N_{0})^{2}]+eAA_{3}[k_{3}(\partial_{2}N_{0})+\beta
k_{3}(\partial_{2}N_{0})^{3}-k_{2}(\partial_{3}N_{0})-\beta
k_{2}(\partial_{3}N_{0})^{3}$ $\displaystyle-eA_{3}k_{2}-\beta
eA_{3}k_{2}(\partial_{3}N_{0})^{2}]-eFA_{3}[k_{0}(\partial_{2}N_{0})+\beta
k_{0}(\partial_{2}N_{0})^{3}-k_{2}(\partial_{0}N_{0})-\beta
k_{2}(\partial_{0}N_{0})^{3}-eA_{0}k_{2}-\beta
eA_{0}k_{2}(\partial_{0}N_{0})^{2}]$ $\displaystyle-m^{2}k_{2}(AD-F^{2})=0$
(13) $\displaystyle[k_{0}(\partial_{0}N_{0})(\partial_{3}N_{0})+\beta
k_{0}(\partial_{0}N_{0})(\partial_{3}N_{0})^{3}+k_{3}(\partial_{0}N_{0})^{2}-\beta
k_{3}(\partial_{0}N_{0})^{4}+eA_{3}k_{0}(\partial_{0}N_{0})+\beta
eA_{3}k_{0}(\partial_{0}N_{0})(\partial_{3}N_{0})^{2}$ $\displaystyle-
eA_{0}k_{3}(\partial_{0}N_{0})-\beta
eA_{0}k_{3}(\partial_{0}N_{0})^{3}]+\frac{A}{B}[k_{1}(\partial_{1}N_{0})(\partial_{3}N_{0})+\beta
k_{1}(\partial_{1}N_{0})(\partial_{3}N_{0})^{3}-k_{3}(\partial_{1}N_{0})^{2}-\beta
k_{3}(\partial_{1}N_{0})^{4}+$ $\displaystyle
eA_{3}k_{1}(\partial_{1}N_{0})+\beta
eA_{3}k_{1}(\partial_{1}N_{0})(\partial_{3}N_{0})^{2}]-\frac{F}{B}[k_{1}(\partial_{0}N_{0})(\partial_{1}N_{0})-\beta
k_{1}(\partial_{0}N_{0})^{3}(\partial_{1}N_{0})-k_{o}(\partial_{1}N_{0})^{2}-\beta
k_{0}(\partial_{1}N_{0})^{4}$
$\displaystyle+eA_{0}k_{1}(\partial_{1}N_{0})+\beta
eA_{0}k_{1}(\partial_{0}N_{0})^{2}(\partial_{1}N_{0})]+\frac{A}{C}[k_{2}(\partial_{2}N_{0})(\partial_{3}N_{0})+\beta
k_{2}(\partial_{2}N_{0})(\partial_{3}N_{0})^{3}-k_{3}(\partial_{2}N_{0})^{2}-\beta
k_{3}(\partial_{2}N_{0})^{4}$
$\displaystyle+eA_{3}k_{2}(\partial_{2}N_{0})+\beta
eA_{3}k_{2}(\partial_{2}N_{0})(\partial_{3}N_{0})^{2}]-\frac{F}{C}[k_{2}(\partial_{0}N_{0})(\partial_{2}N_{0})+\beta
k_{2}(\partial_{0}N_{0})^{2}(\partial_{2}N_{0})-k_{0}(\partial_{2}N_{0})^{2}-\beta
k_{0}(\partial_{2}N_{0})^{4}$
$\displaystyle+eA_{0}k_{2}(\partial_{2}N_{0})+\beta
eA_{0}k_{2}(\partial_{0}N_{0})^{2}(\partial_{2}N_{0})]+eA_{0}[k_{0}(\partial_{3}N_{0})+\beta
k_{0}(\partial_{3}N_{0})^{3}-k_{3}(\partial_{0}N_{0})-\beta
k_{3}(\partial_{0}N_{0})^{3}+eA_{3}k_{0}+$ $\displaystyle\beta
eA_{3}k_{0}(\partial_{3}N_{0})^{2}-eA_{0}k_{3}-\beta
eA_{0}k_{3}(\partial_{0}N_{0})^{2}]+m^{2}[Fk_{0}-Ak_{3}]=0.$ (14)

## References

  * (1) S. W. Hawking, Nature 248, 30(1974); ibid. Commun. Math. Phys. 43, 199(1975); ibid. Phys. Rev. D 13,191(1976).
  * (2) M. Sharif and W. Javed, Eur. Phys. J. C 72, (2012); ibid. Can. J. Phys., 91, 43(2013).
  * (3) K. Jan and H. Gohar.: Astrophys. Space Sci. 350, 279(2014).
  * (4) S. I. Kurglow.: Int. J. Mod. Phys. A 29, 1450118(2014).
  * (5) Z. Feng, et al.: Astrophys. Space Sci. 48, 359(2015).
  * (6) G. R. Chen and Y. C. Huang.: Int. J. Mod. Phys. A 30, 1550083(2015).
  * (7) M. A. Anacleto, et al.: Annal. Phys. 362, 436(2015).
  * (8) G. Li and X. Zu.: J. of Applied Mathematics and Physics 3, 134(2015).
  * (9) X. Q. Li and G. R. Chen.: Phys. Lett. B 751, 34(2015).
  * (10) Z. Feng, Y. Chen and X. Zu.: Astrophys. Space Sci. 48, 359(2015).
  * (11) T. I. Singh, I. A. Meitei and K. Y. Singh.: Astrophys. Space Sci. 361, 103(2016).
  * (12) H. Gursel and I. Sakalli.: Can. J. Phys. 94, 147(2016).
  * (13) A. Övgün and K. Jusufi.: Eur. Phys. J. Plus 131, 177(2016).
  * (14) R. Li and J. K. Zhao.: Commun. Theor. Phys. 65,469(2016).
  * (15) X. Q. Li.: Phys. Lett. B 763, 80(2016).
  * (16) K. Jusufi and A. Övgün.: Astrophys. Space Sci. 361, 207(2016).
  * (17) J. Sadeghi and V. R. Shajiee.: Eur. Phys. J. Plus 132, 132(2017).
  * (18) K. Jusufi and A. Övgün.: Int. J. Theor. Phys. 56, 1725(2017).
  * (19) W. Javed, G. Abbas and R. Ali.: Eur. Phys. J. C 77, 296(2017).
  * (20) W. Javed, R. Ali, R. Babar, A. Övgün, Eur. Phys. J. Plus 134, 511(2019).
  * (21) W. Javed, R. Ali and G. Abbas.: Can. J. Phys. 97, 176(2019).
  * (22) I. Sakalli and A. Övgün., Gen. Rel. Grav., 48, 1(2016).
  * (23) I. Sakalli and A. Övgün, Astrophys. Space Sci., 359, 32(2015).
  * (24) X. M. Kuang, J. Saavedra and A. Övgün, Eur. Phys. J. C 77, 613(2017).
  * (25) X. M. Kuang, B. Liu and A. Övgün, Eur. Phys. J. C 78, 840(2018).
  * (26) W. Javed, R. Babar and A. Övgün, Mod. Phys. Lett. A 34, 1950057(2019).
  * (27) P. A. Gonzalez, A. Övgün, J. Saavedra et al, Gen. Rel. Grav. 50, 62(2018).
  * (28) D. Chen, Eur. Phys. J. C 74, 2687(2014).
  * (29) G. Gecim and Y. Sucu, Phys. Lett. B 773, 391(2017).
  * (30) M. Hossain Ali, Class. Quant. Grav., 24, 5849(2007).
  * (31) M. Hossain Ali, Gen. Rel. Grav., 36, 1171(2004).
  * (32) R. Ali, K. Bamba, S. A. A. Shah, Symmetry. 631, 11(2019).
  * (33) R. Ali, et al, Symmetry. 1165, 12(2020).
  * (34) E. T. Akhmedov, V. Akhmedova and D. Singleton, Phys. Lett. B 642, 124(2006).
  * (35) R. Ali, R. Babar, M. Asgher, and S.A.A. Shah, Annals of Physics, 432, 168572(2021).
  * (36) V. Akhmedova, et al, Phys. Lett. B 666, 269(2008).
  * (37) V. Akhmedova, et al, Phys. Lett. B 673, 227(2009).
  * (38) G. Gecim and Y. Sucu, Gen. Rel. Grav., 50, 152(2018).
  * (39) G. R. Chen, S. Zhou and Y. C. Huang, Int. J. Mod. Phys. D 24, 1550005(2014).
  * (40) X. X. Zeng and S. Z. Yang, Chin. Phys. B 18, 462(2009).
  * (41) S. K. Modak, Phys. Lett. B 671, 167(2009).
  * (42) R. Li and J. R. Ren, Phys. Lett. B 661, 370(2008).
  * (43) Q. Q. Jiang, S. Q. Wu, and X. Cai, Phys. Lett. B 651, 58(2007).
  * (44) X. He and W. Liu, Phys. Lett. B 653, 330(2007).
  * (45) H. L. Li, et al. Phys. Lett. B 641, 139(2006).
  * (46) Z. Feng, Y. Chen, and X. Zu, Astrophys. Space Sci., 48, 359(2015).
  * (47) K. Jusufi and A. Övgün, Astrophys. Space Sci., 361, 207(2016).
  * (48) R. Ali, M. Asgher and M. F. Malik, Mod. Phys. Lett. A 35, 2050225(2020).
  * (49) R. Ali, et al. Int. J. Mod. Phys. D 30, 2150002(2021).
  * (50) R. Ali, R. Babar, P.K. Sahoo, Physics of the Dark Universe 35, 100948(2022).
  * (51) R. Ali and M. Asgher, New Astronomy 93, 101759(2022).
  * (52) R. Ali, et al. Int. J. Geom. Methods Mod. Phys. 19, 2250017(2022).
  * (53) R. Kumar and S. G. Ghosh, Eur. Phys. J. C 78, 750(2018.
  * (54) M. Visser.: Phys. Lett.B 782, 83(2018).
  * (55) Z. Xu, X. Hou, X. Gong and J. Wang, Eur. Phys. J. C 78, 513(2018).
  * (56) R. Kumar, et al. Physics of the Dark Universe 34, 100881(2021).
  * (57) S. G. Ghosh, Eur. Phys. J. C 76, 222(2016).
  * (58) W. Javed, R. Ali, R. Babar, A. Övgün, Chinese Physics C 44, No. 1, 015104(2020).