text:
\begin{document}

\begin{frontmatter}

\title{Enriched Gradient Recovery for Interface Solutions of the Poisson-Boltzmann Equation}

\author{George Borleske, Y. C. Zhou\fnref{myfootnote}}
\address{Department of Mathematics, Colorado State University, Fort Collins, CO 80523}
\fntext[myfootnote]{Corresponding author. Email: yzhou@math.colostate.edu}

\begin{abstract}
Accurate calculation of electrostatic potential and gradient on the molecular surface is highly desirable for the continuum and hybrid modeling of
large scale deformation of biomolecules in solvent. In this article a new numerical method is proposed to calculate these quantities on the dielectric
interface from the numerical solutions of the Poisson-Boltzmann equation. Our method reconstructs a potential field locally in the least square sense
on the polynomial basis enriched with Green's functions, the latter characterize the Coulomb potential induced by charges near the position of reconstruction.
This enrichment resembles the decomposition of electrostatic potential into singular Coulomb component and the regular reaction field in the Generalized Born methods.
Numerical experiments demonstrate that the enrichment recovery produces drastically more accurate and stable potential gradients on molecular surfaces
compared to classical recovery techniques.
\end{abstract}

\begin{keyword}
Biomolecular electrostatics; Poisson-Boltzmann equation; Numerical Solution; Interface methods; Gradient
recovery; High accuracy
\MSC[2014] 53B10 \sep 65D18 \sep 92C15
\end{keyword}

\end{frontmatter}

\section{Introduction}
The Poisson-Boltzmann theory has been widely accepted as a mean-field continuum approximation for electrostatic interactions
in solvated biomolecular systems \cite{Honig1995_science}. The Poisson-Boltzmann theory treats the solute bimolecules as
a singularly charged medium of low dielectric constant ($\epsilon_p = 1 \sim 2$) immersed in a high dielectric ($\epsilon_s = 75 \sim 80$) solvent with a
continuum charge distribution that models the dispersed mobile ions. The large contrast of dielectric constant on
the highly complicated biomolecular surfaces poses a significant computational challenge and affords delicate mathematical and
numerical treatments. The importance of Poisson-Boltzmann theory in biochemistry and biophysics has motivated
extensive mathematical and computational investigations, see, for instance \cite{PBEQ,LuB2007e,RenP2012a,XieD2013a,ChenJ2018a,ZhongY2018a},
and references therein for the development of the subjects in the past decades and the latest overviews.

This paper is concerned with the accurate calculation of the gradient of the electrostatic potential $\nabla \phi$ near the molecular surface
from the numerical solutions of the Poisson-Boltzmann equation (PBE) for modeling solvated biomolecules. Traditionally the potential solutions
of the Poisson-Boltzmann equation are used merely for computing the solvation free energy (and other derived energetic quantities such as pKa value,
binding affinity, etc.), and therefore less attention had been paid on the calculation of the potential gradient. The success of Poisson-Boltzmann
theory in these energetic evaluations promotes the exploration of PBE based electrostatic force calculations for molecular dynamics (MD)
simulations \cite{MacKerellA1998a,LuQ2003a,Prabhu2004a,GengW2011a} and coarse-grained \cite{TangY2006a,MaL2009}
or hybrid models \cite{ChenX2010a,ZhouY2010b,ChengL2013a}. There are three types of electrostatic forces defined by using the free energy
functional derivatives and shape
derivatives \cite{GilsonM1993a,ZhouY2008c,LiB2011a,MikuckiM2014a},
including (i) the body force exerted at each charged atoms of the solution molecules given by $qE = q \nabla \phi$ where $q$ is the local charge density;
(ii) the dielectric boundary force exerted on the dielectric interface that is usually defined by the molecular surface;
and (iii) the ionic pressure, exerted also on the molecular surface.

Perpendicular to the dielectric interface, the dielectric boundary force
\begin{equation}
f_n =-\frac{\epsilon_s}{2} |\nabla \phi^s|^2 + \frac{\epsilon_p}{2} |\nabla \phi^p|^2 - \epsilon_p |\nabla_n \phi^p|^2 + \epsilon_p (\nabla_n \phi^s)
(\nabla_n \phi^p) \label{eqn:surface_force}
\end{equation}
poses significant challenges to the numerical computation because it is defined on the dielectric interface where one usually observes
the peaks of the numerical error in the solutions of the
Poisson-Boltzmann equation \cite{Stillinger1961a,Xiang1995a,Holst1995a,Holst2000a,BoschitschA2002a,Xie2007a,LinH2014a,Geng2017a}.
For boundary element methods, while $\nabla_n \phi^s$ can be obtained directly along with $\phi^s$ from well-conditioned boundary integral formulations,
computing $\nabla \phi^p$ involves additional integrations of the numerical solution of the potential with a supersingular kernel over the entire
surface \cite{LuB2006a,BajajC2011a,GengW2013a,ZhongY2018a}, usually causing a larger error. For finite difference interface methods, one usually extrapolates
the potential solution across the interface by using the interface conditions to compute potential gradients \cite{GengW2011a}. For interface finite elements
methods, gradient recovery techniques are recently developed so that potential gradient could be computed in a subdomain without using the solution values
across the internal interface \cite{GuoH2018a,GuoH2017a,GuoH2018b}.

We will develop an interface gradient recovery method to approximate the screened electrostatic potential from the numerical solution of the Poisson-Boltzmann
equation. For elliptic problems with smooth coefficients, gradient recovery techniques have been well established on structured or unstructed
grids \cite{ZienkiewiczO1992a,ZienkiewiczO1992b,BankR2003a,GuoH2015a,NagaA2005a,XuJ2003a,ZhangZ2005a}. By comparison, there is only a small handful of
numerical techniques concerning the gradient recovery from the primarily computed solutions of the elliptic problems with piecewise smooth
coefficients or singular sources. For 1-D interface problems, special interpolation formulas are constructed to recover flux with high order accuracy
from numerical solutions of linear and quadratic interface finite element methods \cite{ChouS2012a,ChouS2015a}. For 2-D elliptic interface problems on a
specially constructed body-fitted mesh, gradient of the linear finite element solution is shown to be superclose to the gradient of the linear
interpolation of the exact solution \cite{WeiH2014a}. By introducing the jump in the normal derivative as an augmented variable, Li {\it et. al.}
showed that a second-order convergence can be obtained for the solution and its gradient \cite{LiZ2017a}. Recently developed gradient recovery
methods for interface problems are based on the second-order least square reconstruction of the solutions in an individual subdomain \cite{GuoH2018a}.
If the mesh is not conforming to the interface, a practice common in many immersed finite element methods \cite{LiZ1998a,LiuX2000a,HouS2013a,JiH2014a},
the solutions on the interface are first obtained on approperiately subdivided interface elements then supplied to the least square
reconstruction \cite{GuoH2017a}. When Nitsche's method \cite{HanshoP2005a,AnnavarapuC2012a} is used for solving the elliptic interface problems,
this subdivision of the interface element is not necessary, and one can use the solution on the extended fictitious subdomain to carry out the
least square reconstruction \cite{GuoH2018b}. Nevertheless, as we shall demonstrate below, these interface gradient recovery methods are not able
to accurately compute the interface potential gradient on a spatial (domain or surface) discretization that are affordable for realistic
biomolecular simulations.

Our method is based on the least square reconstruction of the numerical solution, with the classical polynomial basis supplemented with
$\{1/r_i\}$ where $r_i$ is the distance between a mesh node and the selected charges in the solute biomolecule. These additional functions have
been indicated in the Green's function for the Poisson equation and particularly in the generalized Born (GB) models for the biomolecular
electrostatics. In GB models, the electrostatics potential at a charge $q_i$ inside the dielectric boundary is given by
\begin{equation}
\phi(x_i) = \sum_{j \ne i} \frac{q_i q_j}{\epsilon_p| x_i - x_j|} + \bfF(x_i,x_j)  \label{eqn:model_GB}
\end{equation}
where $q_i$ is the singular charge at spatial position $\bfx_i$, and $\bfF$ is the reaction field due to the polarization charges induced at
the dielectric boundary \cite{ImW2003a,FeigM2004c}. A large variety of GB models have been developed to implicitly account for the effects of
solvation through various parameterizations of effective Born radius, solvation accessible surfaces,  ionic screening, among others, see
\cite{OnufrievA2019a} for a latest review and references therein. In the solvent domain $\Omega_s$, although an analytical approximation similar
to (\ref{eqn:model_GB}) does not exist, the approximation
\begin{equation}
\phi(x) = \sum_{i} \frac{q_i}{\epsilon_s| x - x_i|} e^{-\kappa |x - x_i|} \label{eqn:PB_bc}
\end{equation}
is widely used to compute the boundary conditions for the Poisson-Boltzmann at a boundary point $\bfx$ far away from the
solute molecules,  where $\kappa$ is the Debye-Huckel inverse screening length. The approximations (\ref{eqn:model_GB})-(\ref{eqn:PB_bc}) motivate
us to choose functions of the form $1/r$ to enrich the classical polynomial basis function in the least square reconstruction.
Numerical experiments demonstrate that our enriched recovery technique is highly accurate and stable for complicated biomolecules.
By construction our recovery technique can be integrated with general interface numerical methods for the Poisson-Boltzmann equation,
regardless of whether the discretization mesh is interface conforming or not.

The rest of the article is organized as follows. In Section \ref{sect:PBE}, the Poisson-Boltzmann equation and
its numerical treatment are briefly reviewed, followed by the introduction of dielectric boundary force and its calculation.
We will present the our gradient recovery technique in Section \ref{sect:recover} and discuss its extension to enforcing the
interface conditions. Extensive numerical experiments will be conducted in Section \ref{sect:numeric} to verify the robustness
and accuracy of our methods on biomolecules of different complexity. The article concludes with a summary in
Section \ref{sect:summary}.

\section{Poisson-Boltzmann Equation and Dielectric Boundary Force} \label{subsect:PBE}
This section presents the energetic theory of biomolecular electrostatics. The electrostatic free energy introduced by Sharp, Honig, Gilson {\it et. al.}
\cite{GilsonM1993a,GilsonM1985a,SharpK1990a,SharpK1990b} provides a unified theory to connect the electrostatic solvation energy and the Poisson-Boltzmann
equation. The inconsistency between their definition of the electrostatic force and the sharp dielectric interface model was recently reconciled through
the introduction of a shape derivative of the solute biomolecules and the application of the Hadamard-Zol\'{e}sio structure
theorem \cite{LiB2011a,MikuckiM2014a}. We will review the regularization techniques for the numerical solution of the PBE, which also sheds light to
our choice of the enriching basis functions in the least square reconstruction for the gradient recovery.

\subsection{Biomolecular electrostatics theory and the Poisson-Boltzmann equation} \label{sect:PBE}
Consider $\Omega \in \bbR^3$ be a domain that encapsulates the solute biomolecules and the aqueous solution in which $M$ species of mobile ions disperse.
The occupation domains of solute and solvent are respectively denoted by $\Omega_p$ and $\Omega_s$, with distinct dielectric constants $\epsilon_p$
and $\epsilon_s$. Denote by $\Gamma$ the dielectric interface separating these two subdomains. Assume dispersion of mobile ions follows the Boltzmann distribution,
then the electrostatic free energy of this solvated system is given by
\begin{equation} \label{eqn:ele_eng}
G(\Gamma;\phi) = \int_{\Omega} \left[ -\frac{\epsilon}{2} | \nabla \phi|^2 + f \phi - \chi_s \beta^{-1} \sum_{j=1}^M c_j^0 (e^{-\beta q_j \phi} -1 ) \right] dx,
\end{equation}
where $c_j^0$ is the bulk concentration of $j^{th}$ species of ions, $f=\sum_{i=1}^N q_i$ is the collection of singular charges of biomolecule atoms, and
$\beta=1/(k_BT)$ is the inverse thermal energy at temperature $T$. The characteristic function $\chi_s=1$ in $\Omega_s$ and vanishes elsewhere.
Equation (\ref{eqn:ele_eng}) highlights the dependence of the electrostatic free energy on the
location of the dielectric interface $\Gamma$. Minimization of $G$ with respect to $\phi$ leads to the nonlinear Poisson-Boltzmann equation
\begin{equation} \label{eqn:PBE}
-\nabla \cdot( \epsilon \nabla \phi) + \chi_s \sum_{j=1}^M c_j^0 q_j  e^{-\beta q_j \phi} = \sum_{i=1}^N q_i \delta(x_i),
\end{equation}
with interface conditions on $\Gamma$:
\begin{equation} \label{eqn:PBE_bc}

\left[ \phi(x) \right ] = 0, \quad \left[ \epsilon \nabla_n \phi \right] = 0,
\end{equation}
where $[\cdot]$ denotes the jump of the enclosed function and $\nabla_n$ is the normal derivative pointing to $\Omega_s$.

To derive the dielectric boundary force we introduce a surface transformation $T_t$ corresponding to the slight change of
domain $\Omega_p$ due to the motion of interface $\Gamma$ in its normal direction. The dielectric boundary force is derived by
taking the variational derivative of the free energy $G$ with respect to  a velocity field $V$ normal to $\Gamma$. Associated with
this velocity fields is a mapping $T_t(X)$
from the original surface coordinate $X \in \Gamma_0$ (material position) to the deformed coordinate $x \in \Gamma_t$ (physical position):
\begin{equation} \label{eqn:deformation}
T_t(X) = x(t,X) = x(0,X) + t \left . \frac{\partial x(t,X)}{\partial t} \right |_{t=0} + \mathcal{O}(t^2).
\end{equation}
Introducing the Jacobian $J_s$ of the surface transformation $T_t$:
\begin{equation} \label{eqn:J_deform}
J_s = \det \left(\nabla T_t(X) \right) \cdot | \nabla T_t^{-T} \cdot n(X)|,
\end{equation}
one can show that
\begin{equation} \label{eqn:V}
\nabla \cdot V = \left. \frac{d J_s}{dt} \right |_{t=0},
\end{equation}
and the shape derivative of $G$
\begin{equation} \label{eqn:var_G}
\delta_{\Gamma} G(\Gamma; \phi) = \lim_{t \rightarrow 0} \frac{G(\Gamma_t;\phi) - G(\Gamma_0; \phi)}{t} = \int_{\Gamma_0} -f_n V \cdot n ds,
\end{equation}
where $f_n$ is the dielectric boundary force given in (\ref{eqn:surface_force}) \cite{MikuckiM2014a}.

Numerical techniques for the PBE are mostly focused on the treatment of the discontinuous dielectric function $\epsilon$ and the singular charge density
$q_i \delta(x_i)$. There are three major regularization schemes with which a direct approximation of the singular delta functions can be avoided:
\begin{enumerate}[(I)]
\item  Subtract the potential
\begin{equation} \label{eqn:phi_singular_1}
\phi^s(x) = \sum_{i=1}^N \frac{q_i}{\epsilon_s |x - x_i|}
\end{equation}
induced by the collection of singular charges from the potential in the entire domain $\Omega$, and thus one only needs to
numerically solve the remaining regular potential $\phi^r = \phi - \phi^s$ \cite{ChenL2007a}. In solvent domain $\Omega_s$, this potential is much
larger than the true potential $\phi = \phi^s + \phi^r$, so is the regular potential $\phi^r$ there. Consequently a small
relative error in the numerical solution of $\phi^r$ will present a larger relative error in $\phi$, constituting an unstable
numerical algorithm \cite{HolstM2012a}.
\item Subtract the singular potential (\ref{eqn:phi_singular_1}) only in the solute domain $\Omega_p$. This will create additional
jumps in the remaining regular components. Therefore a harmonic potential $\phi^h$ defined by
\begin{equation} \label{eqn:phi_harmonic}
\epsilon_p \Delta u^h = 0, ~ x \in \Omega_p; \qquad u^h(x) =-\phi^s, ~ x \in \Gamma,
\end{equation}
is added back to partially compensate the jump and to maintain $C^0$ continuity of the remaining regular potential at the interface $\Gamma$ \cite{ChernY2003a}.
This strategy greatly improves the stability of the numerical method. On the other hand, one needs to compute $\nabla_n \phi^h$ on $\Gamma$
to supply the interface conditions for $\phi^r$:
\begin{equation} \label{eqn:PBE_bc_2}
\left[ \phi^r(x) \right ] = 0, \quad \left[ \epsilon \nabla_n \phi^r \right] = \epsilon_p \nabla_n \phi^h.
\end{equation}
Numerical evaluation of $\partial \phi^r/\partial n$ is nontrivial, and might introduce considerable error to the interface condition for $\phi^r$.
A numerical Dirichlet-Neumann mapping based on boundary integral formulations will make this calculation accurate and efficient \cite{ZhouY2007b}.
\item One could choose not to add back the harmonic potential $\phi^h$ so there are only two components in the solute domain.
The absence of this harmonic component leads to the jumps in the potential and its normal derivative on the dielectric interface
\begin{equation} \label{eqn:PBE_bc_3}
\left[ \phi^r(x) \right ] =\phi^s(x), \quad \left[ \epsilon \nabla_n \phi^r \right] = \epsilon_p \nabla_n \phi^s.
\end{equation}
On the other hand, this regularization scheme helps reduce the complexity and improve the accuracy as one does not need to supply a numerical
computed $\nabla_n \phi^h$ to the related interface condition for the regularized Poisson-Boltzmann equation \cite{Geng2017a}.
\end{enumerate}
Our enriched gradient recovery method is compatible with all these regularization schemes. We shall present the gradient recovery
and force calculation along with the following general form of interface conditions:
\begin{equation} \label{eqn:PBE_bc_4}
\left[ \phi^r(x) \right ] = g(x), \quad \left[ \epsilon \nabla_n \phi^r \right] = h(x),
\end{equation}
with known jumps $g(x),h(x)$.

\subsection{Dielectric boundary force and its calculation} \label{subsect:bc}
To get the dielectric boundary force $f_n$, one shall first compute $\phi_s, \phi_p$ and their respective gradients at selected positions (the centroids of
surface triangles, for example) on the interface $\Gamma$. There are there different procedures to complete this calculation:
\begin{enumerate}[(I)]
\item Independently compute $\phi_s,\phi_p$ and their gradients using proper gradient recovery techniques.
\item Compute only one set of $(\phi_p, \nabla \phi_p)$ or $(\phi_s, \nabla \phi_s)$, and use the interface conditions (\ref{eqn:PBE_bc_4}) to get the other.
\item Compute both sets of $(\phi_p, \nabla \phi_p)$ and $(\phi_s, \nabla \phi_s)$ by enforcing the interface conditions (\ref{eqn:PBE_bc_4}) in the gradient
recovery.
\end{enumerate}
Procedure I is inferior to the others  because the resulting dielectric boundary force is not consistent with the interface conditions.
Procedures II and III will be developed in Section \ref{sect:recover} and implemented in Section \ref{sect:numeric}

\section{Gradient Recovery Techniques Enriched with Green's Functions} \label{sect:recover}
Our gradient recovery techniques for PBE are based on the least square reconstruction of the numerical solution, which is also the
basis of polynomial preserving recovery (PPR) techniques \cite{NagaA2005a,ZhangZ2005a} and their recent variations for elliptic
interface problems \cite{GuoH2017a,GuoH2018a,GuoH2018b}. In this section we will first summarize the standard PPR method and
then present its enrichment with Green's functions. This enriched gradient recovery technique will be integrated
with different regularization schemes described in Section \ref{subsect:PBE} to enforce the respective interface conditions.

\subsection{Polynomial preserving recovery and enrichment with Green's functions}
We consider the molecular surface discretized by triangles and assume that the electrostatic potential $\phi$ has been solved from the
Poisson-Boltzmann equation using a finite difference or finite element interface method. The mesh for the interface method may not be
fitting the interface triangulation. The electrostatic force is usually computed at the centroids of triangles which are not
necessarily the mesh nodes. Our method can be used to calculate the solution derivatives at vertices of surface triangles but one needs
a delicate definition of the normal directions at the corner vertices to finally compute the surface force. For a centroid $o$ with
coordinate $\bfx_o=(x_o,y_o.z_o)$ we define $\calN(o,n)$ be the collection of $n$ mesh nodes in solvent domain that are closest to $o$:
\begin{equation} \label{eqn:node_choice_1}
\calN_s(o,n) = \left \{v: v \in V_s, ~ |\bfx_v - \bfx_o| \le |\bfx_z-\bfx_o| \quad \forall z \notin \calN_s(o,n) \right \},
\end{equation}
where $V_s$ are the mesh nodes in $\Omega_s$ or on the molecular surface $\Gamma$. The numerical solution at this collection of nodes will
be fitted by a polynomial $p_s \in \mathbb{P}_{k}$ in the least square sense sampled at these nodes:
\begin{equation} \label{eqn:LS_1}
p_s(x) = \argmin_{p \in \bbP_{k}} \sum_{v \in \calN_s(o,n)} (\phi_s - p)^2(v).
\end{equation}
We then define the recovered gradient at surface position $o$ as
\begin{equation} \label{eqn:LS_2}
G_o(\phi_s) = \nabla p_s(o).
\end{equation}
Alternatively, one can choose the collection of nodes in the solute domain $\Omega_p$ to fit a polynomial $p_p(x)$
and to compute the recovered gradient
\begin{equation} \label{eqn:LS_3}
G_o(\phi_p) = \nabla p_p(o).
\end{equation}
As in most analysis and applications of PPR, we choose $\bbP_{k}$ to be the space of quadratical polynomials, i.e., $k=2$.
\begin{remark}
The choice of sampling nodes (\ref{eqn:node_choice_1}) is for general interface methods. For finite element interface methods one
can take advantage of mesh structure to choose the mesh nodes in the first $l$ layers around the point $o$.
\end{remark}
\begin{figure}[!ht]
\begin{center}
\includegraphics[width=7cm]{ppr.png}
\caption{2-D illustration of polynomial preserving recovery (PPR) and enriched gradient recovery at a surface point $\mathbf{x}_o$
for the interface solution of the Poisson-Boltzmann equation on a Cartesian grid. PPR locally reconstructs solution using a polynomial
in the least square sense at selected nodes $\mathbf{x}_i$ (black dots) in the solvent region $\Omega_s$ (meshed). Enriched recovery adds Green's
functions induced by the charges $q_j$ at the centers $\mathbf{y}_j$ of selected atoms to the polynomial basis. The solvent accessible surface (SAS; red)
is chosen to the dielectric interface. SAS is the trace of the center of probe sphere (blue) modeling water molecule as it rolls over the atoms (surface in cyan)
of the solute molecule.}
\label{fig:recovery}
\end{center}
\end{figure}

As we will demonstrated in Section \ref{sect:numeric}, the polynomial preserving recovery technique is not able to deliver an accurate surface
gradient on grid sizes affordable to practical biomolecular simulations, mostly because polynomial approximation of singular potential components
have large truncation errors. Since these singular components are induced by charges at known positions, we will introduce the Green's function
as basis functions in addition to the quadratical polynomials in the least square reconstruction, c.f. Fig. \ref{fig:recovery}. \textcolor{red}{We anticipate
that the singular component of the solution will be approximated by these enriching basis functions, while the regular component of the solution will
be approximated in the original polynomial space.} For the surface point $o$ given above, we shall choose $m$ charged atoms of the solute molecule that
are closest to $o$:
\begin{equation} \label{eqn:atom_choose_1}
\calA_1(o,m) = \left \{a: a \in \calA_p, ~ |\bfy_a - \bfx_o| \le |\bfy_b-\bfx_o| \quad \forall b \notin \calA_1(o,m) \right \},
\end{equation}
where $\calA_p$ is the set of all charged atoms of the solute molecule. Alternatively, one may choose all charged atoms within a pre-determined
distance $r_c$ to $o$:
\begin{equation} \label{eqn:atom_choose_2}
\calA_2(o,m) = \left \{a: a \in \calA_p, ~ |\bfy_a - \bfx_o| \le r_c \right \}.
\end{equation}
\textcolor{red}{Since the singular components decay like $1/r$, contributions to the potential and force from charged atoms far away from $o$ are secondary,
and can be well approximated in the original polynomial space.}
By enriching the polynomial basis of degree $2$ with Green's functions centered at a total of $m$ selected charged atoms we define a new space of functions
\begin{equation}
\bbQ_{2}  = \mathrm{span} \left \{ 1, x, y, z, x^2, y^2, z^2, xy, yz, xz, \frac{1}{|r_1|}, \cdots, \frac{1}{|r_{m}|} \right \}
\end{equation}
where $x,y,z$ are the Cartesian coordinates with respect to $o$ and $r_m$ is the distance to $m^{th}$ selected charge. We then
fit the numerical solution at mesh nodes in $\calN_s(o,n)$ with a function $q_s \in \mathbb{Q}_{2}$ such that
\begin{equation} \label{eqn:LS_4}
q_s(x) = \argmin_{q \in \bbQ_{2}} \sum_{v \in \calN_s(o,n)} (\phi_s - q)^2(v).
\end{equation}
This amounts to a least square solution $A \bfv = \bfb$ where each row of matrix $A \in \bbR^{n \times (10 + m)}$
is corresponding to a sampling mesh node in $\Omega_s$ with coordinate $(x,y,z)$ and is given by
\begin{align}

A_{i,1:10} = & \left ( 1, \Delta x_i, \Delta y_i, \Delta z_i, \Delta x_i^2, \Delta y_i^2, \Delta z_i^2, \Delta x_i \Delta y_i, \Delta y_i \Delta z_i,
\Delta x_i \Delta z_i \right ),  \label{eqn:Aij_1}   \\
A_{i,11:10+m} = & \left ( \frac{1}{|\bfx_i - \bfy_1|}, \cdots, \frac{1}{|\bfx_i - \bfy_{m}|} \right ), \label{eqn:Aij_2}
\end{align}
where $\Delta x_i = x_i - x_o, \Delta y_i = y_i - y_o, \Delta z_i = z_i - z_o$ and $\bfy_m$ is the position of the $m^{th}$ selected charge.
Since the functions in (\ref{eqn:Aij_2}) are linearly independent
with each other and with the polynomial basis in (\ref{eqn:Aij_1}), the least square problem (\ref{eqn:LS_4}) is solvable.
\begin{remark} \label{remk:why_Greenfunc}
The potential in $\Omega_s$ does not have a decomposition with an analytical reaction
field $\bfF(x,x_i)$ similar to (\ref{eqn:model_GB}). However,
the approximate solution (\ref{eqn:PB_bc}) suggests that it is possible to define a component of the form $1/r$. The space $\bbQ_2$ is constructed
to approximate this component using $1/r_i$ and the remaining component using the regular polynomial basis.
\end{remark}
\begin{remark} \label{remk:Q2_derivative}
The derivative of the Green's funtion can also be written in terms of the same Green's function. Thus the derivative of the basis functions
in $\bbQ_2$ is also in $\bbQ_2$, allowing us to approximate both the potential and its derivatives in the same space.
\end{remark}
\begin{remark} \label{remk:Q2_bornion}
The approximation of the electrostatic potential in $\bbQ_2$ for the Born ion \cite{RouxB1990a} is exact.
\end{remark}

\subsection{Coupled gradient recovery with interface condition enforcement}
We need potentials and potential gradients on both sides of the dielectric interface to compute the dielectric boundary force.
These quantities independently recovered on either side using the method described above do not satisfy the
interface conditions, although these conditions have been enforced in the interface methods for the numerical solution
of the PBE. Procedure II (c.f. Section (\ref{subsect:bc})) recovers the solution on one side and generate the surface potential and
surface gradient on the other side using the interface conditions. This requires only one least square solution for the gradient recovery
at an interface point $o$, and thus is simpler to implement. It is worth noting that if one chooses to
compute $(\phi_p, \nabla \phi_p)$ instead of $(\phi_s, \nabla \phi_s)$, then the enrichment using Green's function is not necessary,
because the singular potential component has been removed from $\Omega_p$ in all three regularization schemes here. Coupled gradient recovery
on both sides of the interface is also possible by enforcing the interface conditions in least square reconstruction.
To this end we shall define in addition to $\calN_s(o,n)$ another collection of $n$ mesh nodes in the solute domain, $\calN_p(o,n)$, that are closest to $o$:
\begin{equation}
\calN_p(o,n) = \left \{v: v \in V_p, ~ |\bfx_v - \bfx_o| \le |\bfx_z-\bfx_o| \quad \forall z \notin \calN_p(o,n) \right \}.
\end{equation}
Two functions, one polynomial $p_p \in \bbP_2$ and the other $q_s \in \bbQ_2$, respectively sampled on $\calN_p(o,n)$ and $\calN_s(o,n)$, will be sought
from the following problem:
\begin{equation} \label{eqn:LS_5}
(p_p(x),q_s(x)) = \argmin_{p \in \bbP_2, q \in \bbQ_{2}} \left [ \sum_{v \in \calN_p(o,n)} (\phi_p - p)^2(v) + \sum_{u \in \calN_s(o,n)} (\phi_s - q)^2(u) \right ],
\end{equation}
subject to the interface conditions:
\begin{equation} \label{eqn:LS_interfacecondition}
q_s(o) - p_p(o) = g(o), \quad \epsilon_s \nabla_n q_s(o) - \epsilon_p \nabla_n p_p(o) = h(o).
\end{equation}
The solution of (\ref{eqn:LS_5}) leads to a least square problem $A \bfv = \bfb$ where the matrix $A \in \bbR^{(2n+2) \times (20+m)}$ whose first $n$ rows correspond
to the approximation of the solution $\phi_p$ respectively at $n$ sampling nodes $(x_i^p, y_i^p, z_i^p)$ in $\Omega_p$, i.e.,
\begin{align}
A_{i,1:10} = & \left ( 1, \Delta x^p_i, \Delta y^p_i, \Delta z^p_i, (\Delta x^p_i)^2, (\Delta y^p_i)^2, (\Delta z^p_i)^2, \Delta x^p_i \Delta y^p_i,
\Delta y^p_i \Delta z^p_i, \Delta x^p_i \Delta z^p_i \right ),  \label{eqn:Aij_3}   \\
A_{i,11:20+m} = & \left ( 0, \cdots, 0 \right ), \label{eqn:Aij_4}
\end{align}
for $1 \le i\le n$. The next $n$ rows of $A$ are for the approximation at sampling nodes $(x_i^s, y_i^s, z_i^s)$ in $\Omega_s$, i.e.,
\begin{align}
A_{n+i,1:10} = & \left ( 0, \cdots, 0 \right ), \label{eqn:Aij_5} \\
A_{n+i,11:20} = & \left ( 1, \Delta x^s_i, \Delta y^s_i, \Delta z^s_i, (\Delta x^s_i)^2, (\Delta y^s_i)^2, (\Delta z^s_i)^2,
\Delta x^s_i \Delta y^s_i, \Delta y^s_i \Delta z^s_i, \Delta x^s_i \Delta z^s_i \right ),  \label{eqn:Aij_6}   \\
A_{n+i,21:20+m} = & \left ( \frac{1}{|\bfx_i - \bfy_1|}, \cdots, \frac{1}{|\bfx_i - \bfy_{m}|} \right ), \label{eqn:Aij_7}
\end{align}
for $1 \le i \le n$. The last two rows of $A$ come from the approximation of the interface conditions (\ref{eqn:LS_interfacecondition}):
\begin{eqnarray}
A_{2n+1,1}  & = &-1, \label{eqn:Aij_7a} \\
A_{2n+1,11} & = & 1, \label{eqn:Aij_7b}  \\
A_{2n-1,(2:10,12:20)} & = & (0, \cdots, 0), \label{eqn:Aij_8}   \\
A_{2n+1,21:20+m} & = & \left ( \frac{1}{|\bfx_o - \bfy_1|}, \cdots, \frac{1}{|\bfx_o - \bfy_{m}|} \right ). \label{eqn:Aij_9} \\
A_{2n+2,(1,5:10,11,15:20)} &  = &(0,\cdots,0), \label{eqn:Aij_10}   \\
A_{2n+2,2:4} & = & -\epsilon_p (n_x, n_y, n_z), \label{eqn:Aij_10a} \\
A_{2n+2,12:14}  & =  & \epsilon_s(n_x, n_y, n_z), \label{eqn:Aij_11}   \\
A_{2n+2,21:20+m} & = & \epsilon_s \left ( \frac{(\bfx_o - \bfy_1) \cdot n}{|\bfx_o - \bfy_1|^3}, \cdots, \frac{(\bfx_o - \bfy_m) \cdot n}{|\bfx_o - \bfy_{m}|^3} \right ).
\label{eqn:Aij_12}
\end{eqnarray}
where $n=(n_x, n_y, n_z)$ is the unit outer normal on the interface. Solution of this least square problem reconstructs functions that
approximate the potential solutions at sampling mesh nodes on two sides of the interface and the jump conditions at the interface point $o$.

\section{Numerical Experiments} \label{sect:numeric}
In this section we validate the accuracy and stability of our enriched gradient recovery method on a set of biomolecules.
Since there do not exist exact solutions of the electrostatic surface potential and gradient for general molecules with
complex surface geometries, we make up potential
fields for accuracy validation of Procedure II by setting the potential in the solvent domain $\Omega_s$ to be
\begin{equation} \label{eqn:exact_field}
\phi_s(\bfx) = \sum_{i=1}^N \frac{q_i}{\epsilon_s| \bfx - \bfy_i|},
\end{equation}
following the approximation boundary condition (\ref{eqn:PB_bc}). The exponential screen term is neglected because the ion strength $\kappa$
is alway small and makes very slight changes to the surface potential and gradient. Procedure II will be used in enriched recovery with analytical fields.
We will compare the classical polynomial preserving (PPR) and enriched gradient recoveries on molecules with the potential field (\ref{eqn:exact_field}).
We will further solve the Poisson-Boltzmann equation with an interface finite difference method MIB \cite{ZhouY2006a,GengW2007a}, and use the solution at a
fine mesh as the reference to check the accuracy of both PPR and our enriched gradient recoveries. Procedure III will be used since we have solutions
available on both sides of the interface. We will use the following absolute and relative errors for quantifying the
accuracy of recovered surface potential $R(\phi)$, gradient $R(\nabla_n \phi)$, and force $R(f_n)$ on the interface:

$$ e(\phi) = \| R(\phi) - \phi \|_{0, \Gamma}, e(\nabla_n \phi) = \| R(\nabla_n \phi) - \nabla_n \phi \|_{0, \Gamma},
e(f_n) = \| R(f_n) - f_n \|_{0, \Gamma}, $$
$$e_r(\phi) = \frac{e(\phi)}{\| \phi \|_{0, \Gamma}}, e_r(\nabla_n \phi) = \frac{e(\nabla_n \phi)}{\| \nabla_n \phi \|_{0, \Gamma}},
e_r(f_n) = \frac{e(f_n)}{\|f_n \|_{0, \Gamma}}.$$

\subsubsection{Born ion} \label{sect:born}
The potential field (\ref{eqn:exact_field}) is exact for Born ion so our enriched gradient recovery method will give result with machine error,
for that $1/| \bfx - \bfx_i|$ is one of the basis function of $\bbQ_2$. For comparison, we present in Table (\ref{tab:1atm_nor})
the errors in surface potential and gradients generated by PPR. It is seen that PPR has proved optimal rate of
convergence and can give fairly good results when the mesh is sufficiently refined. In particular, a mesh space smaller
than $0.1$\ang~is needed to generate a surface potential gradient with a relative error less than $10\%$, an accuracy barely
sufficient for the real biomolecular simulations. The surface potential and magnitude of the surface gradient recovered using
PPR at $h=0.4$\AA~ are shown in Fig. \ref{fig:1atm_exact_nor}.
\begin{table}[!ht]
\begin{center}
	\begin{tabular}{ r llll } \hline
		$h$ & $e(\phi)$  & $e_r(\phi)$ & $e(\nabla \phi)$ & $e_r(\nabla \phi)$ \\ \hline
		0.4  & 1.33 & 9.95e-2 & 7.97 & 1.18 \\
		0.2  & 2.10e-1 & 1.45e-2 & 3.02 & 2.58e-1 \\
		0.1  & 4.15e-2 & 2.83e-3 & 1.05 & 7.69e-2 \\
		0.05 & 6.69e-3 & 4.55e-4 & 3.16e-1 & 2.19e-2 \\ \hline
	\end{tabular}
	\caption{Errors and relative errors in surface potential and gradient recovered using PPR from the potential field (\ref{eqn:exact_field}) in $\Omega_s$
for Born ion.}
\label{tab:1atm_nor}
\end{center}
\end{table}

\begin{figure}[!ht]
\begin{center}
\includegraphics[height=4.5cm]{1atm_exact_nor_pot_h04_d1000.png} \hspace{1cm}
\includegraphics[height=4.5cm]{1atm_exact_nor_grd_h04_d1000.png}
\caption{Surface potential (left; exact value $4.15$ uniform) and magnitude of the surface gradient (right; exact value $4.15$ uniform) recovered using
PPR for Born ion. Mesh size $h=0.4$\AA.}
\label{fig:1atm_exact_nor}
\end{center}
\end{figure}

\subsubsection{Diatomic model} \label{sect:diatomic}
Our diatomic model consists of two unit spheres centered respectively at $(\pm 1,0.0)$, each carrying a unit positive charge at its center.
We first set an analytical potential field (\ref{eqn:exact_field}) in $\Omega_s$ and recover the potential and potential gradient on surface using PPR and enriched
recovery. If both charges are included for enrichment then recovered results will be accurate with machine error. Instead we choose the charge nearest to a surface
point for the recovery there. It is seen from Table (\ref{tab:2atm_exact_compare}) that the enrichment with a single Green's function significantly improves the accuracy
of the gradient recovery. At $h=0.4$\AA.~the relative error in surface gradient from enriched recovery is less than $8\%$ while the traditional PPR method
gives a surface gradient with an relative error more than $80\%$. This significant difference is also highlighted in Fig. \ref{fig:2atm_exact}, where a smoother
and symmetric surface gradient is delivered by the enriched recovery.

\begin{table}[!ht]
\begin{center}
	\begin{tabular}{ r llllcllll } \hline
                \multirow{2}{*}{$h$} & \multicolumn{4}{c}{PPR} &  & \multicolumn{4}{c}{Enriched Recovery} \\  \cline{2-5} \cline{7-10}
		 & $e(\phi)$  & $e_r(\phi)$ & $e(\nabla \phi)$ & $e_r(\nabla \phi)$ & & $e(\phi)$  & $e_r(\phi)$ & $e(\nabla \phi)$ & $e_r(\nabla \phi)$ \\ \hline
		0.4  & 1.35 & 4.80e-2 & 9.60 & 8.17e-1 &  &1.91e-1 & 6.44e-3 & 1.46 & 7.78e-2 \\
		0.2  & 2.77e-1 & 9.51e-3 & 3.98 & 2.31e-1  & &  1.73e-2 & 5.82e-4 & 2.68e-1 & 1.24e-2 \\
		0.1  & 5.57e-2 & 1.90e-3 & 1.39 & 7.08e-2  & & 3.49e-3 & 1.18e-4 & 8.90e-2 & 4.13e-3 \\
		0.05 & 9.12e-3 & 3.11e-4 & 4.23e-1 & 2.06e-2 & & 5.48e-4 & 1.85e-5 & 2.69e-2 & 1.25e-3 \\ \hline
	\end{tabular}
	\caption{Errors and relative errors in surface potential and gradient recovered using PPR and enriched recovery from the
potential field (\ref{eqn:exact_field}) in $\Omega_s$ for the diatomic molecule. The nearest atom is used in the enrichment.}
\label{tab:2atm_exact_compare}
\end{center}
\end{table}

\begin{figure}[!ht]
\begin{center}
\includegraphics[height=3cm]{2atm_exact_r0_d500_h02_pot.png} \hspace{1cm}
\includegraphics[height=3cm]{2atm_exact_r0_d500_h02_grd.png}
\vspace{2mm}
\includegraphics[height=3cm]{2atm_exact_r1_d500_h02_pot.png} \hspace{1cm}
\includegraphics[height=3cm]{2atm_exact_r1_d500_h02_grd.png}
\caption{Surface potential (left) and magnitude of the surface gradient (right) recovered from the potential field (\ref{eqn:exact_field}) in $\Omega_s$ for
diatomic molecule using PPR (top row) and enriched recovery (bottom). The nearest atom is used in the enrichment. Mesh size $h=0.2$\AA.}
\label{fig:2atm_exact}
\end{center}
\end{figure}

Th enriched recovery is also tested and compared to PPR on the numerical solutions of the PBE for this diatomic molecule. Since there does not exist analytical
values of the surface potential and potential gradient, we would use the results of the enriched recovery at a fine mesh as the reference solution for error
measurement. Both atoms are used in the enrichment for optimal recovery. The improvement of the accuracy in the recovered variables due to the enrichment is
evidenced by the results collected in Table \ref{tab:2atm_mib_compare}. This table and the plots in Fig. \ref{fig:2atm_mib_grd_compare} show that both the
PPR and enriched recovery can deliver accurate surface potential (top row) from the numerical solutions of the PBE but the accuracies of recovered gradients differ
significantly. For all mesh sizes we tested the enriched recovery generates surface gradients with a consistent range. In contrast, the results of PPR vary considerably
as the mesh is refined. At a fine mesh with $h=0.05$\AA~, notable variation is still observed in the surface gradient at the center neck of this diatomic molecule.
The coupled enriched recovery is also tested on the numerical solutions of the PBE and the results are rather similar to the recovery based on the solution in $\Omega_s$
only.
\begin{table}[!ht]
\begin{center}
	\begin{tabular}{ r llllcllll } \hline
                \multirow{2}{*}{$h$} & \multicolumn{4}{c}{PPR} &  & \multicolumn{4}{c}{Enriched Recovery} \\  \cline{2-5} \cline{7-10}
		 & $e(\phi)$  & $e_r(\phi)$ & $e(\nabla \phi)$ & $e_r(\nabla \phi)$ & & $e(\phi)$  & $e_r(\phi)$ & $e(\nabla \phi)$ & $e_r(\nabla \phi)$ \\ \hline
		0.4  &  5.46e-1 & 1.86e-2 & 9.51 & 7.72e-1 & & 8.23e-1    & 3.43e-2 & 1.08 & 7.22e-2 \\
		0.2  &  1.86e-1 & 6.32e-3 & 4.37 & 2.53e-1 & & 2.76e-1 & 9.54e-3 & 4.72e-1 & 2.05e-2 \\
		0.1  &  5.71e-2 & 1.94e-3 & 1.54 & 7.74e-2 & & 2.97e-2 & 9.89e-4 & 2.09e-1 & 9.18e-3 \\
		0.05 &  9.76e-3 & 3.31e-4 & 4.58e-1 & 2.20e-2 & &   &  &   & \\ \hline
	\end{tabular}
	\caption{Errors and relative errors in surface potential and gradient recovered using PPR and enriched recovery from the numerical solutions of the PBE
for solvated diatomic molecule. Both atoms are used in the enrichment. Errors are measured using the enriched recovery at $h=0.05$\AA~as the reference solutions.}
\label{tab:2atm_mib_compare}
\end{center}
\end{table}

\begin{figure}[!ht]
\begin{center}
\includegraphics[height=3cm]{2atm_mib_nor_pot_h02_d500.png} \hspace{1cm}
\includegraphics[height=3cm]{2atm_mib_r_pot_h02_d500.png} \\
\includegraphics[height=3cm]{2atm_mib_nor_grd_h02_d500.png} \hspace{1cm}
\includegraphics[height=3cm]{2atm_mib_r_grd_h02_d500.png} \\
\includegraphics[height=3cm]{2atm_mib_nor_grd_h01_d500.png} \hspace{1cm}
\includegraphics[height=3cm]{2atm_mib_r_grd_h01_d500.png} \\
\includegraphics[height=3cm]{2atm_mib_nor_grd_h005_d500.png} \hspace{1cm}
\includegraphics[height=3cm]{2atm_mib_r_grd_h005_d500.png} \\
\caption{Surface potential (top row) and surface gradient (rest) recovered using PPR (left) and enriched recovery (right) from the numerical solutions of the
PBE for solvated diatomic molecule. Both atoms are used in the enrichment. Mesh size $h=0.2,0.2,0.1$ and $0.05$\AA~from top to bottom.}
\label{fig:2atm_mib_grd_compare}
\end{center}
\end{figure}

\subsubsection{Amino Acids} \label{sect:residues}
The numerical experiments above show that the accuracy of recovered gradients can be increased by enriching Green's functions at all charged atoms. For real biomolecular
simulations we might not be able to afford this global enrichment because of the increasing computational cost. Global enrichment might not be necessary because the Green's functions
decay in $1/r$, which indicates the potential is prone to be dominated by the charges nearby. It is therefore possible only to choose a small set of charged atoms
close to the surface point of recovery. Here we will test the enriched recovery on six amino acids: TRY, ASN, ASP, ARG, LYS and GLU dimer. For each amino acid, we set an analytical
potential field using (\ref{eqn:exact_field}) assuming a unit positive charge at each atom.

As seen in Tables (\ref{tab:TRY_exact},\ref{tab:ASN_exact}), the classical PPR ($m=0$) has to use a very fine mesh $h=0.05$\AA to generate a surface gradient with a relative error
about $2\%$. Recovery enriched with Green's functions at $5$ nearest atoms is able to deliver results of the similar accuracy at $h=0.4$\AA. This number is also identified on
the results for other residues, which are neglected as they show the very similar behavior as in these two tables. Enrichment with more Green's functions will further improve
the accuracy of the recovery but for the cause of efficiency we will choose $m=5$ in the numerical experiments below.
\begin{table}[!ht]
\begin{center}
	\begin{tabular}{ r llcllcll } \hline
                \multirow{2}{*}{$m$} & \multicolumn{2}{c}{$h=0.4$} & & \multicolumn{2}{c}{$h=0.2$} & & \multicolumn{2}{c}{$h=0.05$}\\  \cline{2-3} \cline{5-6} \cline{8-9}
    &  $e_r(\phi)$ & $e_r(\nabla \phi)$ & & $e_r(\phi)$ & $e_r(\nabla \phi)$ & & $e_r(\phi)$ & $e_r(\nabla \phi)$ \\ \hline
	0  & 1.16e-1 & 8.00e-1  & &  2.53e-2 & 2.46e-1 & & 7.12e-4 & 2.10e-2 \\
	1  & 3.66e-2 & 1.62e-1  & &  7.06e-3 & 5.68e-2 & & 1.86e-4 & 5.49e-3 \\
	3  & 1.05e-2 & 4.35e-2  & &  1.41e-3 & 1.17e-2 & & 3.87e-5 & 1.17e-3 \\
	5  & 5.20e-3 & 2.22e-2  & &  6.87e-4 & 5.44e-3 & & 2.00e-5 & 5.71e-4 \\
        8  & 2.28e-3 & 9.67e-3  & &  1.81e-4 & 1.47e-3 & & 2.02e-5 & 5.19e-4 \\
       11  & 8.80e-4 & 3.62e-3  & &  5.90e-5 & 4.99e-4 & & 5.70e-6 & 1.39e-4 \\
       15  & 3.13e-4 & 1.35e-3  & &  1.49e-5 & 1.23e-4 & & 1.08e-6 & 2.84e-5 \\
       20  & 1.70e-5 & 9.04e-5  & &  5.80e-7 & 5.17e-6 & & 5.25e-7 & 1.32e-5 \\ \hline
\end{tabular}
	\caption{Relative errors in the potential and gradient obtained from recovery enriched with different number of Green's functions $m$ for TYR surface.
Exact potential field (\ref{eqn:exact_field}) is used. TYR consists of $21$ atoms.}
\label{tab:TRY_exact}
\end{center}
\end{table}

\begin{table}[!ht]
\begin{center}
	\begin{tabular}{ r llcllcll } \hline
                \multirow{2}{*}{$m$} & \multicolumn{2}{c}{$h=0.4$} & & \multicolumn{2}{c}{$h=0.2$} & & \multicolumn{2}{c}{$h=0.05$}\\  \cline{2-3} \cline{5-6} \cline{8-9}
    &  $e_r(\phi)$ & $e_r(\nabla \phi)$ & & $e_r(\phi)$ & $e_r(\nabla \phi)$ & & $e_r(\phi)$ & $e_r(\nabla \phi)$ \\ \hline
	0  & 1.05e-1 & 7.89e-1  & &  2.34e-2 & 2.45e-1 & & 6.56e-4 & 2.09e-2 \\
	1  & 3.94e-2 & 1.79e-1  & &  8.25e-3 & 6.94e-2 & & 2.07e-4 & 6.44e-3 \\
	3  & 1.12e-2 & 4.86e-2  & &  1.88e-3 & 1.50e-2 & & 3.63e-5 & 1.16e-3 \\
	5  & 5.20e-3 & 2.22e-2  & &  6.85e-4 & 5.57e-3 & & 1.41e-5 & 4.28e-4 \\
        7  & 2.56e-3 & 1.11e-2  & &  4.27e-4 & 3.19e-3 & & 2.02e-5 & 5.97e-4 \\
        9  & 1.09e-3 & 4.72e-3  & &  1.34e-4 & 9.83e-4 & & 1.03e-5 & 2.89e-4 \\
       12  & 1.80e-4 & 8.52e-4  & &  1.73e-5 & 1.38e-4 & & 2.26e-6 & 6.23e-5 \\ \hline
\end{tabular}
	\caption{Relative errors in the potential and gradient obtained from recovery enriched with different number of Green's functions $m$ for ASN surface.
Exact potential field (\ref{eqn:exact_field}) is used. ASN consists of $14$ atoms.}
\label{tab:ASN_exact}
\end{center}
\end{table}

We shall also compare PPR with enriched recovery on numerical solutions of the PBE for these amino acids. The surface potential and gradients recovered with the global
enrichment at $h=0.05$\AA~ are taken as the reference solutions to examine the accuracy of recovery scheme enriched with different number of Green's function and
different $h$. The plots for ARG in Fig. \ref{fig:ARG_mib_compare} and for GLU dimer
in Fig. \ref{fig:GLU2_mib_compare} show that surface potential recovered by the classical PPR has an error more than $10\%$ with respect to enriched recovery
on a coarse mesh with $h=0.4$\AA; the surface gradient delivered by PPR is indeed errant on the central neck of the GLU dimer. Despite the salient errors in the surface quantities
recovered by PPR on coarse meshes, these results converge to the same limit as the enriched recovery when the mesh of the numerical solutions for the PBE is sufficiently refined.
In contrast, the enriched recovery produces very consistent surface potentials and surface gradients on all meshes tested, c.f., the
ranges of these surface quantities in these plots.

\begin{figure}[!ht]
\begin{center}
\includegraphics[height=3cm]{arg_mib_r5_d100_h04_pot.png}  \hspace{1cm}
\includegraphics[height=3cm]{arg_mib_r0_d100_h04_pot.png} \\
\includegraphics[height=3cm]{arg_mib_r0_d100_h01_pot.png}  \hspace{1cm}
\includegraphics[height=3cm]{arg_mib_r5_d100_h01_pot.png} \\
\includegraphics[height=3cm]{arg_mib_r0_d100_h04_grd.png}  \hspace{1cm}
\includegraphics[height=3cm]{arg_mib_r5_d100_h04_grd.png} \\
\includegraphics[height=3cm]{arg_mib_r0_d100_h02_grd.png}   \hspace{1cm}
\includegraphics[height=3cm]{arg_mib_r5_d100_h02_grd.png}  \\
\includegraphics[height=3cm]{arg_mib_r0_d100_h01_grd.png}   \hspace{1cm}
\includegraphics[height=3cm]{arg_mib_r5_d100_h01_grd.png}

\caption{Surface potential (top two rows) and magnitude of the surface gradient (others) recovered using PPR (left) and enriched recovery (right) from the numerical solutions of the
PBE for amino acid ARG. Five closest charged atoms are used in the enrichment. Mesh size $h=0.4,0.1,0.4,0.2$ and $0.1$\AA~from top to bottom.}
\label{fig:ARG_mib_compare}
\end{center}
\end{figure}

\begin{figure}[!ht]
\begin{center}
\includegraphics[height=3cm]{glu2_mib_r0_d50_h04_pot.png}  \hspace{1cm}
\includegraphics[height=3cm]{glu2_mib_r5_d50_h04_pot.png} \\
\includegraphics[height=3cm]{glu2_mib_r0_d50_h01_pot.png}  \hspace{1cm}
\includegraphics[height=3cm]{glu2_mib_r5_d50_h01_pot.png} \\
\includegraphics[height=3cm]{glu2_mib_r0_d50_h04_grd.png}  \hspace{1cm}
\includegraphics[height=3cm]{glu2_mib_r5_d50_h04_grd.png} \\
\includegraphics[height=3cm]{glu2_mib_r0_d50_h02_grd.png}   \hspace{1cm}
\includegraphics[height=3cm]{glu2_mib_r5_d50_h02_grd.png}  \\
\includegraphics[height=3cm]{glu2_mib_r0_d50_h01_grd.png}   \hspace{1cm}
\includegraphics[height=3cm]{glu2_mib_r5_d50_h01_grd.png}
\caption{Surface potential (top two rows) and magnitude of the surface gradient (others) recovered using PPR (left) and enriched recovery (right) from the numerical solutions of the
PBE for a GLU dimer. Five closest charged atoms are used in the enrichment. Mesh size $h=0.4,0.1,0.4,0.2$ and $0.1$\AA~from top to bottom.}
\label{fig:GLU2_mib_compare}
\end{center}
\end{figure}

\subsubsection{Proteins} \label{sect:proteins}
Our last numerical experiments are on a low-density lipoprotein receptor 1AJJ which has 37 residues and 519 atoms. The accuracy
of the recovered surface potential and gradient from the numerical solutions of the PBE is documented in Tab. \ref{tab:1ajj_mib_compare},
where the results of the enriched recovery on the uniform mesh with $h=0.05$\AA~is used as the reference. This table provides evidence of the
higher accuracy of the enriched recovery over the classical PPR. For $h=0.4$\AA, a relative error over $40\%$ is found in the surface gradient recovered by PPR.
This error is reduced to about $4\%$ on a very fine mesh with $h=0.1$\AA, which is still larger than the enriched recovery on a coarse mesh with
$h=0.4$\AA. The accuracy improvement due to the enrichment is also indicated in Fig. \ref{fig:1ajj_mib_compare}. Consistent with the observations on model
molecules and amino acids above, enriched recovery produces more stable ranges of surface potential and gradient than PPR.

\begin{table}[!ht]
\begin{center}
	\begin{tabular}{ r llllcllll } \hline
                \multirow{2}{*}{$h$} & \multicolumn{4}{c}{PPR} &  & \multicolumn{4}{c}{Enriched Recovery} \\  \cline{2-5} \cline{7-10}
		 & $e(\phi)$  & $e_r(\phi)$ & $e(\nabla \phi)$ & $e_r(\nabla \phi)$ & & $e(\phi)$  & $e_r(\phi)$ & $e(\nabla \phi)$ & $e_r(\nabla \phi)$ \\ \hline
		0.4  &  1.61 & 1.87e-2 & 11.0 & 4.65e-1 & & 1.52e-1 & 1.76e-3 & 9.49e-1 & 3.03e-2 \\
		0.2  &  3.59e-1 & 4.15e-3 & 4.50 & 1.60e-1 & & 1.73e-2 & 2.00e-4 & 2.18e-1 & 6.96e-3 \\
		0.1  &  6.31e-2 & 7.30e-4 & 1.48 & 4.90e-2 & & 2.70e-3 & 3.13e-5 & 6.28e-2 & 2.01e-3 \\ \hline
	\end{tabular}
	\caption{Errors and relative errors in surface potential and gradient recovered using PPR and enriched recovery from the numerical solutions of the PBE
for solvated protein (PDB ID: 1AJJ). At most five closest charged atoms within a distance of $8$\AA~are used in the enrichment.}
\label{tab:1ajj_mib_compare}
\end{center}
\end{table}

\begin{figure}[!ht]
\begin{center}
\includegraphics[height=3cm]{1ajj_mib_r0_d5_h04_pot.png}  \hspace{1cm}
\includegraphics[height=3cm]{1ajj_mib_r5_d5_h04_pot.png} \\
\includegraphics[height=3cm]{1ajj_mib_r0_d5_h01_pot.png} \hspace{1cm}
\includegraphics[height=3cm]{1ajj_mib_r5_d5_h01_pot.png} \\
\includegraphics[height=3cm]{1ajj_mib_r0_d5_h04_grd.png} \hspace{1cm}
\includegraphics[height=3cm]{1ajj_mib_r5_d5_h04_grd.png} \\
\includegraphics[height=3cm]{1ajj_mib_r0_d5_h02_grd.png} \hspace{1cm}
\includegraphics[height=3cm]{1ajj_mib_r5_d5_h02_grd.png} \\
\includegraphics[height=3cm]{1ajj_mib_r0_d5_h01_grd.png} \hspace{1cm}
\includegraphics[height=3cm]{1ajj_mib_r5_d5_h01_grd.png}
\caption{Surface potential (top two rows) and magnitude of the surface gradient (others) recovered using PPR (left) and enriched recovery (right) from the numerical solutions of the
PBE for solvated protein (PDB ID: 1AJJ). At most five closest charged atoms with a distance of $8$\AA~are used in the enrichment. Mesh size $h=0.4,0.1,0.4,0.2$ and $0.1$\AA~from top to bottom.}
\label{fig:1ajj_mib_compare}
\end{center}
\end{figure}

\section{Concluding Remarks} \label{sect:summary}
We present in this paper a novel numerical method for computing the potential and its gradient on the molecular surfaces from the numerical solutions
of the Poisson-Boltzmann equation, which is widely used in modeling the electrostatic interactions for biomolecules in solvent. Our method reconstructs
a solution in the least square sense locally in the solvent region. In contrast to the classical polynomial preserving recovery (PPR) method that is based
only on the polynomial basis functions, our reconstruction enriches the basis with Green's functions modeling the electrostatic potential field induced
by selected charges. This enrichment is motivated by the Generalized Born method which decomposes the electrostatic potential inside the solute molecules
into a Coulomb potential and a reaction field, and also by the approximate boundary condition for the Poisson-Boltzmann equation. Extensive numerical methods
on Born atom, model diatomic molecule, amino acids, and proteins demonstrate that our enriched recovery is more accurate than PPR on the same mesh, more
stable along the mesh refinement, and indeed provides a stable reference to which the PPR solutions will converge.

Our enriched recovery is easy to implement, and can be readily integrated with other interface finite difference or finite element methods for the
Poisson-Boltzmann equation. As the reconstruction is independently computed at each surface point of interest, parallelization of this enriched recovery is
straightforward. We are currently working to supply the recovered surface gradient to compute the dielectric boundary force and use the force to
drive the large scale continuum molecular deformation as induced due to varying solvation states such as binding of molecules or ligands, changes in
ion concentration, or protonation.

\section*{Acknowledgements}
This work has been partially supported by National Institutes of Health through the grant R01GM117593 as
part of the joint DMS/NIGMS initiative to support research at the interface of the biological and mathematical sciences.


TITLE: conversion from per year to per decade for time-series data
QUESTION [1 upvotes]: I was wondering if this logic makes sense.
I have a time series where the slope is -0.19 over 21 years, so the rate would be -0.19 per year. Now if I want a rate per decade, would it logically make sense to multiple the yearly rate by 10 to achieve a rate per decade?
For this example the rate would become -1.9 per decade. Does this logic make sense or do I need to make additional calculations to perform this conversion?
random generated time series 2000-2020

REPLY [1 votes]: Yes, you can show this explicitly by using "dimensional analysis" with a conversion factor between years and decades.
$$\require{cancel}
-0.19 \frac{\text{units}}{\text{year}}
= -0.19 \frac{\text{units}}{\cancel{\text{year}}} \times 10 \frac{\cancel{\text{years}}}{\text{decade}}
= -1.9 \frac{\text{units}}{\text{decade}}$$


TITLE: Riemann extension theorem with $Re(g)$ bounded
QUESTION [0 upvotes]: Let $D$ be an open subset of the complex plane,  a point $a$ of  $D$ and $g$ a holomorphic function defined on the set $D$ \ ${a}$,if $Re(g)$ is bounded from above,how to show that $g$ can extends to $D$. I want to use Riemann extension Theorem but can't start. I think $a$ is a removable singularity in this case.

REPLY [0 votes]: Begin by introducing $f=e^g$. Details are below.

The function $|f|=e^{\operatorname{Re} g}$ is bounded. By the Riemann extension theorem, $f$ has a holomorphic extension to $D$.  If $f(a)\ne 0$, then we can take a holomorphic branch of $\log f$ in a neighborhood of $a$ and observe that $\log f - g$ is constant (a multiple of $2\pi i$). This provides an extension of $g$.
If $f(a)=0$, write $f(z)=(z-a)^m F(z)$ where $F(a)\ne 0$. Now $\log F$ has a holomorphic branch in a neighborhood of $a$. Hence, $\log F-g$ is holomorphic in a punctured neighborhood of $a$. But locally in this neighborhood, $\log F-\log g = - m\log (z-a)$ up to a multiple of $2\pi i$. This is a contradiction because $m\log (z-a)$ does not have a single-valued holomorphic branch in a  punctured neighborhood of $a$.


TITLE: How can I get unblocked on learning Calculus?
QUESTION [6 upvotes]: I love math and science. In fact paid for $1/2$ my college tuition by tutoring algebra and trigonometry. But when it came to calculus, I became blocked. I understand the concepts of speed and rate of change, but when I start seeing all the symbols like $f'(x)$ and $dx$, my mind goes all fuzzy.
My learning style is that if I can picture something and comprehend it in real life, then I get it in symbols. Algebra, geometry, and trig are all easy enough to picture. But what does it look like in real life to 'take the derivative' or the integral of something? Does it look like a sphere becoming a circle?
Does anyone else's brain work like this? How did you 'get' calculus?

REPLY [1 votes]: I can understand why your mind goes fuzzy when you see things like $x^2 dx$ and $\frac{dy}{dx}$, because these are tricky concepts that are difficult to formalize and frequently used incorrectly, and their geometric meaning is pretty complicated. In fact, even after 5 or so years of university mathematics, I honestly still don't get the geometric meaning of $\frac{dy}{dx}$. However, the meaning of the notation $f'(x)$ should be perfectly clear, and you should be able to visualize it with only a little bit of thought. If you can't do so at the moment, perhaps the problem is that maybe you are missing is the concept of a function. I won't try to explain it here, but I suggest searching online for an explanation of the function concept. Some key words:

Set, function
Domain, codomain; make sure you know what the notation $f : X \rightarrow Y$ means, where $X$ and $Y$ denote sets.
Injective (one-to-one), surjective (onto)
Higher-order function
Lambda abstraction

Another important point is that, given a function $f : \mathbb{R} \rightarrow \mathbb{R},$

$f$ can usually be visualized as a curve in the plane.
the notation $f(x)$ ("$f$ evaluated at $x$") can be visualized as the height of $f$ at $x$.

Make sure you understand this.
Once you've got these concepts, calculus shouldn't make your mind go numb anymore. A simple way to understand the notation $f'(x)$ is that it really means the slope of $f$ at $x$. In other words, its the derivative of $f$, evaluated at $x$. Try to think of derivatives as higher-order functions; the notation $f'(x)$ really means something more like $$(D(f))(x),$$ where $D$ is the derivative function (which is higher-order). The expression $(D(f))(x)$ means: start with the function $f$, then apply $D$ to it, thereby obtaining the corresponding "slope function" $D(f)$, and then evaluate this new function $D(f)$ at $x$; you can visualize this as the height of the slope function $D(f)$ at $x$.


TITLE: Extrema of quartic functions $f(x)=ax^4+bx^3+cx^2+dx$ with $a>0$, $c<0$
QUESTION [3 upvotes]: I am looking for some general statement how the extrema (or only the maximum) of a quartic function looks like.
My quartic polynomial is bounded from below (fourth order coefficient is positive), the second order coefficient is negative (can look like a "mexican hat" potential), and $f(0)=0 $ (no constant term).
To summarize:
$f(x)=ax^4+bx^3+cx^2+dx$ with $a>0, c<0,(e=0)$
What I want to know: Is there any possibility that the function value at the local maximum is negative? If yes, can we say something about relations of the coefficients?
What I have tried so far: If the maximum becomes negative, we are left with two real roots. There is a condition for this to happen, which I can use to derive some relations
(see http://en.wikipedia.org/wiki/Quartic_function#Nature_of_the_roots, which becomes much simpler in my case since $e=0$). The problem is, that these conditions are necessary but not sufficient. It does not exclude the case where the function looks effectively quadratic or only has a saddle point. I still need the maximum. Unfortunately, I wasn't able to go on from here, since the equations become messy for me.
Any hint on literature or general theorems are appreciated.

REPLY [1 votes]: There is no such polynomial. If there is a negative local maximum at a point $x_0$, then
we can write $f(x_0)=-p$, $f'(x_0)=0$ and $f''(x_0)=-2q$ with $p,q \geq 0$. Using
Taylor expansion,
$$
\begin{array}{lcl}
f(x) &=& f(x_0)+f'(x_0)(x-x_0)+\frac{f''(x_0)}{2}(x-x_0)^2+\frac{f^{(3)}(x_0)}{6}(x-x_0)^2+\frac{f^{(4)}(x_0)}{4}(x-x_0)^4 \\
&=& -p-q(x-x_0)^2+w(x-x_0)^3+a(x-x_0)^4 \ (\text{where} \ w=\frac{f^{(3)}(x_0)}{6})\\
&=& ax^4+(w - 4ax_0)x^3 + (6ax_0^2 - 3wx_0 -q)x^2 + (-4ax_0^3 + 3wx_0^2+qx_0)x + (ax_0^4 -wx_0^3 - qx_0^2 - p)
\end{array}
$$
We see then that $f(0)=ax_0^4 -wx_0^3 - qx_0^2 - p$ and $c=6ax_0^2 - 3wx_0 -q$. Now, consider
$$
t_1=3ax_0^4,\ t_2=3p, \ t_3=2qx_0^2, \ t_4=-cx_0^2, \ t_5=3f(0)
$$
Then, by hypothesis, all the $t_i$ are nonnegative. But
$$
t_1+t_2+t_3+t_4+t_5+t_6=
3ax_0^4+3p+2qx_0^2+(-6ax_0^4 + 3wx_0^3 +qx_0^2) +(3ax_0^4 -3wx_0^3 - 3qx_0^2 - 3p)=0.
$$
So all the $t_i$ must in fact be zero. In particular $t_1=0$, whence  $x_0=0$. And also $p=0$ since $t_2=0$. Then the local maximum $f(x_0)$
is zero and is not negative.


text:
\begin{document}

\title{Commutative Bezout domains of stable range 1.5\\
{\normalsize \it Dedicated to the 70-th birthday  of Professor V.V.~Sergeichuk}}

\author[bov]{Victor A. Bovdi\corref{cor}} \ead{vbovdi@gmail.com}
\address[bov]{United Arab Emirates University, Al Ain, UAE}

\author[sh]{Volodymyr P.~ Shchedryk} \ead{shchedrykv@ukr.net}
\address[sh]{Pidstryhach Institute for Applied Problems of Mechanics
and Mathematics,  National Academy of Sciences of Ukraine, Lviv,  Ukraine}

\cortext[cor]{Corresponding author.}

\begin{abstract}
A ring $R$ is said to be of stable range 1.5  if for each  $a, b\in R$ and $0\neq c\in R$  satisfying  $aR + bR + cR = R$ there exists   $r \in R$ such that $(a + br)R + cR = R$.  Let $R$ be a commutative domain in which  all   finitely generated ideals are principal, and let $R$ be of stable range 1.5. Then each matrix $A$ over $R$ is reduced to Smith's canonical form by transformations $PAQ$ in which $P$ and $Q$ are invertible and at least one of them  can be chosen  to be a product of elementary matrices. We generalize Helmer's theorem about the   greatest common divisor  of entries of $A$  over $R$.
\end{abstract}

\begin{keyword}
Commutative Bezout domain\sep Elementary divisor ring\sep Adequate ring\sep Stable range of a ring
\MSC 13C05, 13G05, 16U10
\end{keyword}

\maketitle

\section{Introduction and main  results}

The problem of finding  canonical forms of a matrix  up  to equivalency  is classical. The rings over which matrices are equivalent to certain diagonal matrices have been studied extensively. In the present paper we investigate such question for some classes of  commutative Bezout domains.

The matrix  ${\rm diag}(d_1, d_2,\ldots)$ means  a (possibly rectangular) matrix
having $d_1, d_2,\ldots $ on  main diagonal and zeros elsewhere (by the main
diagonal we mean the one beginning at the upper left corner). We use the following  notations  of commutative rings: $(a_1,\ldots, a_n)$ denotes   the greatest common divisor  of elements $a_1,\ldots, a_n$ and   $a|b$  means that  $a$ is a divisor of  $b$. The set of all matrices of size $n\times m$ over a
 ring $R$ is denoted by $R^{n\times m}$.

An associative (not necessary commutative) ring $R$ is called an  elementary divisor ring (introduced by  I.~Kaplansky in \cite{Kaplansky})  if every (not necessary square)  matrix $A$  over $R$  admits a diagonal reduction, that is, there exist   invertible matrices $P$ and $Q$ over the ring $R$ such that
\begin{equation}\label{EEE:1}
PAQ={\rm diag}(\varphi_1, \; \ldots , \; \varphi_k, \; 0, \; \ldots , \; 0) =\Phi,
\end{equation}
in which  each element $\varphi_i$ is a total divisor of $\varphi_{i+1}$ for $ i=1, \ldots , k-1$ (i.e.
\[
R\varphi_{i+1}R\subseteq \varphi_i R\cap R\varphi_i,
\]
which  is equivalent to $\varphi_{i}| \varphi_{i+1}$ when  $R$ is commutative).
The matrix $\Phi$ is called the {\it Smith normal form}  and
$\varphi_1,\ldots, \varphi_k$ are the   {\it  invariant factors} of the matrix $A$. Examples  of such rings are the ring of integers $\mathbb{Z}$ (see \cite{Smith}), Euclidean rings and  principal ideal rings  (see \cite{Jacobson, Wedderburn}).

A ring $R$ is a  {\it Bezout} ring if  each of its   finitely generated ideals is  principal.
Each matrix from  $R^{1\times n}$ and   $R^{n\times 1}$ over an  elementary divisor ring $R$ admits a diagonal reduction. This is equivalent to the condition that  each   finitely generated ideal in $R$ is principal. Hence an  elementary divisor ring  is a Bezout ring.
Gilman and  Henriksen constructed an example of a commutative Bezout ring which is not an  elementary divisor ring (see \cite[Example 4.11, p.\,382]{Gillman_Henriksen_2}). This raises the problem  whether  an arbitrary  commutative Bezout domain is an  elementary divisor ring. Euclidean rings and principal ideal rings  satisfy the ascending  chain condition on ideals.
However  Helmer  \cite{Helmer} showed that this condition  can be replaced by the less restrictive hypothesis that $R$ is   adequate.

A commutative Bezout domain $R$ is  {\it adequate} if  for  $a,b\in R$ with  $a\not=0$,  there exist $r,s\in R$ such that $a=rs$,  in which   $(r,b)=1$ and if $s'$ is a non-unit divisor of $s$, then  $ (s',b)\neq 1$.  Commutative principal ideal domains and commutative regular rings with identity  are adequate rings; see \cite[Theorem 11, p.\,365]{Gillman_Henriksen} (see also \cite{Chen_Sheibani}).

The proof of the fact that  an adequate ring is an  elementary divisor ring (see \cite[Theorem 3, p.\,234]{Helmer})  was based on  \cite[Theorem 1, p.\,228]{Helmer} which says that  if   $A\in R^{n \times m}$ has  maximal rank over   an adequate ring $R$, then  there is a row   $u =[1, {u_{2}}, {\ldots},   {u_{n}}]\in R^{1\times n}$ such that  the    ${\rm g.c.d.}$ of the entries of $uA$  and the  ${\rm g.c.d.}$ of entries of the matrix $A$  coincide.  It means that there is an invertible matrix
$
U= {\left[\begin{smallmatrix}
  1 & u_2 & \ldots & u_{n} \\ 0 & 1 &
 \ldots & 0 \\ \vdots  & & \ddots & \\ 0 & 0 & \ldots & 1
\end{smallmatrix}\right] }$
such that ${\rm g.c.d.}$ of entries of the matrix $A$ and ${\rm g.c.d.}$ of the elements of the first  row of matrix   $UA$ coincide. This result was generalized for matrices with rank greater then one by Petrychkovych in \cite[Lemma 3.1,p.\,71] {Petrichkovich}.

Deep studies of the theory of   elementary divisor rings increasingly suggest that   methods of  pure  ring theory are insufficient.  Promising  studies were based on the concept of  stable range of rings, introduced by  Bass \cite{Bass_book} as  an  important $K$-theory invariant.

According to \cite[Property $(7.2)_n$, p.\,106]{Bass_Milnor_Serre} (see also   the definition after Lemma 1 in \cite{Vaserstein}),   the {\it stable range} of a ring $R$ is  the smallest positive integer  $n$ such that the following condition holds:

\qquad $(*)_n$\quad  for each  $a_1, \ldots, a_{n + 1} \in R$ satisfying
$a_1R + \cdots + a_{n+1}R = R$
there exist  $b_1, \ldots, b_n \in R$ for which
\[
(a_1 +a_{n+1}b_1)R+ \cdots +(a_n +a_{n+1}b_n)R = R.
\]
If such $n$ does not exist, then the  stable range of $R$ is infinity.

The concept of the  stable range  of a ring turned out to be useful in the study of elementary divisor rings.
In particular, Zabavsky \cite[Theorem 1, p.\,666]{Zabavsky_2} proves that each   elementary divisor ring has stable range $\le 2$. His survey  \cite{Zabavsky_survey} contains results on the problem when  a  commutative Bezout domain is an  elementary divisor ring.

We say that associative  ring $R$ has  \emph{stable range 1.5}  if for each  $a, b\in R$ and $0\neq c\in R$  satisfying  $aR + bR + cR = R$ there exists   $r \in R$ with
\[
(a + br)R + cR = R.
\]
This notion  was  introduced by the second author  \cite{Shchedryk_3} and studied in \cite{Shchedrik_4, Shchedryk}. Commutative principal ideal domains, adequate rings (see   \cite[Prorositions 3.15 and 3.14]{McGovern} and  \cite[Proposition 4]{Anderson_Juett}),
rings of $2\times 2$ matrices over rings listed before  (see \cite[Theorem 5, p.\,856]{Shchedryk_3}) has stable range 1.5.

Evidently each ring  with stable range $1.5$ has Bass stable range $2$. The converse is not always true. For instance, the subring $\mathbb{Z} +x \mathbb{Q}[[x]]$ of the ring of formal power series  $\mathbb{Q}[[x]]$ over the field of rational numbers $\mathbb{Q}$ (see \cite[Example 1, p.\,160]{Henriksen}) has stable range $2$ but not $1.5$ (see \cite[ Example 1.1, p.\,22]{Shchedryk_book}).  This shows that the rings of stable range 1.5 are between the rings of stable range 1 and 2, respectively.

Note that the notion of  stable range 1.5 is closely related to the concept  of almost stable range 1, introduced by McGovern \cite{McGovern}.
A ring $R$ has {\it  almost  stable range} 1 if each proper homomorphic image of $R$ has
stable range 1.
In such rings if $aR + bR + cR = R$,  where $c$ does not belong to the Jacobson radical  ${\mathfrak{J}}(R)$ of $R$,  then  there exists   $r \in R$  (see \cite[Theorem 3.6]{McGovern}) such that
\[
(a + br)R + cR = R.
\]
In commutative rings of almost stable range 1 the condition $c \notin {\mathfrak{J}}(R)$ can be replaced to  $c\neq 0 $  (see \cite[Proposition 4]{Anderson_Juett}). Using this result and \cite[Theorem 3.7]{McGovern} we conclude that each   commutative Bezout domain of stable range  1.5 is an  elementary divisor ring.

Several authors define and study rings of  idempotent stable range \cite{Chen, McGovern}, rings of unit stable range  \cite{Goodearl_Menal},  rings of  neat stable range  \cite{Zabavsky_3}, and rings of square stable range \cite{Khurana_Lam_Wang}.

Our first  result is a  generalization of      Helmer's result  \cite[Theorem 1]{Helmer}.

\begin{theorem}\label{Th:1}
If $R$  is a commutative Bezout domain, then
the following conditions are equivalent:
\begin{itemize}
\item[(i)] $R$ has stable range $1.5;$

\item[(ii)] for each  $A\in R^{n \times m}$  with ${\rm rank} (A)>1$, there exists
 $u =[1,  {u_{2}},   {\ldots},  {u_{n}}]\in R^{1\times n}$
such that
$uA = [{b_{1}}, {b_{2}},  {\ldots},  {b_{m}}]$,
in which  $(b_{1},   b_{2}, \ldots, b_{m} )$ coincides with the  {\rm g.c.d.} of entries of the matrix   $A$.
\end{itemize}
 \end{theorem}

The following example shows that the condition $\rank(A)>1$ in  Theorem $1(ii)$ is essential.
Indeed, let $A:=\left[ \begin{smallmatrix}
5& 0\\
7& 0
\end{smallmatrix}
\right]\in \mathbb{Z}^{2\times 2}$ and  $u=[1,u_2]\in \mathbb{Z}^{1\times 2}$. Then
$uA=[5+7u_2, 0]$. The  {\rm g.c.d.} of entries of matrix $A$ is equal to 1. But
\[
(5+7u_2, 0)=5+7u_2\not=\pm 1
\]
 for any $u_2\in\mathbb{Z}$.
However the ring $\mathbb{Z}$ has stable range $1.5$.

Matrices $P$ and $Q$ satisfying  equality (\ref{EEE:1}) are called {\it transforming matrices} of the matrix $A$.  The set of all transforming matrices $P$ and $Q$ are denoted by $\mathrm{T}_l(A)$ and  $\mathrm{T}_r(A)$, respectively.  An {\it elementary matrix} is a matrix which is obtained from the identity matrix by elementary  transformations.

\begin{theorem}\label{Th:2}
Let $R$  be a commutative Bezout domain of stable range  $1.5$. If
$A=[a_{ij}] \in R^{n\times m}$ with    ${\rm rank}(A)>1$, then  both of the sets  $\mathrm{T}_l(A)$,   $\mathrm{T}_r(A)$  contain an elementary matrices.
\end{theorem}

\section{Proofs}

If  $A, B$ are   matrices such that $A=UBV$ for some  invertible matrices $U$ and $V$, then we say $A\sim B$. We use the following result proved in \cite[Property 6, p.\,50]{Shchedryk_3}.

\begin{lemma}\label{L:1}
Let $R$ be  a commutative Bezout domain of stable range  1.5. Let  $a_1, \ldots, a_n$ be a collection of relatively prime elements in $R$  and $0\ne\psi\in R$. Then there exist $u_1, \ldots, u_n\in R$    such that
\begin{itemize}
\item[\rm (i)] $ u_1a_1 + \cdots + u_na_n = 1$;
\item[\rm (ii)] $ (u_1, \ldots, u_k )= (\psi, u_k)=1$ for each  fixed $2 \leq k\leq n$.
\end{itemize}
\end{lemma}

\begin{proof}[Proof of Theorem \ref{Th:1}]  $(i) \Longrightarrow (ii)$.  Let   $A\in R^{n \times m}$  has   rank greater than 1. Without loss of generality, assume that   $m \geq n$. A ring $R$ is an  elementary divisor ring. Therefore the equation  \eqref{EEE:1} holds for some invertible  $P$ and  $Q$.  Since $\rank(A)> 1$, we have $k\geq 2$ in  \eqref{EEE:1}. Consider
\[
U:=\left[   \begin{smallmatrix}
0 & 0 & 1 & 0 \\
0 & I_{k-2} & 0 & 0 \\
1& 0 & 0 & 0 \\
0 & 0 & 0 & I_{n-k} \\
 \end{smallmatrix}\right]\quad  \text{and}\quad  V:=\left[   \begin{smallmatrix}
0 & 0 & 1 & 0 \\
0 & I_{k-2} & 0 & 0 \\
1& 0 & 0 & 0 \\
0 & 0 & 0 & I_{m-k} \\
 \end{smallmatrix}\right]
\]
in which $I_s$ ($s\geq 1$) is the  identity $s\times s$ matrix and $I_0$ is an empty matrix.
It is easy to check that
\[
 \Phi':= (UP)A(QV)={\rm diag}(\varphi_k, \; \varphi_2,\;
  \varphi_3,\; \ldots , \; \varphi_{k-1}, \;\varphi_1, \;0,\ldots,0).
\]
If  $P_1:=(UP)\det(UP)^{-1}$ and  $Q_1:=(QV)\det(UP)$, then  $P_1AQ_1=\Phi'$ and from     $\det(P_1^{-1})=1$ we obtain that $\det(P_1^{-1})=\sum_{i=1}^n(-1)^{i+1}p_{1i}\Delta_{i}=1$,
where $P_1^{-1}:=[p_{ij}]$ and $\Delta_{1}, \ldots, \Delta_{n} $ are the  corresponding minors.   There exist   $s_{11}, \ldots , s_{1n}$  by Lemma \ref{L:1}, such that
$\sum_{i=1}^n  s_{1i}\Delta_{i} = 1$ and
\begin{equation}\label{EEE:2}
 (s_{11}, s_{12}, \ldots , s_{1k})= (\varphi_k, s_{1k})=1.
\end{equation}
Moreover (see \cite[Property 3, p.\,48]{Shchedryk_3})
\begin{equation}\label{EEE:3}
 [s_{11}, s_{12},  \ldots,   s_{1n}]=[1,  t_2,  \ldots,  t_{n}] P_1^{-1},
\end{equation}
for some   $t_2, \ldots ,  t_{n} \in R$.
The equality $P_1AQ_1=\Phi'$ implies  $AQ_1=P_1^{-1}\Phi'$, and  so by  \eqref{EEE:3} we have
\begin{equation}\label{EEE:4}
\begin{split}
 {\left[\begin{smallmatrix}
  1 & t_2 & \ldots & t_{n} \\
  0 & 1 & \ldots & 0 \\
 \vdots  & & \ddots & \\
 0 & 0 & \ldots & 1
\end{smallmatrix} \right]}
AQ_1&=
 {\left[\begin{smallmatrix}
 s_{11} &  s_{12} & \ldots &  s_{1n} \\
   p_{21} &  p_{22} & \ldots &  p_{2n}\\
  \dots &\dots & \dots &\dots \\
   p_{n1} &  p_{n2} & \ldots &  p_{nn}
\end{smallmatrix} \right]} \Phi'\\

&= \left[\begin{smallmatrix}
 s_{11}\varphi_k &  s_{12}\varphi_2 & \ldots  &   s_{1.k-1}\varphi_{k-1}
 & s_{1k}\varphi_1 & 0 &  {\ldots} & 0  \\
  p_{21}\varphi_k &  p_{22}\varphi_2 & \ldots  &   p_{2.k-1}\varphi_{k-1}
  & p_{2k}\varphi_1 & 0 &  {\ldots} & 0 \\
 \ldots &\ldots & \ldots & \ldots & \ldots & \ldots &\ldots &\ldots \\
  p_{n1}\varphi_k &  p_{n2}\varphi_2 & \ldots  &   p_{n.k-1}\varphi_{k-1}
   & p_{nk}\varphi_1 & 0 &  {\ldots} & 0
\end{smallmatrix}\right].
\end{split}
\end{equation}
According to \eqref{EEE:2},  the  greatest common divisor $\tau$ of elements of the  first row   of the last  matrix in  \eqref{EEE:4} is equal to
\[
\begin{split}
\tau\;=&\; (s_{11}\varphi_k,  s_{12}\varphi_2, \ldots  ,   s_{1.k-1}\varphi_{k-1}, s_{1k}\varphi_1 )\\
=&\varphi_1\left(s_{11}\frac{\varphi_k}{\varphi_1}, s_{12}\frac{\varphi_2}{\varphi_1},\ldots,
\frac{\varphi_{k-1}}{\varphi_1}s_{1.k-1},   s_{1k} \right)\\
=&\varphi_1\left(\left(s_{1k}, s_{11}\frac{\varphi_k}{\varphi_1}\right),
\left(s_{1k}, s_{12}\frac{\varphi_2}{\varphi_1}\right), \ldots  ,   \left(s_{1k}, s_{1.k-1}\frac{\varphi_{k-1}}{\varphi_1}\right) \right)\\
=&\varphi_1(s_{11}, s_{12}, \ldots , s_{1k})=\varphi_1.
\end{split}
\]
Using \eqref{EEE:4} and reasoning as above,   we obtain that
\[
\begin{split}
(1, {t_{2}},  {\ldots},  t_{n}) AQ_1 &=(s_{11},  s_{12},   {\ldots},  s_{1n})\Phi'\\
&=(s_{11}\varphi_k,   s_{12}\varphi_2, \ldots,    s_{1.k-1}\varphi_{k-1},  s_{1k}\varphi_1,   0,\ldots,0)\\
&\sim (\varphi_1 , 0,   {\ldots},  0 ),
\end{split}
\]
so $(1,  {t_{2}},  {\ldots}, t_{n})A
\sim (\varphi_1, 0, {\ldots},  0)$,  in which $\varphi_1$ is the  ${\rm g.c.d.}$  of entries of  $A$.

\medskip

$(i) \Longleftarrow (ii)$.  Let $A:=\left[ \begin{smallmatrix}
a & c \\
b & 0
\end{smallmatrix}\right]\in R^{2\times 2}$,  where    $a,b\in R$,    $0\neq c\in R$, and $ (a,b,c)=1$.  There is $[1, r]\in R^{1\times 2}$ such that $[1, r]A=[b_1,  b_2]$, where
\[
 (b_1, b_2)= (a,b,c)=1.
\]
Hence $ (a+br,c)=1$ and so  $R$ has  stable range  1.5. \end{proof}

\begin{proof}[Proof of Theorem \ref{Th:2}]   By Theorem \ref{Th:1},   there exists $u =[1, {u_{2}},  {\ldots},  {u_{n}}]\in R^{1\times n}$
such that $uA =[{b_{1}},   {\ldots},  {b_{n}}]$,  in which  $\varphi_1:= (b_{1}, \ldots, b_{m} )$ is equal to the   ${\rm g.c.d.}$  of entries of $A$. Thus  $\varphi_1 $ is the first invariant factor of  $A$ (see \eqref{EEE:1}).  Clearly
\[
U_1A=
\textstyle
\begin{bmatrix}
b_{1} &  b_{2} & \ldots &  b_{m} \\
   a_{21} &  a_{22} & \ldots &  a_{2m}\\
  \dots &\dots & \dots &\dots \\
   a_{n1} &  a_{n2} & \ldots &  a_{nm}
\end{bmatrix}, \quad \text{where } U_1:=\textstyle
\begin{bmatrix}
  1 & u_2 & \ldots & u_{n} \\ 0 & 1 &
 \ldots & 0 \\ \vdots  & & \ddots & \\ 0 & 0 & \ldots & 1
\end{bmatrix}.
\]
There exists an invertible  $V_1$  such that
$U_1AV_1=
 \left[\begin{smallmatrix}
\varphi_{1} &  0 & \ldots &  0 \\
   a'_{21} &  a'_{22} & \ldots &  a'_{2m}\\
  \ldots &\ldots & \ldots &\ldots \\
   a'_{n1} &  a'_{n2} & \ldots &  a'_{nm}
\end{smallmatrix}\right]$.

Then   $\varphi_{1}$ is also the   ${\rm g.c.d.}$  of entries of $A_1:=U_1AV_1$.  Thus,
  $\varphi_{1}| a'_{i1}$ for  $i= 2, \ldots , n$ and
\[
U'_1A_1=
\textstyle
\begin{bmatrix}
\varphi_{1} & 0 & \ldots &  0 \\
  0 &  c_{22} & \ldots &  c_{2m}\\
  \dots &\dots & \dots &\dots \\
  0 &  c_{n2} & \ldots &  c_{nm}
\end{bmatrix},  \quad \text{where }
 U'_1:=\textstyle
\begin{bmatrix}
  1 & 0 & \ldots & 0 \\
  -\frac{a'_{21}}{\varphi_{1}} & 1 & \ldots & 0 \\
   \vdots  & & \ddots & \\
   -\frac{a'_{n1}}{\varphi_{1}} & 0 & \ldots & 1
\end{bmatrix}.
\]
Note that $U_1$ and $U'_1$ are elementary matrices.

Consider the  submatrix $B:=\left[\begin{smallmatrix}
    c_{22} & \ldots &  c_{2m}\\
   \dots & \dots &\dots \\
    c_{n2} & \ldots &  c_{nm}
\end{smallmatrix}\right]\in R^{(n-1)\times (m-1)}$ of $U'_1A_1$.
Using the same technique as above, we can find  an  elementary matrix $Z$ and an invertible matrix $W$ such that
$ZB_2W=\left[\begin{smallmatrix}
\varphi_{2} & 0 & \ldots &  0 \\
  0 &  d_{33} & \ldots &  d_{3m}\\
  \dots &\dots & \dots &\dots \\
  0 &  d_{n3} & \ldots &  d_{nm}
\end{smallmatrix}\right]$,
in which  $\varphi_{2}$ is the second invariant factor of $A$.

Evidently $U_2:=1\oplus Z$ is elementary,  $V_2:=1\oplus W\in GL_m(R)$ and
\[
(U_2U'_1U_1)A(V_1V_2)={\rm diag} (\varphi_{1}, \varphi_{2})\oplus F, \qquad (F\in R^{(n-2)\times (m-2)}).
\]
Continuing  this process we obtain that  there exist $P\in GL_n(R)$ and $Q\in GL_m(R)$ such that $PAQ={\rm diag} (\varphi_{1},\ldots,  \varphi_{k})$ in which  $P$ is a product of elementary matrices.

Taking $A^T$ instead of $A$ and applying the same reduction, we construct ''new'' matrices $P$ and $Q$ such that $\Phi:=PA^TQ$ is the Smith canonical matrix and $P$  is a product of  elementary matrices. Since $\Phi$ is symmetric, $\Phi=\Phi^T=Q^TAP^T$, where  $P^T$  is a product of  elementary matrices.
\end{proof}

Note that from Theorem \ref{Th:2}  does not imply  that in (\ref{EEE:1}) both of  $P$ and   $Q$ are elementary.

\section*{ Acknowledgement}
The authors are grateful for the referee's valuable  remarks and suggestions.
The research  was supported  by the UAEU UPAR grant G00002160.

config: arxiv


TITLE: Whether Product of a Negative Definite Matrix and Positive Definite Matrix is Negative Definite?
QUESTION [1 upvotes]: This question arose out of curiosity. Consider a positive definite hermitian matrix $A$ and a negative definite hermitian matrix $B$. Consider the product matrix $C=AB$. I read some where that positive definite matrices are some what like the positive real numbers and similarly negative definite matrices are like negative numbers. So does that kind of similarity carry over here, in detail, is $C$ negative definite?. A quick simulation tells me it is true for randomly generated examples, but I was looking for a mathematical proof if it is true.
EDIT--
I just figured out that this is equivalent to proving product of two positive definite matrices is positive definite.

REPLY [3 votes]: $A=\begin{pmatrix}1 & 1\cr 1 & 4\end{pmatrix}$
$B=\begin{pmatrix}1 & 1\cr 1 & 5\end{pmatrix}$
the above two matrix is positive definte matrix, but $AB$ is not a symmetrical matrix. if we add $AB=BA$,then product of two positive definite matrices is positive definite.


TITLE: Stability of a recurrence relation
QUESTION [1 upvotes]: What I can say about the stability of the recursive formula
$$S_{r+1}(t)=at^{r}-bS_{r}(t), \ \ \ \ r\geq0$$
 for the determination of $S_{r}(t)$?
Assume: $S_{0}(t)=0$ and a,b>0 are constant.
Thanks

REPLY [1 votes]: If you consider $t$ as a fixed parameter, this is just a linear recurrence of first order, with constant coefficients. You can solve it explicitly. Define the generating function $g(z) = \sum_{n \ge 0} S_r(t) z^r$, multiply the recurrence by $z^r$, add up over $r \ge 0$ and recognize some sums:
$\begin{align*}
S_{r + 1}(t)
  &= a t^r - b S_r(t) \\
\sum_{r \ge 0} S_{r + 1}(t) z^r
  &= a \sum_{r \ge 0} (t z)^r - b \sum_{r \ge 0} S_r(t) z^r \\
\frac{g(z) - S_0(t)}{z}
  &= \frac{a}{1 - t z} - b g(z) \\
\end{align*}$
Solve for $g(z)$:
$\begin{align*}
g(z)
  &= \frac{S_0(t) - (t S_0(t) - a) z}{1 - (t - b) z + b t z^2} \\
  &= \frac{S_0(t) t + S_0(t) b - a}{t + b} \cdot \frac{1}{1 + b z}
       + \frac{a}{t + b} \cdot \frac{1}{1 - t z}
\end{align*}$
From here you get your solution as the coefficient of $z^r$, which is easy as it is just two geometric series:
$\begin{align*}
S_r(t)
  &= [z^r] g(z) \\
  &= \frac{S_0(t) t + S_0(t) b - a}{t + b} \cdot (-1)^r b^r
       + \frac{a}{t + b} \cdot t^r
\end{align*}$