# Gradient estimates of solutions to the insulated conductivity problem in
dimension greater than two

YanYan Li111Department of Mathematics, Rutgers University, 110 Frelinghuysen
Rd, Piscataway, NJ 08854, USA. Email<EMAIL_ADDRESS>222Partially
supported by NSF Grants DMS-1501004, DMS-2000261, and Simons Fellows Award
677077. and Zhuolun Yang333Department of Mathematics, Rutgers University, 110
Frelinghuysen Rd, Piscataway, NJ 08854, USA. Email<EMAIL_ADDRESS>444Partially supported by NSF Grants DMS-1501004 and DMS-2000261.

###### Abstract

We study the insulated conductivity problem with inclusions embedded in a
bounded domain in $\mathbb{R}^{n}$. The gradient of solutions may blow up as
$\varepsilon$, the distance between inclusions, approaches to $0$. An upper
bound for the blow up rate was proved to be of order $\varepsilon^{-1/2}$. The
upper bound was known to be sharp in dimension $n=2$. However, whether this
upper bound is sharp in dimension $n\geq 3$ has remained open. In this paper,
we improve the upper bound in dimension $n\geq 3$ to be of order
$\varepsilon^{-1/2+\beta}$, for some $\beta>0$.

##  1 Introduction and main result

Let $\Omega$ be a bounded domain in $\mathbb{R}^{n}$ with $C^{2}$ boundary,
and let $D_{1}^{*}$ and $D_{2}^{*}$ be two open sets whose closure belongs to
$\Omega$, touching only at the origin with the inner normal vector of
$\partial{D}_{1}^{*}$ pointing in the positive $x_{n}$-direction. Translating
$D_{1}^{*}$ and $D_{2}^{*}$ by $\frac{\varepsilon}{2}$ along $x_{n}$-axis, we
obtain

$D_{1}^{\varepsilon}:=D_{1}^{*}+(0^{\prime},\frac{\varepsilon}{2}),\quad\mbox{and}\quad\,D_{2}^{\varepsilon}:=D_{2}^{*}-(0^{\prime},\frac{\varepsilon}{2}).$

When there is no confusion, we drop the superscripts $\varepsilon$ and denote
$D_{1}:=D_{1}^{\varepsilon}$ and $D_{2}:=D_{2}^{\varepsilon}$. Denote
$\widetilde{\Omega}:=\Omega\setminus\overline{(D_{1}\cup D_{2})}$, we consider
the following elliptic equation with Dirichlet boundary data:

$\begin{cases}\mathrm{div}\Big{(}a_{k}(x)\nabla{u}_{k}\Big{)}=0&\mbox{in}~{}\Omega,\\\
u_{k}=\varphi(x)&\mbox{on}~{}\partial\Omega,\end{cases}$ (1.1)

where $\varphi\in{C}^{2}(\partial\Omega)$ is given, and

$a_{k}(x)=\begin{cases}k\in(0,\infty)&\mbox{in}~{}D_{1}\cup{D}_{2},\\\
1&\mbox{in}~{}\widetilde{\Omega}.\end{cases}$

The equation above can be considered as a simple model for electric
conduction, where $a_{k}$ refers to conductivities, which can be assumed to be
1 in the matrix after normalization, and the solution $u_{k}$ gives the
voltage potential. From an engineering point of view, it is very important to
estimate $\nabla u_{k}$, which represents the electric fields, in the narrow
region between the inclusions. This problem is analogous to a linear system of
elasticity studied by Babuška, Andersson, Smith and Levin [5], where they
analyzed numerically that, when the ellitpicity constants are bounded away
from $0$ and infinity, gradient of solutions remain bounded independent of
$\varepsilon$, the distance between inclusions. Bonnetier and Vogelius [12]
proved that for a fixed $k$, $|\nabla u_{k}|$ remains bounded as $\varepsilon$
goes to 0, for circular inclusions $D_{1}$ and $D_{2}$ in dimension $n=2$.
This result was extended by Li and Vogelius [27] to general second order
elliptic equation of divergence form with piecewise Hölder coefficients and
general shape of inclusions $D_{1}$ and $D_{2}$ in any dimension. Furthermore,
they established a stronger $C^{1,\alpha}$ control of $u_{k}$, which is
independent of $\varepsilon$, in each region. Li and Nirenberg [26] further
extended this $C^{1,\alpha}$ result to general second order elliptic systems
of divergence form.

When $k$ equals to $\infty$ (perfect conductor) or $0$ (insulator), it was
shown in [21, 13, 30] that the gradient of solutions generally becomes
unbounded, as $\varepsilon\to 0$. When $k$ goes to $\infty$, $u_{k}$ converges
to the solution of the following perfect conductivity problem:

$\begin{cases}\Delta{u}=0&\mbox{in}~{}\widetilde{\Omega},\\\ u=C_{i}\mbox{
(Constants)}&\mbox{on}~{}\partial{D}_{i},~{}i=1,2,\\\
\int_{\partial{D}_{i}}\frac{\partial{u}}{\partial\nu}=0&i=1,2,\\\
u=\varphi(x)&\mbox{on}~{}\partial\Omega.\end{cases}$ (1.2)

When $k$ goes to $0$, $u_{k}$ converges to the solution of the following
insulated conductivity problem:

$\left\\{\begin{aligned} -\Delta u&=0\quad\mbox{in }\widetilde{\Omega},\\\
\frac{\partial u}{\partial\nu}&=0\quad\mbox{on}~{}\partial{D}_{i},~{}i=1,2,\\\
u&=\varphi\quad\mbox{on }\partial\Omega.\end{aligned}\right.$ (1.3)

See, e.g., Appendix of [6] and [7] for derivations of the above equations.
Here $\nu$ denotes the inward unit normal vectors on $\partial D_{i}$,
$i=1,2$.

Ammari et al. proved in [3] and [4], among other things, the following. Let
$D_{1}^{*}$ and $D_{2}^{*}$ be unit balls in $\mathbb{R}^{2}$, and let $H$ be
a harmonic function in $\mathbb{R}^{2}$. They considered the perfect and
insulated conductivity problems in $\mathbb{R}^{2}$:

$\begin{cases}\Delta{u}=0&\mbox{in}~{}\mathbb{R}^{2}\setminus\overline{(D_{1}\cup
D_{2})},\\\ u=C_{i}\mbox{
(Constants)}&\mbox{on}~{}\partial{D}_{i},~{}i=1,2,\\\
\int_{\partial{D}_{i}}\frac{\partial{u}}{\partial\nu}=0&i=1,2,\\\
u(x)-H(x)=O(|x|^{-1})&\mbox{as}~{}|x|\to\infty,\end{cases}$

and

$\begin{cases}\Delta{u}=0&\mbox{in}~{}\mathbb{R}^{2}\setminus\overline{(D_{1}\cup
D_{2})},\\\ \frac{\partial
u}{\partial\nu}=0&\mbox{on}~{}\partial{D}_{i},~{}i=1,2,\\\
u(x)-H(x)=O(|x|^{-1})&\mbox{as}~{}|x|\to\infty.\end{cases}$

In both cases, they proved that for some $C$ independent of $\varepsilon$,

$\|\nabla u\|_{L^{\infty}(B_{4})}\leq C\varepsilon^{-1/2}.$

They also showed that the upper bounds are optimal in the sense that for
appropriate $H$,

$\|\nabla u\|_{L^{\infty}(B_{4})}\geq\varepsilon^{-1/2}/C.$

Yun extended in [31] and [32] the results allowing $D_{1}^{*}$ and $D_{2}^{*}$
to be any bounded strictly convex smooth domains.

The above gradient estimates were localized and extended to higher dimensions
by Bao, Li and Yin in [6] and [7]. For the perfect conductor case, they
considered problem (1.2) and proved in [6] that

$\begin{cases}\|\nabla u\|_{L^{\infty}(\widetilde{\Omega})}\leq
C\varepsilon^{-1/2}\|\varphi\|_{C^{2}(\partial\Omega)}&\mbox{when}~{}n=2,\\\
\|\nabla u\|_{L^{\infty}(\widetilde{\Omega})}\leq
C|\varepsilon\ln\varepsilon|^{-1}\|\varphi\|_{C^{2}(\partial\Omega)}&\mbox{when}~{}n=3,\\\
\|\nabla u\|_{L^{\infty}(\widetilde{\Omega})}\leq
C\varepsilon^{-1}\|\varphi\|_{C^{2}(\partial\Omega)}&\mbox{when}~{}n\geq
4.\end{cases}$

The above bounds were shown to be optimal in the paper. For further works on
the perfect conductivity problem and closely related ones, see e.g. [1, 10,
11, 15, 19, 20, 23, 24, 25, 8, 9, 16, 22, 14, 2, 18, 29] and the references
therein.

For the insulated problem (1.3), it was proved in [7] that

$\|\nabla u\|_{L^{\infty}(\widetilde{\Omega})}\leq
C\varepsilon^{-1/2}\|\varphi\|_{C^{2}(\partial\Omega)}\quad\mbox{when}~{}n\geq
2.$ (1.4)

The upper bound is optimal for $n=2$ as mentioned earlier, while it was not
known whether it is optimal in dimensions $n\geq 3$.

Yun [33] considered the following insulated problem in $\mathbb{R}^{3}$ minus
2 balls: Let $H$ be a harmonic function in $\mathbb{R}^{3}$,
$D_{1}=B_{1}\left(0,0,1+\frac{\varepsilon}{2}\right)$, and
$D_{2}=B_{1}\left(0,0,-1-\frac{\varepsilon}{2}\right)$,

$\begin{cases}\Delta{u}=0&\mbox{in}~{}\mathbb{R}^{3}\setminus\overline{(D_{1}\cup
D_{2})},\\\ \frac{\partial
u}{\partial\nu}=0&\mbox{on}~{}\partial{D}_{i},~{}i=1,2,\\\
u(x)-H(x)=O(|x|^{-2})&\mbox{as}~{}|x|\to\infty.\end{cases}$

He proved that for some positive constant $C$ independent of $\varepsilon$,

$\max_{|x_{3}|\leq\varepsilon/2}|\nabla u(0,0,x_{3})|\leq
C\varepsilon^{\frac{\sqrt{2}-2}{2}}.$

He also showed that this upper bound of $|\nabla u|$ on the
$\varepsilon$-segment connecting $D_{1}$ and $D_{2}$ is optimal for
$H(x)\equiv x_{1}$.

In this paper, we focus on the insulated conductivity problem (1.3) in
dimension $n\geq 3$, and improve the upper bound (1.4) to the rate
$\varepsilon^{-1/2+\beta}$, for some $\beta>0$. Analogous questions for
elliptic system are still open, and we give some discussions in Section 4. We
point out that the insulator case for Lamé systems in dimension $n=2$ was
studied by Lim and Yu [28].

From now on, we assume that $\partial{D}_{1}^{*}$ and $\partial{D}_{2}^{*}$
are $C^{2}$, and they are relatively convex near the origin. That is, for some
positive constants $R_{0},\kappa$, we assume that when $0<|x^{\prime}|<R_{0}$,
$\partial{D}_{1}^{*}$ and $\partial{D}_{2}^{*}$ are respectively the graphs of
two $C^{2}$ functions $f$ and $g$ in terms of $x^{\prime}$, and

$f(x^{\prime})>g(x^{\prime}),\quad\mbox{for}~{}~{}0<|x^{\prime}|<R_{0},$
$f(0^{\prime})=g(0^{\prime})=0,\quad\nabla_{x^{\prime}}f(0^{\prime})=\nabla_{x^{\prime}}g(0^{\prime})=0,$
(1.5) $\nabla^{2}_{x^{\prime}}(f-g)(x^{\prime})\geq\kappa
I_{n-1},\quad\mbox{for}~{}~{}0<|x^{\prime}|<R_{0},$ (1.6)

where $I_{n-1}$ denotes the $(n-1)\times(n-1)$ identity matrix. Let $a(x)\in
C^{\alpha}(\overline{\widetilde{\Omega}})$, for some $\alpha\in(0,1)$, be a
symmetric, positive definite matrix function satisfying

$\lambda\leq a(x)\leq\Lambda,\quad\mbox{for }x\in\widetilde{\Omega},$

for some positive constants $\lambda,\Lambda$. Let
$\nu=(\nu_{1},\cdots,\nu_{n})$ denote the unit normal vector on $\partial
D_{1}$ and $\partial D_{2}$, pointing towards the interior of $D_{1}$ and
$D_{2}$. We consider the following insulated conductivity problem:

$\left\\{\begin{aligned} -\partial_{i}(a^{ij}\partial_{j}u)&=0\quad\mbox{in
}\widetilde{\Omega},\\\ a^{ij}\partial_{j}u\nu_{i}&=0\quad\mbox{on
}\partial(D_{1}\cup D_{2}),\\\ u&=\varphi\quad\mbox{on
}\partial\Omega,\end{aligned}\right.$ (1.7)

where $\varphi\in C^{2}(\partial\Omega)$ is given.

For $0<\,r\leq\,R_{0}$, we denote

$\displaystyle\Omega_{x_{0},r}:=\left\\{(x^{\prime},x_{n})\in\widetilde{\Omega}~{}\big{|}~{}-\frac{\varepsilon}{2}\right.$
$\displaystyle\left.+g(x^{\prime})<x_{n}<\frac{\varepsilon}{2}+f(x^{\prime}),~{}|x^{\prime}-x_{0}^{\prime}|<r\right\\},$
$\displaystyle\Gamma_{+}:=$
$\displaystyle\left\\{x_{n}=\frac{\varepsilon}{2}+f(x^{\prime}),~{}|x^{\prime}|<R_{0}\right\\},$
(1.8) $\displaystyle\Gamma_{-}:=$
$\displaystyle\left\\{x_{n}=-\frac{\varepsilon}{2}+g(x^{\prime}),~{}|x^{\prime}|<R_{0}\right\\}.$

Since the blow-up of gradient can only occur in the narrow region between
$D_{1}$ and $D_{2}$, we will focus on the following problem near the origin:

$\left\\{\begin{aligned} -\partial_{i}(a^{ij}\partial_{j}u)&=0\quad\mbox{in
}\Omega_{0,R_{0}},\\\ a^{ij}\partial_{j}u\nu_{i}&=0\quad\mbox{on
}\Gamma_{+}\cup\Gamma_{-},\\\ \end{aligned}\right.$ (1.9)

where $\nu=(\nu_{1},\cdots,\nu_{n})$ denotes the unit normal vector on
$\Gamma_{+}$ and $\Gamma_{-}$, pointing upward and downward respectively.

Here is our main result in the paper.

######  Theorem 1.1.

Let $f,g,a,\alpha$ be as above, and let $u\in H^{1}(\Omega_{0,R_{0}})$ be a
solution of (1.9) in dimension $n\geq 3$. There exist positive constants
$r_{0},\beta$ and $C$ depending only on $n$, $\lambda$, $\Lambda$, $R_{0}$,
$\kappa$, $\alpha$, $\|a\|_{C^{\alpha}(\Omega_{0,R_{0}})}$,
$\|f\|_{C^{2}(\\{|x^{\prime}|\leq R_{0}\\})}$ and
$\|g\|_{C^{2}(\\{|x^{\prime}|\leq R_{0}\\})},$ such that

$|\nabla u(x_{0})|\leq
C\|u\|_{L^{\infty}(\Omega_{0,R_{0}})}\left(\varepsilon+|x_{0}^{\prime}|^{2}\right)^{-1/2+\beta},$
(1.10)

for all $x_{0}\in\Omega_{0,r_{0}}$ and $\varepsilon\in(0,1)$.

Let $u\in H^{1}(\widetilde{\Omega})$ be a weak solution of (1.7). By the
maximum principle and the gradient estimates of solutions of elliptic
equations,

$\|u\|_{L^{\infty}(\widetilde{\Omega})}\leq\|\varphi\|_{L^{\infty}(\partial\Omega)},$
(1.11)

and

$\|\nabla u\|_{L^{\infty}(\widetilde{\Omega}\setminus\Omega_{0,r_{0}})}\leq
C\|\varphi\|_{C^{2}(\partial\Omega)}.$

Therefore, a corollary of Theorem 1.1 is as follows.

######  Corollary 1.2.

Let $u\in H^{1}(\widetilde{\Omega})$ be a weak solution of (1.7) in dimension
$n\geq 3$. There exist positive constants $\beta$ and $C$ depending only on
$n$, $\lambda$, $\Lambda$, $R_{0}$, $\kappa$, $\|a\|_{C^{\alpha}}$,
$\|f\|_{C^{2}}$ and $\|g\|_{C^{2}},$ such that

$\|\nabla u\|_{L^{\infty}(\widetilde{\Omega})}\leq
C\|\varphi\|_{C^{2}(\partial\Omega)}\varepsilon^{-\frac{1}{2}+\beta}.$ (1.12)

######  Remark 1.3.

If there are more than two inclusions, estimate (1.12) still holds, with
$\varepsilon$ being the minimal distance between inclusions.

The rest of this paper will be organized as follows. In section 2, we prove a
lemma which is used in the proof of Theorem 1.1. Theorem 1.1 is proved in
Section 3. In section 4, we give a gradient estimate to a problem for elliptic
systems analogous to problem (1.9).

##  2 A regularity lemma

In this section, we give a regularity lemma for elliptic systems (elliptic
equations when $N=1$). Let us first describe the nature of domains and
operators. We define $S$ to be a cylinder

$S=\\{(x^{\prime},x_{n})\in\mathbb{R}^{n}~{}\big{|}~{}|x^{\prime}|<1,|x_{n}|<1\\},$

and some constants $c_{m}$, with $0\leq m\leq l$, such that

$-1=c_{0}<c_{1}<\cdots<c_{l}=1,$

and denote the integer $m_{0}$ to be the integer such that

$c_{m_{0}-1}\leq 0<c_{m_{0}}.$

We divide the domain $S$ into $l$ parts by setting

$\Omega_{m}=\\{x\in S~{}\big{|}~{}c_{m-1}<x_{n}<c_{m}\\},\quad\mbox{for }1\leq
m\leq l.$

For $1\leq\alpha,\beta\leq n,1\leq i,j\leq N$, let $A^{\alpha\beta}_{ij}(x)$
be a function such that

$\|A^{\alpha\beta}_{ij}\|_{L^{\infty}(S)}\leq\Lambda,$
$\int_{S}A^{\alpha\beta}_{ij}(x)\partial_{\alpha}\varphi_{i}(x)\partial_{\beta}\varphi_{j}(x)\geq\lambda\int_{S}|\nabla\varphi|^{2},\quad\forall\varphi\in
H_{0}^{1}(S;\mathbb{R}^{N}),$

for some $\lambda,\Lambda>0$, and for each $1\leq m\leq l$,
$A^{\alpha\beta}_{ij}(x)\in C^{\mu}(\overline{\Omega}_{m})$, for some
$0<\mu<1$. We denote $(A^{\alpha\beta}_{ij}(x))$ by $A(x)$.

For $1\leq\alpha\leq n,1\leq i\leq N$, let

$\displaystyle H(x)$ $\displaystyle=\\{H_{i}\\}\in L^{\infty}(S),$
$\displaystyle G(x)$ $\displaystyle=\\{G_{i}^{\alpha}\\}\in
C^{\mu}(\overline{\Omega}_{m}),$

for all $m=1,\cdots,l$. Then we have the following interior gradient estimate.

######  Lemma 2.1.

Let $A(x)$, $H(x)$ and $G(x)$ be as above. There exists a positive constant
$C$, depending only on $n,\mu,\lambda,\Lambda$ and an upper bound of
$\\{\|A\|_{C^{\mu}(\overline{\Omega}_{m})}\\}_{m=1}^{l}$, such that if $u\in
H^{1}(S;\mathbb{R}^{N})$ is a weak solution to

$\partial_{\alpha}(A^{\alpha\beta}_{ij}(x)\partial_{\beta}u_{j})=H_{i}+\partial_{\alpha}G^{\alpha}_{i}\quad\mbox{in
}S,$

then

$\|u\|_{L^{\infty}(\frac{1}{2}S)}+\|\nabla u\|_{L^{\infty}(\frac{1}{2}S)}\leq
C\left(\|u\|_{L^{2}(S)}+\|H\|_{L^{\infty}(S)}+\max_{1\leq m\leq
l}\|G\|_{C^{\mu}(\overline{\Omega}_{m})}\right).$

######  Remark 2.2.

We point out that the constant $C$ in the Lemma is independent of $l$.

###### Proof.

The proof of Lemma 2.1 is a modification of the proof of Proposition 4.1 in
[26]. Even though the constant $C$ in [26, Proposition 4.1] depends on $l$,
the number of subdomains we divide in the domain $S$, this dependence only
enters in estimating the quantities $\|A-\bar{A}\|_{Y^{1+\alpha,2}}$,
$\|G-\bar{G}\|_{Y^{1+\alpha,2}}$, and $\|H-\bar{H}\|_{Y^{\alpha,2}}$ which
will be defined below. We will show that such quantity is independent of $l$
due to the nature of our domain $S$, and hence the constant $C$ in Lemma 2.1
is independent of $l$.
∎

For $s>0,1<p<\infty$, we define the norm

$\|f\|_{Y^{s,p}}:=\sup_{0<r\leq
1}r^{1-s}\left({\int\diagup}_{rS}|f|^{p}\right)^{1/p}.$

We define a piecewise-constant coefficients $\bar{A}$ associated to $A$ by
setting

$\bar{A}(x):=\left\\{\begin{aligned}
\lim_{x\in\Omega_{m},x\to(0^{\prime},c_{m-1})}&A(x),&&\mbox{if
}x\in\Omega_{m},m>m_{0};\\\ &A(0),&&\mbox{if }x\in\Omega_{m_{0}};\\\
\lim_{x\in\Omega_{m},x\to(0^{\prime},c_{m})}&A(x),&&\mbox{if
}x\in\Omega_{m},m<m_{0}.\end{aligned}\right.$

Similarly, we can define piecewise-constant tensor $\bar{G}$ associated to
$G$. We also define a constant vector $\bar{H}$ associated to $H$ by

$\bar{H}:={\int\diagup}_{S}H.$

######  Lemma 2.3.

Let $A,\bar{A},H,\bar{H},G,\bar{G}$ be as above. Then there exists a positive
constant $C$, depending only on $n$, such that

$\displaystyle\|A-\bar{A}\|_{Y^{1+\mu,2}}$ $\displaystyle\leq C\max_{1\leq
m\leq l}\|A\|_{C^{\alpha}(\overline{S}_{m})},$
$\displaystyle\|G-\bar{G}\|_{Y^{1+\mu,2}}$ $\displaystyle\leq C\max_{1\leq
m\leq l}\|G\|_{C^{\alpha}(\overline{S}_{m})},$
$\displaystyle\|H-\bar{H}\|_{Y^{\mu,2}}$ $\displaystyle\leq
C\|H\|_{L^{\infty}(S)}.$

###### Proof.

The last inequality follows immediately from the definition of $Y^{\mu,2}$ and
$\bar{H}$:

$\displaystyle\|H-\bar{H}\|_{Y^{\mu,2}}\leq\sup_{0<r\leq
1}r^{1-\mu}\left({\int\diagup}_{rS}|H-\bar{H}|^{2}\right)^{1/2}\leq
C\|H\|_{L^{\infty}(S)}.$

By a direct computation, we will have

$\displaystyle\left({\int\diagup}_{rS}|A-\bar{A}|^{2}\right)^{1/2}$
$\displaystyle\leq\left(\frac{1}{|rS|}\sum_{m=1}^{l}\int_{rS\cap
S_{m}}|A(x)-\bar{A}(x)|^{2}\,dx\right)^{1/2}$
$\displaystyle\leq\left[\frac{1}{|rS|}\left(\sum_{m=1}^{m_{0}-1}\|A\|^{2}_{C^{\mu}(\overline{S}_{m})}\int_{rS\cap
S_{m}}|x-(0^{\prime},c_{m})|^{2\mu}\,dx\right.\right.$
$\displaystyle+\|A\|^{2}_{C^{\mu}(\overline{S}_{m_{0}})}\int_{rS\cap
S_{m_{0}}}|x|^{2\mu}\,dx$ $\displaystyle+\sum_{m=m_{0}+1}^{l}$
$\displaystyle\left.\left.\|A\|^{2}_{C^{\mu}(\overline{S}_{m-1})}\int_{rS\cap
S_{m-1}}|x-(0^{\prime},c_{m-1})|^{2\mu}\,dx\right)\right]^{1/2}$
$\displaystyle\leq\max_{1\leq m\leq
l}\|A\|_{C^{\mu}(\overline{S}_{m})}\left({\int\diagup}_{rS}|x|^{2\mu}\,dx\right)^{1/2}$
$\displaystyle\leq C\max_{1\leq m\leq
l}\|A\|_{C^{\mu}(\overline{S}_{m})}r^{\mu}.$

This proves the first inequality. The second inequality follows similarly.

∎

##  3 Proof of Theorem 1.1

In this section, we prove Theorem 1.1. For a small $r_{0}$ independent of
$\varepsilon$, and any $x_{0}\in\Omega_{0,r_{0}}$, we estimate $|\nabla
u(x_{0})|$ as follows: First we establish a Harnack inequality in
$\Omega_{x_{0},r}\setminus\Omega_{x_{0},r/2}$, for $r>0$ in a suitable range.
Together with the maximum principle, this gives the oscillation of $u$ in
$\Omega_{x_{0},\delta}$ a decay $\delta^{2\beta}$, for some positive
$\varepsilon$-independent $\beta$, where

$\delta:=(\varepsilon+|x_{0}^{\prime}|^{2})^{1/2}.$

Then we perform a suitable change of variables in $\Omega_{x_{0},\delta/4}$,
and apply Lemma 2.1 to obtain the desired estimate on $|\nabla u(x_{0})|$.

We fix a $\gamma\in(0,1)$, and let $r_{0}>0$ denote a constant depending only
on $n$, $\kappa$, $\gamma$, $R_{0}$, $\|f\|_{C^{2}}$ and $\|g\|_{C^{2}}$,
whose value will be fixed in the proof. We will always consider
$0<\varepsilon\leq r_{0}^{2}$. First, we require $r_{0}$ small so that for
$|x_{0}^{\prime}|<r_{0}$,

$10\delta<\delta^{1-\gamma}<R_{0}/4.$

######  Lemma 3.1.

There exists a small $r_{0}$, depending only on
$n,\kappa,\gamma,R_{0},\|f\|_{C^{2}}$ and $\|g\|_{C^{2}}$, such that for any
$x_{0}\in\Omega_{0,r_{0}}$, $5|x_{0}^{\prime}|<r<\delta^{1-\gamma}$, if $u\in
H^{1}(\Omega_{x_{0},2r}\setminus\Omega_{x_{0},r/4})$ is a positive solution to
the equation

$\left\\{\begin{aligned}
-\partial_{i}(a^{ij}(x)\partial_{j}u(x))&=0\quad\mbox{in
}\Omega_{x_{0},2r}\setminus\Omega_{x_{0},r/4},\\\
a^{ij}(x)\partial_{j}u(x)\nu_{i}(x)&=0\quad\mbox{on
}(\Gamma_{+}\cup\Gamma_{-})\cap\overline{\Omega_{x_{0},2r}\setminus\Omega_{x_{0},r/4}},\\\
\end{aligned}\right.$

then,

$\sup_{\Omega_{x_{0},r}\setminus\Omega_{x_{0},r/2}}u\leq
C\inf_{\Omega_{x_{0},r}\setminus\Omega_{x_{0},r/2}}u,$ (3.1)

for some constant $C>0$ depending only on
$n,\kappa,\lambda,\Lambda,R_{0},\|f\|_{C^{2}}$ and $\|g\|_{C^{2}},$ but
independent of $r$ and $u$.

###### Proof.

We only need to prove (3.1) for $|x_{0}^{\prime}|>0$, since the
$|x_{0}^{\prime}|=0$ case follows from the result for $|x_{0}^{\prime}|>0$ and
then sending $|x_{0}^{\prime}|$ to $0$. We denote

$h_{r}:=\varepsilon+f\left(x_{0}^{\prime}-\frac{r}{4}\frac{x_{0}^{\prime}}{|x_{0}^{\prime}|}\right)-g\left(x_{0}^{\prime}-\frac{r}{4}\frac{x_{0}^{\prime}}{|x_{0}^{\prime}|}\right),$

and perform a change of variables by setting

$\left\\{\begin{aligned} y^{\prime}&=x^{\prime}-x_{0}^{\prime},\\\
y_{n}&=2h_{r}\left(\frac{x_{n}-g(x^{\prime})+\varepsilon/2}{\varepsilon+f(x^{\prime})-g(x^{\prime})}-\frac{1}{2}\right),\end{aligned}\right.\quad(x^{\prime},x_{n})\in\Omega_{x_{0},2r}\setminus\Omega_{x_{0},r/4}.$
(3.2)

This change of variables maps the domain
$\Omega_{x_{0},2r}\setminus\Omega_{x_{0},r/4}$ to an annular cylinder of
height $h_{r}$, denoted by $Q_{2r,h_{r}}\setminus Q_{r/4,h_{r}}$, where

$Q_{s,t}:=\\{y=(y^{\prime},y_{n})\in\mathbb{R}^{n}~{}\big{|}~{}|y^{\prime}|<s,|y_{n}|<t\\},$
(3.3)

for $s,t>0$. We will show that the Jacobian matrix of the change of variables
(3.2), denoted by $\partial_{x}y$, and its inverse matrix $\partial_{y}x$
satisfy

$|(\partial_{x}y)^{ij}|\leq C,\quad|(\partial_{y}x)^{ij}|\leq C,\quad\mbox{for
}y\in Q_{2r,h_{r}}\setminus Q_{r/4,h_{r}},$ (3.4)

where $C>0$ depends only on $n,\kappa,R_{0},\|f\|_{C^{2}}$ and
$\|g\|_{C^{2}}$.

Let $v(y)=u(x)$, then $v$ satisfies

$\left\\{\begin{aligned}
-\partial_{i}(b^{ij}(y)\partial_{j}v(y))&=0\quad\mbox{in
}Q_{2r,h_{r}}\setminus Q_{r/4,h_{r}},\\\
b^{nj}(y)\partial_{j}v(y)&=0\quad\mbox{on
}\\{y_{n}=-h_{r}\\}\cup\\{y_{n}=h_{r}\\},\end{aligned}\right.$ (3.5)

where the matrix $(b^{ij}(y))$ is given by

$(b^{ij}(y))=\frac{(\partial_{x}y)(a^{ij})(\partial_{x}y)^{t}}{\det(\partial_{x}y)},$
(3.6)

$(\partial_{x}y)^{t}$ is the transpose of $\partial_{x}y$.

It is easy to see that (3.4) implies, using $\lambda\leq(a^{ij})\leq\Lambda$,

$\frac{\lambda}{C}\leq(b^{ij}(y))\leq C\Lambda,\quad\mbox{for }y\in
Q_{2r,h_{r}}\setminus Q_{r/4,h_{r}},$ (3.7)

for some constant $C>0$ depending only on $n,R_{0},\kappa,\|f\|_{C^{2}}$ and
$\|g\|_{C^{2}}$.

In the following and throughout this section, we will denote $A\sim B$, if
there exists a positive universal constant $C$, which might depend on
$n,\lambda,\Lambda,R_{0},\kappa,\|f\|_{C^{2}}$, and $\|g\|_{C^{2}},$ but not
depend on $\varepsilon$, such that $C^{-1}B\leq A\leq CB$.

From (3.2), one can compute that

$\displaystyle(\partial_{x}y)^{ii}$ $\displaystyle=1,\quad\mbox{for }1\leq
i\leq n-1,$ $\displaystyle(\partial_{x}y)^{nn}$
$\displaystyle=\frac{2h_{r}}{\varepsilon+f(x_{0}^{\prime}+y^{\prime})-g(x_{0}^{\prime}+y^{\prime})},$
$\displaystyle(\partial_{x}y)^{ni}$
$\displaystyle=-\frac{2h_{r}\partial_{i}g(x_{0}^{\prime}+y^{\prime})+2y_{n}[\partial_{i}f(x_{0}^{\prime}+y^{\prime})-\partial_{i}g(x_{0}^{\prime}+y^{\prime})]}{\varepsilon+f(x_{0}^{\prime}+y^{\prime})-g(x_{0}^{\prime}+y^{\prime})},\quad\mbox{for
}1\leq i\leq n-1,$ $\displaystyle(\partial_{x}y)^{ij}$
$\displaystyle=0,\quad\mbox{for }1\leq i\leq n-1,j\neq i.$

By (1.5) and (1.6), one can see that

$h_{r}\sim\varepsilon+\left|x_{0}^{\prime}-\frac{r}{4}\frac{x_{0}^{\prime}}{|x_{0}^{\prime}|}\right|^{2}.$

Since $|y_{n}|\leq h_{r}$, by using (1.5) and (1.6), we have that, for $1\leq
i\leq n-1$,

$\displaystyle\left|(\partial_{x}y)^{ni}\right|$ $\displaystyle\leq
C\frac{h_{r}|\partial_{i}g(x_{0}^{\prime}+y^{\prime})|+h_{r}[|\partial_{i}f(x_{0}^{\prime}+y^{\prime})|+|\partial_{i}g(x_{0}^{\prime}+y^{\prime})|]}{\varepsilon+f(x_{0}^{\prime}+y^{\prime})-g(x_{0}^{\prime}+y^{\prime})}$
$\displaystyle\leq
C\frac{h_{r}}{\varepsilon+f(x_{0}^{\prime}+y^{\prime})-g(x_{0}^{\prime}+y^{\prime})}\left[|\partial_{i}f(x_{0}^{\prime}+y^{\prime})|+|\partial_{i}g(x_{0}^{\prime}+y^{\prime})|\right]$
$\displaystyle\leq
C\frac{\varepsilon+\left|x_{0}^{\prime}-\frac{r}{4}\frac{x_{0}^{\prime}}{|x_{0}^{\prime}|}\right|^{2}}{\varepsilon+|x_{0}^{\prime}+y^{\prime}|^{2}}|x_{0}^{\prime}+y^{\prime}|,$

Since $r/4<|y^{\prime}|<2r<2\delta^{1-\gamma}$ and $|x_{0}^{\prime}|<\delta$,
we can estimate

$\displaystyle\left|(\partial_{x}y)^{ni}\right|\leq
C|x_{0}^{\prime}+y^{\prime}|\leq C(|x_{0}^{\prime}|+|y^{\prime}|)\leq
C\delta^{1-\gamma}.$

Next, we will show that

$(\partial_{x}y)^{nn}\sim 1,\quad\mbox{for }y\in Q_{2r,h_{r}}\setminus
Q_{r/4,h_{r}}.$ (3.8)

Indeed, by (1.5) and (1.6), we have

$(\partial_{x}y)^{nn}=\frac{2h_{r}}{\varepsilon+f(x_{0}^{\prime}+y^{\prime})-g(x_{0}^{\prime}+y^{\prime})}\sim\frac{\varepsilon+\left|x_{0}^{\prime}-\frac{r}{4}\frac{x_{0}^{\prime}}{|x_{0}^{\prime}|}\right|^{2}}{\varepsilon+|x_{0}^{\prime}+y^{\prime}|^{2}}.$

Since $|y^{\prime}|>r/4$, it is easy to see

$(\partial_{x}y)^{nn}\leq
C\frac{\varepsilon+\left|x_{0}^{\prime}-\frac{r}{4}\frac{x_{0}^{\prime}}{|x_{0}^{\prime}|}\right|^{2}}{\varepsilon+|x_{0}^{\prime}+y^{\prime}|^{2}}\leq
C.$

On the other hand, since $|y^{\prime}|<2r$ and $|x_{0}^{\prime}|<r/5$, we have

$\displaystyle\varepsilon+\left|x_{0}^{\prime}-\frac{r}{4}\frac{x_{0}^{\prime}}{|x_{0}^{\prime}|}\right|^{2}$
$\displaystyle\geq\varepsilon+\left(\left|\frac{r}{4}\frac{x_{0}^{\prime}}{|x_{0}^{\prime}|}\right|-|x_{0}^{\prime}|\right)^{2}\geq\varepsilon+\left(\frac{r}{4}-\frac{r}{5}\right)^{2}=\varepsilon+\frac{1}{400}r^{2},$

and

$\displaystyle\varepsilon+|x_{0}^{\prime}+y^{\prime}|^{2}$
$\displaystyle\leq\varepsilon+2|x_{0}^{\prime}|^{2}+2|y^{\prime}|^{2}\leq\varepsilon+\frac{2}{25}r^{2}+8r^{2}<\varepsilon+9r^{2}.$

Therefore,

$(\partial_{x}y)^{nn}\geq\frac{1}{C}\frac{\varepsilon+\left|x_{0}^{\prime}-\frac{r}{4}\frac{x_{0}^{\prime}}{|x_{0}^{\prime}|}\right|^{2}}{\varepsilon+|x_{0}^{\prime}+y^{\prime}|^{2}}\geq\frac{1}{C}\frac{\varepsilon+r^{2}/400}{\varepsilon+9r^{2}}\geq\frac{1}{C},$

and (3.8) is verified.

We have shown $(\partial_{x}y)^{ii}\sim 1$, for all $i=1,\cdots,n$, and
$|(\partial_{x}y)^{ij}|\leq C\delta^{1-\gamma}$, for $i\neq j$. We further
require $r_{0}$ to be small enough so that off-diagonal entries of
$\partial_{x}y$ are small. Therefore (3.4) follows. As mentioned earlier,
(3.7) follows from (3.4).

Now we define, for any integer $l$,

$A_{l}:=\left\\{y\in\mathbb{R}^{n}~{}\big{|}~{}\frac{r}{4}<|y^{\prime}|<2r,~{}(l-1)h_{r}<z_{n}<(l+1)h_{r}\right\\}.$

Note that $A_{0}=Q_{2r,h_{r}}\setminus Q_{r/4,h_{r}}.$ For any
$l\in\mathbb{Z}$, we define a new function $\tilde{v}$ by

$\tilde{v}(y):=v\left(y^{\prime},(-1)^{l}\left(y_{n}-2lh_{r}\right)\right),\quad\forall
y\in A_{l}.$

We also define the corresponding coefficients, for $k=1,2,\cdots,n-1$,

$\tilde{b}^{nk}(y)=\tilde{b}^{kn}(y):=(-1)^{l}b^{nk}\left(y^{\prime},(-1)^{l}\left(y_{n}-2lh_{r}\right)\right),\quad\forall
y\in A_{l},$

and for other indices,

$\tilde{b}^{ij}(y):=b^{ij}\left(y^{\prime},(-1)^{l}\left(y_{n}-2lh_{r}\right)\right),\quad\forall
y\in A_{l}.$

Therefore, $\tilde{v}(y)$ and $\tilde{b}^{ij}(y)$ are defined in the infinite
cylinder shell $Q_{2r,\infty}\setminus Q_{r/4,\infty}$. By (3.5),
$\tilde{v}\in H^{1}(Q_{2r,\infty}\setminus Q_{r/4,\infty})$ satisfies

$-\partial_{i}(\tilde{b}^{ij}(y)\partial_{j}\tilde{v}(y))=0\quad\mbox{in
}Q_{2r,\infty}\setminus Q_{r/4,\infty}.$

Note that for any $l\in\mathbb{Z}$ and $y\in A_{l}$,
$\tilde{b}(y)=(\tilde{b}^{ij}(y))$ is orthogonally conjugated to
$b\left(y^{\prime},(-1)^{l}\left(y_{n}-2lh_{r}\right)\right)$. Hence, by
(3.7), we have

$\frac{\lambda}{C}\leq\tilde{b}(y)\leq C\Lambda,\quad\mbox{for }y\in
Q_{2r,\infty}\setminus Q_{r/4,\infty}.$

We restrict the domain to be $Q_{2r,r}\setminus Q_{r/4,r}$, and make the
change of variables $z=y/r$. Set
$\bar{v}(z)=\tilde{v}(y),\bar{b}^{ij}(z)=\tilde{b}^{ij}(y)$, we have

$-\partial_{i}(\bar{b}^{ij}(z)\partial_{j}\bar{v}(z))=0\quad\mbox{in
}Q_{2,1}\setminus Q_{1/4,1},$

and

$\frac{\lambda}{C}\leq\bar{b}(z)\leq C\Lambda,\quad\mbox{for }z\in
Q_{2,1}\setminus Q_{1/4,1}.$

Then by the Harnack inequality for uniformly elliptic equations of divergence
form, see e.g. [17, Theorem 8.20], there exists a constant $C$ depending only
on $n,\kappa,\lambda,\Lambda,R_{0},\|f\|_{C^{2}}$ and $\|g\|_{C^{2}},$ such
that

$\sup_{Q_{1,1/2}\setminus Q_{1/2,1/2}}\bar{v}\leq C\inf_{Q_{1,1/2}\setminus
Q_{1/2,1/2}}\bar{v}.$

In particular, we have

$\sup_{Q_{1,h_{r}/r}\setminus Q_{1/2,h_{r}/r}}\bar{v}\leq
C\inf_{Q_{1,h_{r}/r}\setminus Q_{1/2,h_{r}/r}}\bar{v},$

which is (3.1) after reversing the change of variables.
∎

######  Remark 3.2.

If dimension $n=2$, Lemma 3.1 fails since $Q_{2,1}\setminus
Q_{1/4,1}\subset\mathbb{R}^{2}$ is the union of two disjoint rectangular
domains, and the Harnack inequality cannot be applied on it. In fact, in our
proof of Theorem 1.1, Lemma 3.1 is the only ingredient where dimension $n\geq
3$ is used. As mentioned above, the conclusion of Theorem 1.1 does not hold in
dimension $n=2$.

For any domain $A\subset\widetilde{\Omega}$, we denote the oscillation of $u$
in A by $\mbox{osc}_{A}u:=\sup_{A}u-\inf_{A}u$. Using Lemma 3.1, we obtain a
decay of $\mbox{osc}_{\Omega_{x_{0},\delta}}u$ in $\delta$ as follows.

######  Lemma 3.3.

Let $u$ be a solution of (1.9). For any $x_{0}\in\Omega_{0,r_{0}}$, where
$r_{0}$ is as in Lemma 3.1, there exist positive constants $\sigma$ and $C$,
depending only on $n,R_{0},\kappa,\|f\|_{C^{2}}$ and $\|g\|_{C^{2}},$ such
that

$\mbox{osc}_{\Omega_{x_{0},\delta}}u\leq
C\|u\|_{L^{\infty}(\Omega_{x_{0},\delta^{1-\gamma}})}\delta^{\gamma\sigma}.$
(3.9)

###### Proof.

For simplicity, we drop the $x_{0}$ subscript and denote
$\Omega_{r}=\Omega_{x_{0},r}$ in this proof. Let
$5|x_{0}^{\prime}|<r<\delta^{1-\gamma}$ and
$u_{1}=\sup_{\Omega_{2r}}u-u,u_{2}=u-\inf_{\Omega_{2r}}u.$ By Lemma 3.1, we
have

$\displaystyle\sup_{\Omega_{r}\setminus\Omega_{r/2}}u_{1}$ $\displaystyle\leq
C_{1}\inf_{\Omega_{r}\setminus\Omega_{r/2}}u_{1},$
$\displaystyle\sup_{\Omega_{r}\setminus\Omega_{r/2}}u_{2}$ $\displaystyle\leq
C_{1}\inf_{\Omega_{r}\setminus\Omega_{r/2}}u_{2},$

where $C_{1}>1$ is a constant independent of $r$. Since both $u_{1}$ and
$u_{2}$ satisfy equation (1.9), by the maximum principle,

$\displaystyle\sup_{\Omega_{r}\setminus\Omega_{r/2}}u_{i}=\sup_{\Omega_{r}}u_{i},\quad\inf_{\Omega_{r}\setminus\Omega_{r/2}}u_{i}=\inf_{\Omega_{r}}u_{i},$

for $i=1,2$. Therefore,

$\displaystyle\sup_{\Omega_{r}}u_{1}$ $\displaystyle\leq
C_{1}\inf_{\Omega_{r}}u_{1},$ $\displaystyle\sup_{\Omega_{r}}u_{2}$
$\displaystyle\leq C_{1}\inf_{\Omega_{r}}u_{2}.$

Adding up the above two inequalities, we have

$\mbox{osc}_{\Omega_{r}}u\leq\left(\frac{C_{1}-1}{C_{1}+1}\right)\mbox{osc}_{\Omega_{2r}}u.$

Now we take $\sigma>0$ such that $2^{-\sigma}=\frac{C_{1}-1}{C_{1}+1}$, then

$\mbox{osc}_{\Omega_{r}}u\leq 2^{-\sigma}\mbox{osc}_{\Omega_{2r}}u.$ (3.10)

We start with $r=r_{0}=\delta^{1-\gamma}/2$, and set $r_{i+1}=r_{i}/2$. Keep
iterating (3.10) $k+1$ times, where $k$ satisfies $5\delta\leq
r_{k}<10\delta$, we will have

$\mbox{osc}_{\Omega_{\delta}}u\leq\mbox{osc}_{\Omega_{r_{k}}}u\leq
2^{-(k+1)\sigma}\mbox{osc}_{\Omega_{2r_{0}}}u\leq
2^{1-(k+1)\sigma}\|u\|_{L^{\infty}(\Omega_{\delta^{1-\gamma}})}.$

Since $10\delta>r^{k}=2^{-k}r_{0}=2^{-(k+1)}\delta^{1-\gamma},$ we have
$2^{-(k+1)}<10\delta^{\gamma}$, and hence (3.9) follows immediately.
∎

###### Proof of Theorem 1.1.

Let $u\in H^{1}(\Omega_{0,R_{0}})$ be a solution of (1.9). For
$x_{0}\in\Omega_{0,r_{0}}$, we have, using Lemma 3.3,

$\|u-u_{0}\|_{L^{\infty}(\Omega_{x_{0},\delta})}\leq
C\|u\|_{L^{\infty}(\Omega_{x_{0},\delta^{1-\gamma}})}\delta^{\gamma\sigma},$
(3.11)

for some constant $u_{0}$. We denote $v:=u-u_{0}$, and $v$ satisfies the same
equation (1.9). We work on the domain $\Omega_{x_{0},\delta/4}$, and perform a
change of variables by setting

$\begin{cases}y^{\prime}=\delta^{-1}(x^{\prime}-x_{0}^{\prime}),\\\
y_{n}=\delta^{-1}x_{n}.\end{cases}$ (3.12)

The domain $\Omega_{x_{0},\delta/4}$ becomes

$\displaystyle\left\\{y\in\mathbb{R}^{n}~{}\big{|}~{}|y^{\prime}|\leq\frac{1}{4},\delta^{-1}\left(-\frac{1}{2}\varepsilon+g(x_{0}^{\prime}+\delta
y^{\prime})\right)<y_{n}<\delta^{-1}\left(\frac{1}{2}\varepsilon+f(x_{0}^{\prime}+\delta
y^{\prime})\right)\right\\}.$

We make a change of variables again by

$\begin{cases}z^{\prime}=4y^{\prime},\\\ z_{n}=2\delta\left(\frac{\delta
y_{n}-g(x_{0}^{\prime}+\delta
y^{\prime})+\varepsilon/2}{\varepsilon+f(x_{0}^{\prime}+\delta
y^{\prime})-g(x_{0}^{\prime}+\delta
y^{\prime})}-\frac{1}{2}\right).\end{cases}$ (3.13)

Now the domain in $z$-variables becomes a thin plate $Q_{1,\delta}$, where
$Q_{s,t}$ is defined as in (3.3). Let $w(z)=v(x)$, then $w$ satisfies

$\left\\{\begin{aligned}
-\partial_{i}(b^{ij}(z)\partial_{j}w(z))&=0\quad\mbox{in }Q_{1,\delta},\\\
b^{nj}(z)\partial_{j}w(z)&=0\quad\mbox{on
}\\{z_{n}=-\delta\\}\cup\\{z_{n}=\delta\\},\end{aligned}\right.$ (3.14)

where the matrix $b(z)=(b^{ij}(z))$ is given by

$(b^{ij}(z))=\frac{(\partial_{y}z)(a^{ij})(\partial_{y}z)^{t}}{\det(\partial_{y}z)}.$
(3.15)

Similar to the proof of Lemma 3.1, we will show that the Jacobian matrix of
the change of variables (3.13), denoted by $\partial_{y}z$, and its inverse
matrix $\partial_{z}y$ satisfy

$|(\partial_{y}z)^{ij}|\leq C,\quad|(\partial_{z}y)^{ij}|\leq C,\quad\mbox{for
}z\in Q_{1,\delta},$ (3.16)

where $C>0$ depends only on $n,\kappa,R_{0},\|f\|_{C^{2}}$ and
$\|g\|_{C^{2}}$. This leads to

$\frac{\lambda}{C}\leq b(z)\leq C\Lambda,\quad\mbox{for }z\in Q_{1,\delta}.$
(3.17)

From (3.13), one can compute that

$\displaystyle(\partial_{y}z)^{ii}$ $\displaystyle=4,\quad\mbox{for }1\leq
i\leq n-1,$ $\displaystyle(\partial_{y}z)^{nn}$
$\displaystyle=\frac{2\delta^{2}}{\varepsilon+f(x_{0}^{\prime}+\delta
z^{\prime}/4)-g(x_{0}^{\prime}+\delta z^{\prime}/4)},$
$\displaystyle(\partial_{y}z)^{ni}$
$\displaystyle=-\frac{2\delta^{2}\partial_{i}g(x_{0}^{\prime}+\delta
z^{\prime}/4)+2z_{n}\delta[\partial_{i}f(x_{0}^{\prime}+\delta
z^{\prime}/4)-\partial_{i}g(x_{0}^{\prime}+\delta
z^{\prime}/4)]}{\varepsilon+f(x_{0}^{\prime}+\delta
z^{\prime}/4)-g(x_{0}^{\prime}+\delta z^{\prime}/4)}$
$\displaystyle~{}\quad\mbox{for }1\leq i\leq n-1,$
$\displaystyle(\partial_{y}z)^{ij}$ $\displaystyle=0,\quad\mbox{for }1\leq
i\leq n-1,j\neq i.$

First we will show that

$(\partial_{y}z)^{nn}\sim 1,\quad\mbox{for }z\in Q_{1,\delta}.$ (3.18)

Since $|z^{\prime}|<1$ and $|x_{0}^{\prime}|<\delta$, it is easy to see that

$(\partial_{y}z)^{nn}\sim\frac{\delta^{2}}{\varepsilon+|x_{0}^{\prime}+\delta
z^{\prime}/4|^{2}}\geq\frac{\delta^{2}}{\varepsilon+C\delta^{2}}\geq\frac{1}{C},\quad\mbox{for
}z\in Q_{1,\delta},$

due to (1.5) and (1.6). On the other hand,

$\displaystyle(\partial_{y}z)^{nn}$
$\displaystyle\sim\frac{\delta^{2}}{\varepsilon+|x_{0}^{\prime}+\delta
z^{\prime}/4|^{2}}$
$\displaystyle=\frac{\delta^{2}}{\varepsilon+|x_{0}^{\prime}|^{2}+(1/4)^{2}\delta^{2}|z^{\prime}|^{2}+\delta
x_{0}^{\prime}\cdot z^{\prime}/2}$
$\displaystyle\leq\frac{\delta^{2}}{\delta^{2}+(1/4)^{2}|z^{\prime}|^{2}\delta^{2}-|z^{\prime}||x_{0}^{\prime}|\delta/2}$
$\displaystyle\leq\frac{\delta^{2}}{(1+(1/4)^{2}|z^{\prime}|^{2}-1/2)\delta^{2}}\leq
C,\quad\mbox{for }z\in Q_{1,\delta}.$

Therefore, (3.18) is verified.

Since $|z_{n}|<\delta$, $|z^{\prime}|<1$ and $|x_{0}^{\prime}|<\delta$, by
(1.5) and (1.6), for $1\leq i\leq n-1$,

$\displaystyle|(\partial_{y}z)^{ni}|$
$\displaystyle\leq\frac{2\delta^{2}|\partial_{i}g(x_{0}^{\prime}+\delta
z^{\prime}/4)|+2\delta^{2}[|\partial_{i}f(x_{0}^{\prime}+\delta
z^{\prime}/4)|+|\partial_{i}g(x_{0}^{\prime}+\delta
z^{\prime}/4)|]}{\varepsilon+f(x_{0}^{\prime}+\delta
z^{\prime}/4)-g(x_{0}^{\prime}+\delta z^{\prime}/4)}$
$\displaystyle\leq\frac{C\delta^{2}}{\varepsilon+f(x_{0}^{\prime}+\delta
z^{\prime}/4)-g(x_{0}^{\prime}+\delta
z^{\prime}/4)}[|\partial_{i}f(x_{0}^{\prime}+\delta
z^{\prime}/4)|+|\partial_{i}g(x_{0}^{\prime}+\delta z^{\prime}/4)|]$
$\displaystyle\leq C\frac{\delta^{2}}{\varepsilon+|x_{0}^{\prime}+\delta
z^{\prime}/4|^{2}}|x_{0}^{\prime}+\delta z^{\prime}/4|$ $\displaystyle\leq
C(|x_{0}^{\prime}|+\delta|z^{\prime}|)\leq C\delta,$

where in the last line, we have used the same arguments in showing
$(\partial_{y}z)^{nn}\leq C$ earlier.

We have shown $(\partial_{y}z)^{ii}\sim 1$, for all $i=1,\cdots,n$, and
$|(\partial_{y}z)^{ij}|\leq C\delta$, for $i\neq j$. We further require
$r_{0}$ to be small enough so that off-diagonal entries are small. Therefore
(3.16) follows. As mentioned earlier, (3.17) follows from (3.16).

Next, we will show

$\|b\|_{C^{\alpha}(\overline{Q}_{1,\delta})}\leq C,$ (3.19)

for some $C>0$ depending only on
$n,\kappa,R_{0},\|a\|_{C^{\alpha}},\|f\|_{C^{2}}$ and $\|g\|_{C^{2}}$, by
showing

$|\nabla_{z}(\partial_{y}z)^{ij}(z)|\leq
C,\quad\left|\nabla_{z}\frac{1}{\det(\partial_{y}z)}\right|\leq
C,\quad\mbox{for }z\in Q_{1,\delta}.$ (3.20)

Then (3.19) follows from (3.20), (3.15), and $\|a\|_{C^{\alpha}}\leq C$.

By a straightforward computation, we have, for any $i=1,\cdots,n-1$,

$\displaystyle\left|\partial_{z_{i}}\frac{1}{\det(\partial_{y}z)}\right|$
$\displaystyle=\left|\partial_{z_{i}}\left(\frac{\varepsilon+f(x_{0}^{\prime}+\delta
z^{\prime}/4)-g(x_{0}^{\prime}+\delta z^{\prime}/4)}{2\cdot
4^{n-1}\delta^{2}}\right)\right|$
$\displaystyle=\left|\frac{\delta[\partial_{i}f(x_{0}^{\prime}+\delta
z^{\prime}/4)-\partial_{i}g(x_{0}^{\prime}+\delta z^{\prime}/4)]}{2\cdot
4^{n-1}\delta^{2}}\right|$
$\displaystyle\leq\frac{C}{\delta}[|\partial_{i}f(x_{0}^{\prime}+\delta
z^{\prime}/4)|+|\partial_{i}g(x_{0}^{\prime}+\delta z^{\prime}/4)|]$
$\displaystyle\leq\frac{C}{\delta}|x_{0}^{\prime}+\delta z^{\prime}/4|\leq
C,\quad\mbox{for }z\in Q_{1,\delta},$

where in the last inequality, (1.5) and (1.6) have been used. For any
$i=1,\cdots,n-1$,

$\displaystyle|\partial_{z_{i}}(\partial_{y}z)^{nn}|$
$\displaystyle=\left|\frac{2\delta^{3}[\partial_{i}f(x_{0}^{\prime}+\delta
z^{\prime}/4)-\partial_{i}g(x_{0}^{\prime}+\delta
z^{\prime}/4)]}{(\varepsilon+f(x_{0}^{\prime}+\delta
z^{\prime}/4)-g(x_{0}^{\prime}+\delta z^{\prime}/4))^{2}}\right|$
$\displaystyle\leq\frac{C\delta^{3}}{(\varepsilon+|x_{0}^{\prime}+\delta
z^{\prime}/4|^{2})^{2}}|x_{0}^{\prime}+\delta z^{\prime}/4|$
$\displaystyle\leq\frac{C\delta^{3}}{\delta^{4}}(|x_{0}^{\prime}|+|\delta
z^{\prime}|)\leq C,\quad\mbox{for }z\in Q_{1,\delta},$

where in the last line, we have used the same arguments in showing
$(\partial_{y}z)^{nn}\leq C$ earlier. Similar computations apply to
$\partial_{z_{i}}(\partial_{y}z)^{ni}$, for $i=1,\cdots,n-1$, and we have

$|\partial_{z_{i}}(\partial_{y}z)^{ni}|\leq C,\quad\mbox{for }z\in
Q_{1,\delta}.$

Finally, we compute, for $i=1,\cdots,n-1$,

$\displaystyle|\partial_{z_{n}}(\partial_{y}z)^{ni}|$
$\displaystyle=\left|\frac{2\delta[\partial_{i}f(x_{0}^{\prime}+\delta
z^{\prime}/4)-\partial_{i}g(x_{0}^{\prime}+\delta
z^{\prime}/4)]}{\varepsilon+f(x_{0}^{\prime}+\delta
z^{\prime}/4)-g(x_{0}^{\prime}+\delta z^{\prime}/4)}\right|$
$\displaystyle\leq\frac{C\delta|x_{0}^{\prime}+\delta
z^{\prime}/4|}{\varepsilon+|x_{0}^{\prime}+\delta z^{\prime}/4|^{2}}\leq
C,\quad\mbox{for }z\in Q_{1,\delta}.$

Therefore, (3.20) is verified, and hence (3.19) follows as mentioned above.

Now we define

$S_{l}:=\left\\{z\in\mathbb{R}^{n}~{}\big{|}~{}|z^{\prime}|<1,~{}(l-1)\delta<z_{n}<(l+1)\delta\right\\}$

for any integer $l$, and

$S:=\left\\{z\in\mathbb{R}^{n}~{}\big{|}~{}|z^{\prime}|<1,~{}|z_{n}|<1\right\\}.$

Note that $Q_{1,\delta}=S_{0}$. As in the proof of Lemma 3.1, we define, for
any $l\in\mathbb{Z}$, a new function $\tilde{w}$ by setting

$\tilde{w}(z):=w\left(z^{\prime},(-1)^{l}\left(z_{n}-2l\delta\right)\right),\quad\forall
z\in S_{l}.$

We also define the corresponding coefficients, for $k=1,2,\cdots,n-1$,

$\tilde{b}^{nk}(z)=\tilde{b}^{kn}(z):=(-1)^{l}b^{nk}\left(z^{\prime},(-1)^{l}\left(z_{n}-2l\delta\right)\right),\quad\forall
z\in S_{l},$

and for other indices,

$\tilde{b}^{ij}(z):=b^{ij}\left(z^{\prime},(-1)^{l}\left(z_{n}-2l\delta\right)\right),\quad\forall
y\in S_{l}.$

Then $\tilde{w}$ and $\tilde{b}^{ij}$ are defined in the infinite cylinder
$Q_{1,\infty}$. By (3.14), $\tilde{w}$ satisfies the equation

$-\partial_{i}(\tilde{b}^{ij}\partial_{j}\tilde{w})=0,\quad\mbox{in
}Q_{1,\infty}.$

Note that for any $l\in\mathbb{Z}$, $\tilde{b}(z)$ is orthogonally conjugated
to $b\left(z^{\prime},(-1)^{l}\left(z_{n}-2l\delta\right)\right),$ for $z\in
S_{l}$. Hence, by (3.17), we have

$\frac{\lambda}{C}\leq\tilde{b}(z)\leq C\Lambda,\quad\mbox{for }z\in
Q_{1,\infty},$

and, by (3.19),

$\|\tilde{b}\|_{C^{\alpha}(\overline{S}_{l})}\leq C,\quad\forall
l\in\mathbb{Z}.$

Apply Lemma 2.1 on $S$ with $N=1$, we have

$\|\nabla\tilde{w}\|_{L^{\infty}(\frac{1}{2}S)}\leq
C\|\tilde{w}\|_{L^{2}(S)}.$

It follows that

$\|\nabla
w\|_{L^{\infty}(Q_{1/2,\delta})}\leq\frac{C}{\delta}\|w\|_{L^{2}(Q_{1,\delta})}\leq
C\|w\|_{L^{\infty}(Q_{1,\delta})},$

for some positive constant $C$, depending only on
$n,\alpha,R_{0},\kappa,\lambda,\Lambda,\|a\|_{C^{\alpha}},\|f\|_{C^{2}}$ and
$\|g\|_{C^{2}}$.

Since $\|(\partial_{z}y)\|_{L^{\infty}(Q_{1,\delta})}\leq C$ by (3.16), where
$C$ depends only on $R_{0},\kappa,\|f\|_{C^{2}}$ and $\|g\|_{C^{2}}$, and in
particular, is independent of $\varepsilon$ and $\delta$. Reversing the change
of variables (3.13) and (3.12), we have

$\delta\|\nabla v\|_{L^{\infty}(\Omega_{x_{0},\delta/8})}\leq
C\|v\|_{L^{\infty}(\Omega_{x_{0},\delta/4})}\leq
C\|u\|_{L^{\infty}(\Omega_{x_{0},\delta^{1-\gamma}})}\delta^{\gamma\sigma}$

by (3.11). In particular, this implies

$|\nabla u(x_{0})|\leq
C\|u\|_{L^{\infty}(\Omega_{x_{0},\delta^{1-\gamma}})}\delta^{-1+\gamma\sigma},$

and it concludes the proof of Theorem 1.1 after taking $\beta=\gamma\sigma/2$.
∎

##  4 Gradient estimates of elliptic systems

A natural question is whether the estimate in Theorem 1.1 can be extended to
elliptic systems of divergence form. We tend to believe that the answer to
this question is affirmative, and plan to pursue this in a subsequent paper.
Following closely the proof of (1.4) in [7], we give a preliminary gradient
estimate of elliptic systems in this section.

We consider the vector-valued function $u=(u_{1},\cdots,u_{N})$, and for
$1\leq\alpha,\beta\leq n,1\leq i,j\leq N$, let $A^{\alpha\beta}_{ij}(x)$ be a
function such that

$\|A^{\alpha\beta}_{ij}\|_{L^{\infty}(\Omega_{0,R_{0}})}\leq\Lambda,$
$\int_{\Omega_{0,R_{0}}}A^{\alpha\beta}_{ij}(x)\partial_{\alpha}\varphi_{i}(x)\partial_{\beta}\varphi_{j}(x)\geq\lambda\int_{\Omega_{0,R_{0}}}|\nabla\varphi|^{2},\quad\forall\varphi\in
H_{0}^{1}(\Omega_{0,R_{0}};\mathbb{R}^{N}),$

for some $\lambda,\Lambda>0$, where $\Omega_{0,R_{0}}$ is defined as in (1).
We assume $A^{\alpha\beta}_{ij}(x)\in C^{\mu}(\Omega_{0,R_{0}})$ for some
$\mu\in(0,1)$, and consider the system

$\left\\{\begin{aligned}
-\partial_{\alpha}\left(A^{\alpha\beta}_{ij}(x)\partial_{\beta}u_{j}(x)\right)&=0\quad\mbox{in
}\Omega_{0,R_{0}},\\\
A^{\alpha\beta}_{ij}(x)\partial_{\beta}u_{j}(x)\nu_{\alpha}(x)&=0\quad\mbox{on
}\Gamma_{+}\cup\Gamma_{-},\\\ \end{aligned}\right.$ (4.1)

for $i=1,\cdots,N$, where $\Gamma_{+},\Gamma_{-}$ are defined as in (1),
$\nu=(\nu_{1},\cdots,\nu_{n})$ denotes the unit normal vector on $\Gamma_{+}$
and $\Gamma_{-}$, pointing upward and downward respectively. We have the
following gradient estimate by essentially following the proof of Theorem 1.2
in [7].

######  Theorem 4.1.

Let $u\in H_{0}^{1}(\Omega_{0,R_{0}};\mathbb{R}^{N})$ be a solution to (4.1)
in dimension $n\geq 2$, with the coefficient $A^{\alpha\beta}_{ij}$ defined as
above. There exist positive constants $r_{0}$ and $C$ depending only on $n$,
$\lambda$, $\Lambda$, $R_{0}$, $\kappa$, $\mu$,
$\|A\|_{C^{\mu}(\Omega_{0,R_{0}})}$, $\|f\|_{C^{2}(\\{|x^{\prime}|\leq
R_{0}\\})}$ and $\|g\|_{C^{2}(\\{|x^{\prime}|\leq R_{0}\\})},$ such that

$|\nabla u(x_{0})|\leq
C\|u\|_{L^{\infty}(\Omega_{0,R_{0}})}\left(\varepsilon+|x_{0}^{\prime}|^{2}\right)^{-1/2},$
(4.2)

for all $\varepsilon\in(0,1)$, $x_{0}\in\Omega_{0,r_{0}}$.

######  Remark 4.2.

The elliptic systems we have considered include the linear systems of
elasticity: $n=N$, and the coefficients $A^{\alpha\beta}_{ij}$ satisfy

$A^{\alpha\beta}_{ij}=A^{\beta\alpha}_{ji}=A^{i\beta}_{\alpha j},$

and for all $n\times n$ symmetric matrices $\\{\xi_{\alpha}^{i}\\}$,

$\lambda|\xi|^{2}\leq
A^{\alpha\beta}_{ij}\xi_{\alpha}^{i}\xi_{\beta}^{j}\leq\Lambda|\xi|^{2}.$

###### Proof of Theorem 4.1.

Let $u\in H^{1}(\Omega_{0,R_{0}};\mathbb{R}^{N})$ be a solution to (4.1). We
perform the changes of variables (3.12) and (3.13). For any $1\leq i,j\leq N$,
we define

$B_{ij}^{\alpha\beta}(z)=\frac{(\partial_{y}z)(A^{\alpha\beta}_{ij})(\partial_{y}z)^{t}}{\det(\partial_{y}z)},$

and let $v(z)=u(x)$. Then $v$ satisfies

$\left\\{\begin{aligned}
-\partial_{\alpha}\left(B^{\alpha\beta}_{ij}(z)\partial_{\beta}v_{j}(z)\right)&=0\quad\mbox{in
}Q_{1,\delta},\\\ B^{n\beta}_{ij}(z)\partial_{\beta}v_{j}(z)&=0\quad\mbox{on
}\\{z_{n}=-\delta\\}\cup\\{z_{n}=\delta\\},\end{aligned}\right.$

for $i=1,\cdots,N$, where $Q_{s,t}$ is defined as in (3.3). As in the proof of
Theorem 1.1, we can show that

$\|B^{\alpha\beta}_{ij}\|_{L^{\infty}(Q_{1,\delta})}\leq
C\Lambda,\quad\|B^{\alpha\beta}_{ij}\|_{C^{\mu}(\bar{Q}_{1,\delta})}\leq C,$
$\int_{Q_{1,\delta}}B^{\alpha\beta}_{ij}(z)\partial_{\alpha}\varphi_{i}(z)\partial_{\beta}\varphi_{j}(z)\geq\frac{\lambda}{C}\int_{Q_{1,\delta}}|\nabla\varphi|^{2},\quad\forall\varphi\in
H_{0}^{1}(Q_{1,\delta};\mathbb{R}^{N}),$

where $C$ is a positive constant that depends only on $n$, $N$, $\mu$,
$R_{0}$, $\kappa$, $\lambda$, $\Lambda$, $\|A\|_{C^{\mu}}$, $\|f\|_{C^{2}}$
and $\|g\|_{C^{2}}$. Then we argue as in the proof of Theorem 1.1 to obtain

$|\nabla v(0)|\leq C\|v\|_{L^{\infty}(Q_{1,\delta})},$

which is (4.2) after reversing the changes of variables (3.12) and (3.13).
∎

## References

  * [1] H. Ammari, G. Ciraolo, H. Kang, H. Lee, and K. Yun.  Spectral analysis of the Neumann-Poincaré operator and characterization of the stress concentration in anti-plane elasticity.  Arch. Ration. Mech. Anal., 208(1):275–304, 2013.
  * [2] H. Ammari, B. Davies, and S. Yu.  Close-to-touching acoustic subwavelength resonators: eigenfrequency separation and gradient blow-up.  Multiscale Model. Simul., 18(3):1299–1317, 2020.
  * [3] H. Ammari, H. Kang, H. Lee, J. Lee, and M. Lim.  Optimal estimates for the electric field in two dimensions.  J. Math. Pures Appl. (9), 88(4):307–324, 2007.
  * [4] H. Ammari, H. Kang, and M. Lim.  Gradient estimates for solutions to the conductivity problem.  Math. Ann., 332(2):277–286, 2005.
  * [5] I. Babuška, B. Andersson, P.J. Smith, and K. Levin.  Damage analysis of fiber composites. I. Statistical analysis on fiber scale.  Comput. Methods Appl. Mech. Engrg., 172(1-4):27–77, 1999.
  * [6] E. Bao, Y.Y. Li, and B. Yin.  Gradient estimates for the perfect conductivity problem.  Arch. Ration. Mech. Anal., 193(1):195–226, 2009.
  * [7] E. Bao, Y.Y. Li, and B. Yin.  Gradient estimates for the perfect and insulated conductivity problems with multiple inclusions.  Comm. Partial Differential Equations, 35(11):1982–2006, 2010.
  * [8] J.G. Bao, H.G. Li, and Y.Y. Li.  Gradient estimates for solutions of the Lamé system with partially infinite coefficients.  Arch. Ration. Mech. Anal., 215(1):307–351, 2015.
  * [9] J.G. Bao, H.G. Li, and Y.Y. Li.  Gradient estimates for solutions of the Lamé system with partially infinite coefficients in dimensions greater than two.  Adv. Math., 305:298–338, 2017.
  * [10] E. Bonnetier and F. Triki.  Pointwise bounds on the gradient and the spectrum of the Neumann-Poincaré operator: the case of 2 discs.  In Multi-scale and high-contrast PDE: from modelling, to mathematical analysis, to inversion, volume 577 of Contemp. Math., pages 81–91. Amer. Math. Soc., Providence, RI, 2012.
  * [11] E. Bonnetier and F. Triki.  On the spectrum of the Poincaré variational problem for two close-to-touching inclusions in 2D.  Arch. Ration. Mech. Anal., 209(2):541–567, 2013.
  * [12] E. Bonnetier and M. Vogelius.  An elliptic regularity result for a composite medium with “touching” fibers of circular cross-section.  SIAM J. Math. Anal., 31(3):651–677, 2000.
  * [13] B. Budiansky and G. F. Carrier.  High Shear Stresses in Stiff-Fiber Composites.  Journal of Applied Mechanics, 51(4):733–735, 1984.
  * [14] Yves Capdeboscq and S.C. Yang Ong.  Quantitative Jacobian determinant bounds for the conductivity equation in high contrast composite media.  Discrete Contin. Dyn. Syst. Ser. B, 25(10):3857–3887, 2020.
  * [15] H. Dong and H.G. Li.  Optimal estimates for the conductivity problem by Green’s function method.  Arch. Ration. Mech. Anal., 231(3):1427–1453, 2019.
  * [16] H. Dong and H. Zhang.  On an elliptic equation arising from composite materials.  Arch. Ration. Mech. Anal., 222(1):47–89, 2016.
  * [17] D. Gilbarg and N. Trudinger.  Elliptic partial differential equations of second order.  Classics in Mathematics. Springer-Verlag, Berlin, 2001.
  * [18] Y. Gorb.  Singular behavior of electric field of high-contrast concentrated composites.  Multiscale Model. Simul., 13(4):1312–1326, 2015.
  * [19] H. Kang, M. Lim, and K. Yun.  Asymptotics and computation of the solution to the conductivity equation in the presence of adjacent inclusions with extreme conductivities.  J. Math. Pures Appl. (9), 99(2):234–249, 2013.
  * [20] H. Kang, M. Lim, and K. Yun.  Characterization of the electric field concentration between two adjacent spherical perfect conductors.  SIAM J. Appl. Math., 74(1):125–146, 2014.
  * [21] J. B. Keller.  Stresses in Narrow Regions.  Journal of Applied Mechanics, 60(4):1054–1056, 1993.
  * [22] J. Kim and M. Lim.  Electric field concentration in the presence of an inclusion with eccentric core-shell geometry.  Math. Ann., 373(1-2):517–551, 2019.
  * [23] H.G. Li.  Asymptotics for the Electric Field Concentration in the Perfect Conductivity Problem.  SIAM J. Math. Anal., 52(4):3350–3375, 2020.
  * [24] H.G. Li, Y.Y. Li, and Z. Yang.  Asymptotics of the gradient of solutions to the perfect conductivity problem.  Multiscale Model. Simul., 17(3):899–925, 2019.
  * [25] H.G. Li, F. Wang, and L. Xu.  Characterization of electric fields between two spherical perfect conductors with general radii in 3D.  J. Differential Equations, 267(11):6644–6690, 2019.
  * [26] Y.Y. Li and L. Nirenberg.  Estimates for elliptic systems from composite material.  Comm. Pure Appl. Math., 56(7):892–925, 2003.
  * [27] Y.Y. Li and M. Vogelius.  Gradient estimates for solutions to divergence form elliptic equations with discontinuous coefficients.  Arch. Ration. Mech. Anal., 153(2):91–151, 2000.
  * [28] M. Lim and S. Yu.  Stress concentration for two nearly touching circular holes.  arXiv:1705.10400.
  * [29] M. Lim and K. Yun.  Blow-up of electric fields between closely spaced spherical perfect conductors.  Comm. Partial Differential Equations, 34(10-12):1287–1315, 2009\.
  * [30] X. Markenscoff.  Stress amplification in vanishingly small geometries.  Computational Mechanics, 19(1):77–83, 1996.
  * [31] K. Yun.  Estimates for electric fields blown up between closely adjacent conductors with arbitrary shape.  SIAM J. Appl. Math., 67(3):714–730, 2007.
  * [32] K. Yun.  Optimal bound on high stresses occurring between stiff fibers with arbitrary shaped cross-sections.  J. Math. Anal. Appl., 350(1):306–312, 2009.
  * [33] K. Yun.  An optimal estimate for electric fields on the shortest line segment between two spherical insulators in three dimensions.  J. Differential Equations, 261(1):148–188, 2016.


# Optimization of Hydrogen Blending in Natural Gas Networks for Carbon
Emissions Reduction

Mo Sodwatana1, Saif R. Kazi2, Kaarthik Sundar3, and Anatoly Zlotnik2 *This
study was funded by the U.S. Department of Energy’s Advanced Grid Modeling
(AGM) projects “Joint Power System and Natural Gas Pipeline Optimal Expansion”
and “Dynamical Modeling, Estimation, and Optimal Control of Electrical Grid-
Natural Gas Transmission Systems”. Research conducted at Los Alamos National
Laboratory is done under the auspices of the National Nuclear Security
Administration of the U.S. Department of Energy under Contract No.
89233218CNA000001.1Mo Sodwatana is a Ph.D. student in the Department of Energy
Science and Engineering, Stanford University, Stanford, CA 94305, USA
<EMAIL_ADDRESS>R. Kazi and Anatoly Zlotnik are in the Applied
Mathematics & Plasma Physics Group, and Kaarthik Sundar is in the Information
Systems & Modeling Group, Los Alamos National Laboratory, Los Alamos, NM
87545, USA<EMAIL_ADDRESS>
###### Abstract

We present an economic optimization problem for allocating the flow of natural
gas and hydrogen blends through a large-scale transportation pipeline network.
Physical flow of the gas mixture is modeled using a steady-state relation
between pressure decrease and flow rate, which depends on mass concentration
of the constituents as it varies by location in the network. The objective
reflects the economic value provided by the system, accounting for delivered
energy in withdrawn flows, the cost of natural gas and hydrogen injections,
and avoided carbon emissions. The problem is solved subject to physical flow
equations, nodal balance and mixing laws, and engineering inequality
constraints. The desired energy delivery rate and minimum hydrogen
concentration can be specified as upper and lower bound values, respectively,
of inequality constraints, and we examine the sensitivity of the physical
pressure and flow solution to these parameters for two test networks. The
results confirm that increasing hydrogen concentration requires greater energy
expended for compression to deliver the same energy content, and the
formulation could be used for valuation of the resulting mitigation of carbon
emissions.

##  I INTRODUCTION

The United States bulk power grid increasingly relies on natural gas (NG)
fueled generation and renewable energy, while planning for reduced dependence
on fossil fuels. The blending of hydrogen (H2) generated using clean energy
into natural gas delivery systems is proposed to support this energy
transition while using the capital investments in existing pipelines for their
entire lifespan [1]. The use of NG-H2 mixtures has the potential to reduce
greenhouse gas emissions from power generation as well as residential end-use
[2]. The ability to inject hydrogen produced using renewable energy into gas
pipelines will also provide operational flexibility and storage capacity for
the power grid, and takes advantage of an existing infrastructure that would
otherwise be stranded with increasing electrification [3]. There are however
system integration issues and associated costs of such hydrogen utilization
that affect pipeline transport, which must be considered using nontrivial
extensions of methods developed for homogeneous gas transport [4].

Blends of hydrogen and natural gas exhibit different physical flow properties
than pure hydrogen or natural gas. Several studies have examined the effects
of the properties of the blended gas, such as density, viscosity, phase
interactions, and energy densities, on the pipeline network and end-use
applications [5, 6]. In particular, when hydrogen is blended into natural gas
injected into a pipeline and without changes to operating setpoints, there is
a reduction in transported gas pressure and energy delivered downstream.
Whereas hydrogen has a higher calorific value by mass than methane, it has a
lower density at given pressures. In a recent case study, the energy quantity
transported for the same pressure ratio was reduced by 4% and 14% for a 10%
and a 40% hydrogen blend, respectively, assuming equipment compatibility [7].
Moreover, due to the pressure reduction, the energy required for compression
increases by 7% and 30% for the respective blend percentages. Other recent
technical studies focus on the feasibility and safety of injection and
pressurization [8].

For natural gas transport, canonical problems utilize steady-state
optimization to evaluate capacity [9] and determine economically optimal
allocation [10]. In these problems, the gas is assumed to be chemically
homogeneous. However, composition tracking may be important when a pipeline
has multiple receipt points for species with different calorific values [11].
Assuming that technical issues related to sealing, compression, and end-use
appliances can be resolved and that engineering limitations on quantities such
as minimum and maximum pressures and compressor energy can be specified [12],
the design, operation, and economics of H2-NG blend pipelines must be
considered. Energy from hydrogen gas produced by electrolysis is much more
costly today than prevailing natural gas prices [13]. Quantifying the cost of
avoided carbon emissions resulting from H2 blending into NG pipelines, as well
as the resulting change in energy delivery capacity, is therefore of
particular interest.

In this study, we formulate an economic optimization problem for determining a
feasible flow allocation that maximizes economic benefit for users of a
pipeline system that transports a blend of two significantly different gases.
The formulation is effectively a single auction market mechanism in which
suppliers offer natural gas or hydrogen at a given price per mass flow rate,
and consumers bid for deliveries at a given price per energy content. In
addition, consumers can individually provide a bid for carbon mitigation at a
given price per mass of carbon dioxide not emitted at that location. This
gives a specific value to the amount of energy delivered in the form of
hydrogen that otherwise would arrive in the form of natural gas. Control
variables in the formulation include optimized injections of pure NG or H2,
withdrawal rates of the mixture, and compressor operating setpoints. The
optimization is solved subject to physical flow equations, nodal balance and
mixing laws, and engineering inequality constraints.

The rest of the manuscript is structured as follows. In Section II, we specify
network modeling and physical flow equations for a pipeline that transports a
spatially inhomogeneous mixture of two gases arising from distributed
injections of these constituents. In Section III, we define an objective
function and additional inequality constraints that include engineering
limitations, and formulate the optimization problem. We present the results of
sensitivity analyses performed for two case studies in Section IV, discuss the
implications and potential follow-on studies in Section V, and briefly
conclude in Section VI.

##  II Heterogeneous Gas Pipeline Network Model

We consider a gas pipeline network that is represented using a directed graph
with junctions $j\in\mathcal{V}$ that are connected by pipes
$(i,j)\in\mathcal{E}$. Another set $\mathcal{C}$ of node-connecting elements
is the collection of compressors, which are used to boost gas pressure.
Following previously developed notation [10], we also specify a set
$\mathcal{G}$ of gNodes, where each gNode $m\in\mathcal{G}$ represents a user
of the pipeline that is associated to a physical node $j(m)\in\mathcal{V}$.
Physical nodes with associated gNodes can be either supply or withdrawal
nodes, because there cannot, for example, be one gNode that injects hydrogen
while another makes a withdrawal at the same physical node. The sets of
components that constitute the network are denoted by

$j\in\mathcal{V}$ | set of all physical nodes,
---|---
$(i,j)\in\mathcal{E}$ | set of edges representing pipes,
$(i,j)\in\mathcal{C}$ | set of edges with compressors,
$j\in\mathcal{V}_{s}$ | set of slack physical nodes, subset of $\mathcal{V}$,
$m\in\mathcal{G}_{s}^{H_{2}}$ | set of gNodes that inject hydrogen,
$m\in\mathcal{G}_{s}^{NG}$ | set of gNodes that inject natural gas,
$m\in\mathcal{G}_{d}$ | set of gNodes that withdraw gas.

The physical state of the network is defined by the gas pressure $P_{j}$ and
mass fraction of hydrogen $\gamma_{j}$ for each junction $j\in\mathcal{V}$,
and the total mass flow $\phi_{ij}$ and mass fraction of hydrogen
$\gamma_{ij}$ on each pipe $(i,j)\in\mathcal{E}$. Our key structural
assumption is that flow directions are fixed _a priori_. We suppose that the
control variables available to the operator and users of the pipeline system
are the compressor ratio $\alpha_{ij}$ of each compressor
$(i,j)\in\mathcal{C}$, the hydrogen gas supply $s_{m}^{H_{2}}$ at gNode
$m\in\mathcal{G}_{s}^{H_{2}}$, the natural gas supply $s_{m}^{NG}$ at gNode
$m\in\mathcal{G}_{s}^{NG}$, and the withdrawal flow $d_{m}$ of the gas mixture
from gNode $m\in\mathcal{G}_{d}$. Diagrams of test networks used in our study
are found in Figs. 1 and 6. We now describe the relations between these
quantities, and modeling parameters, which represent the function of the
network.

###  II-A Pipe Equations

We suppose that in steady-state flow, the hydrogen concentration is uniform
along each edge $(i,j)$. We use the Weymouth equation to model the relation
between pressures at the endpoints of the pipe and the flow through it [14].
This equation is derived from the momentum conservation law in the Euler
equations for one-dimensional turbulent flow in a pipe. The relation is given
as

$P_{i}^{2}-P_{j}^{2}=\frac{\lambda_{ij}L_{ij}}{D_{ij}A_{ij}^{2}}V_{ij}\phi_{ij}\left|\phi_{ij}\right|\quad\forall(i,j)\in\mathcal{E}$
(1a)

where here $V_{ij}$ is the speed of sound of the blended gas,

$V_{ij}=\gamma_{ij}a_{H_{2}}^{2}+(1-\gamma_{ij})a_{NG}^{2}\quad\forall(i,j)\in\mathcal{E}.$
(1b)

We approximate the state equation using the ideal gas law, so that partial
pressures are additive. The sound speed of the blend is approximated as the
linear combination of the wave speeds $a_{NG}$ and $a_{H_{2}}$ [m/s] of
natural gas and hydrogen, respectively. We use $a_{NG}\approx 370m/s$ and
$a_{H_{2}}\approx 1090m/s$ in our case studies. Using the ideal gas
approximation significantly simplifies our exposition of the economic
optimization presented here. The assumptions used to obtain equations
(1a)-(1b) can in principle be relaxed to extend the results to the regime of
non-ideal gases, which better approximates the conditions of gas
transportation pipelines [15].

###  II-B Nodal Equations

At every physical node $j$, the net mass flow through the node and the net
injection into the node must be balanced. We impose mass balance equations on
natural gas and hydrogen, given by

$\begin{split}(1-\gamma_{j})\sum_{k\in\partial_{j}^{-}}\phi_{jk}-\sum_{i\in\partial_{j}^{+}}(1-\gamma_{ij})\phi_{ij}\qquad\qquad\\\
\qquad\qquad=\sum_{m\in\partial_{j}^{g}}s_{m}^{NG}-(1-\gamma_{j})\sum_{m\in\partial_{j}^{g}}d_{m},\end{split}\quad\forall
j\in\mathcal{V},$ (2a)
$\\!\\!\gamma_{j}\sum_{k\in\partial_{j}^{-}}\phi_{jk}-\sum_{i\in\partial_{j}^{+}}\gamma_{ij}\phi_{ij}=\\!\sum_{m\in\partial_{j}^{g}}s_{m}^{H_{2}}-\gamma_{j}\\!\sum_{m\in\partial_{j}^{g}}d_{m}\quad\forall
j\in\mathcal{V}.$ (2b)

Here $\partial_{j}^{+}$ and $\partial_{j}^{-}$ are the sets of nodes connected
to node $j$ by outgoing and incoming edges, respectively. Adding together the
two equations (2a) and (2b) imposes the total nodal mass balance. To enforce
appropriate continuity in the concentration from the a node to an outgoing
edge, the following constraint is imposed at the node-to-edge interface to
enforce that the hydrogen concentration at the edge $(i,j)$ leaving node $i$
equals the concentration at node $i$:

$\gamma_{i}=\gamma_{ij}\quad\forall(i,j)\in\mathcal{E}.$ (3)

Finally, we suppose that at each slack node $j\in\mathcal{V}_{s}$, we specify
the pressure to be a nominal value $\sigma_{j}$:

$P_{j}=\sigma_{j}\quad\forall j\in\mathcal{V}_{s}.$ (4)

Standard boundary conditions for flow network simulation, including gas
pipelines, require at least one slack node for well-posedness, and we use this
convention here.

###  II-C Compressor Modeling

Gas transmission pipelines are constructed with compressor stations, which are
complex facilities that may have multiple compressor machines. For the purpose
of large-scale system modeling, we suppose that the action of a compressor
station $(i,j)\in\mathcal{C}$ is aggregated as a pressure boost ratio
$\alpha_{ij}$, which acts as

$P_{j}^{2}=\alpha_{ij}^{2}P_{i}^{2}\quad\forall(i,j)\in\mathcal{C}.$ (5)

We consider the power used for gas compression as an important factor in our
study, because increasing hydrogen fraction increases the amount of
compression work required to transport a given amount of energy in the form of
the gas blend. Following standard practice [16], the power $W_{c}$ used to
drive the compressor is formulated as

$W_{c}=\left(\frac{286.76\cdot(\kappa_{ij}-1)\cdot
T}{G_{ij}{\kappa{ij}}}\right)\left(\alpha_{ij}^{m}-1\right)\left|\phi_{ij}\right|,\quad\forall(i,j)\in\mathcal{C}$
(6)

where $m=(\kappa_{ij}-1)/\kappa_{ij}$. Here we use $\kappa_{ij}$ to denote the
specific heat capacity ratio for the mixed gas and $G_{ij}$ to denote the
specific gravity ratio. We approximate these ratios for the blend using linear
combinations of the specific ratios of each gas by

$\displaystyle\kappa_{ij}=\kappa_{H_{2}}\gamma_{ij}+\kappa_{NG}(1-\gamma_{ij})\quad\forall(i,j)\in\mathcal{C},$
(7a) $\displaystyle
G_{ij}=G_{H_{2}}\gamma_{ij}+G_{NG}(1-\gamma_{ij})\quad\forall(i,j)\in\mathcal{C}.$
(7b)

In this study, we use $\kappa_{NG}=1.304$, $\kappa_{H_{2}}=1.405$,
$G_{NG}=0.5537$, and $G_{H_{2}}=0.0696$. $T$ [K] is the compressor suction
temperature at which gas is being compressed, and is kept at 288.7 K in our
simulations. Here $\phi_{ij}$ [kg/s] is the flow rate along the compressor
edge $(i,j)$.

###  II-D Carbon Emissions Offset

A key innovation of our study is to include carbon emissions mitigation based
on the value of carbon displacement that is submitted as part of a bid by each
consumer of energy to the optimization-based auction market mechanism. The
emissions offset $E_{m}$ of a specific consumer $m\in\mathcal{G}_{d}$ is
formulated as

$E_{m}=d_{m}\gamma_{j(m)}\cdot\frac{R_{H_{2}}}{R_{NG}}\cdot\zeta_{NG},$ (8)

where $\zeta_{NG}$ is the ratio of molecular weight of carbon dioxide to
natural gas, and is approximately 44/18. Here the constants $R_{NG}$ and
$R_{H_{2}}$ are the calorific values for natural gas and hydrogen,
respectively, for which we use values of $R_{H_{2}}=141.8MJ/kg$ and
$R_{NG}=44.2MJ/kg$. The emissions term (8) denotes the amount of carbon
dioxide [kg/s] that was not produced because hydrogen was burned instead of
natural gas to produce a given amount of energy using the delivered flow
$d_{m}$ [kg/s] at concentration $\gamma_{j(m)}$. Here, we examine the
emissions from the energy replaced, rather than the mass replaced, from
burning hydrogen, and therefore we include the ratio of calorific value of
hydrogen to natural gas.

##  III Optimization Formulation

We now describe the economic optimization formulation for constrained
heterogeneous gas transport. In addition to the physical and network modeling
in Section II, we define an objective function, inequality constraints that
specify engineering limitations, and parameters that constitute an economic
bid structure.

###  III-A Economic Value Objective Function

Our objective is to maximize the economic value of transporting gas between
suppliers of constituent gases as commodities and consumers who purchase
energy and carbon emissions offsets. The economic value produced by the
pipeline is the sum over all gNodes $m$ of payments by consumers for delivered
blended gas minus purchases from suppliers of received gas constituents.
Suppliers place offer prices for natural gas, $c_{m}^{NG}$ [$/kg], and
hydrogen, $c_{m}^{H_{2}}$ [$/kg], at the supplying gNodes, while off-takers
place bids for the energy contente of the blended gas, $c_{m}^{d}$ [$/MJ], at
the off-taking gNodes. Off-takers can also place a value on the carbon
emissions avoided, $c_{m}^{CO_{2}}$ [$/kg]. The offer price is in terms of
mass while the bid price is in energy units, which represents how suppliers
and off-takers consider the value of gas differently in a market that includes
the cost of emissions. We also include the work done to compress gas by each
compressor $c\in\mathcal{C}$, which is denoted by $W_{c}$ in the economic
value formulation. The economic value objective function is then expressed as

$\begin{split}J_{EV}=\sum_{m\in
G}\bigg{(}c_{m}^{d}d_{m}(R_{H_{2}}\gamma_{j(m)}+R_{NG}(1-\gamma_{j(m)}))-\\\
c_{m}^{H_{2}s}s_{m}^{H_{2}}-c_{m}^{NGs}s_{m}^{NG}+c_{m}^{CO_{2}}E_{m}\bigg{)}-\eta\sum_{c\in
C}W_{c}\end{split}$ (9)

where $\eta$ [$/kw-s] is a conversion factor used to define the economic cost
of applied compressor power. $R_{NG}$ and $R_{H_{2}}$ are the calorific values
for natural gas and hydrogen, respectively. The calorific value of the blended
gas is the linear combination of the calorific values of the two gases with
respect to $\gamma_{j(m)}$, the fractional concentration of hydrogen at $j$.
In our computational case studies, $\eta=$$0.13/3600kw-s, and
$R_{H_{2}}=141.8MJ$/kg and $R_{NG}=44.2MJ$/kg as specified in Section II-D.
$d_{m}$ [kg/s] is the mass flow rate of delivered blended gas while
$s_{m}^{NG}$ and $s_{m}^{H_{2}}$ [kg/s] are the respective mass flow rates of
natural gas and hydrogen at the supply gNodes.

###  III-B Pressure, Compressor, and Concentration Limits

We suppose that minimum pressure limits as well as minimum and maximum
hydrogen concentration limits may be specified at each node:

$\displaystyle P_{j}^{min}\leq P_{j}\quad$ $\displaystyle\forall
j\in\mathcal{V},$ (10a)
$\displaystyle\gamma_{j}^{min}\leq\gamma_{j}\leq\gamma_{j}^{max}\quad$
$\displaystyle\forall j\in\mathcal{V}.$ (10b)

In addition, the discharge node of each compressor has a maximum allowable
operating pressure, and each compressor has a maximum boost ratio:

$\displaystyle\alpha_{ij}P_{i}\leq P_{ij}^{max}\quad$
$\displaystyle\forall(i,j)\in\mathcal{C},$ (11a) $\displaystyle
1\leq\alpha_{ij}\leq\alpha^{max}_{ij}\quad$
$\displaystyle\forall(i,j)\in\mathcal{C}.$ (11b)

###  III-C Supply and Demand Limits

The supplies of natural gas and hydrogen at gNodes $m$ are positive, and are
constrained by upper bounds of the form

$\displaystyle 0\leq s_{m}^{NG}\leq s_{m}^{max,NG}\quad$ $\displaystyle\forall
m\in\mathcal{G}_{s}^{NG},$ (12a) $\displaystyle 0\leq s_{m}^{H_{2}}\leq
s_{m}^{max,H_{2}}\quad$ $\displaystyle\forall m\in\mathcal{G}_{s}^{H_{2}}.$
(12b)

We suppose that outflows $d_{m}$ from the network are positive. The upper
bound of the demand is in terms of energy, with calorific values as conversion
from mass flow to energy flow. Energy deliveries are either optimized or
fixed, corresponding to gNode sets $\mathcal{G}_{d,o}$ and
$\mathcal{G}_{d,f}$, with constraints of form

$\displaystyle
d_{m}\left(R_{H_{2}}\gamma_{j(m)}+R_{NG}(1-\gamma_{j(m)})\right)\leq
g_{m}^{max}\quad$ $\displaystyle\forall m\in\mathcal{G}_{d,o},$ (13a)
$\displaystyle
d_{m}\left(R_{H_{2}}\gamma_{j(m)}+R_{NG}(1-\gamma_{j(m)})\right)=\bar{g}_{m}\quad$
$\displaystyle\forall m\in\mathcal{G}_{d,f}.$ (13b)

Combining equations (1a)-(13) yields the economic heterogeneous gas transport
optimization problem:

$\begin{array}[]{ll}\\!\\!\\!\\!\mathrm{max}&J_{EV}\triangleq\text{max
economic value objective }\eqref{eq:obj}\\\ \\!\\!\\!\\!\text{s.t.}&\text{NG
flow balance }\eqref{eq:ngflowbalance}\\\ &\text{H\textsubscript{2} flow
balance }\eqref{eq:h2flowbalance}\\\ &\text{pressure balance
}\eqref{eq:pipe_equations}\\\ &\text{compression }\eqref{eq:comp_boost}\\\
&\text{slack pressure }\eqref{eq:slack_pressure}\\\ &\text{fixed consumption
}\eqref{eq:demand_fixed}\\\ &\text{concentration continuity
}\eqref{eq:continuity}\\\ &\text{pressure limits
}\eqref{eq:minpressure},\eqref{eq:maxpressure}\\\ &\text{compressor boost
limits }\eqref{eq:complimits}\\\ &\text{H\textsubscript{2} concentration
limits }\eqref{eq:conclimits}\\\ &\text{supply limits
}\eqref{eq:supplylimits}\\\ &\text{demand limits
}\eqref{eq:demand_opt}\end{array}$ (14)

Observe that parameters
$(g_{m}^{max},c_{m}^{d},c_{m}^{CO_{2}},\gamma_{j(m)}^{min})$ can be
interpreted as a market bid by gas consumers, to include a quantity (MJ) and
price ($/MJ) of delivered energy, a price ($/Kg) of carbon offset, and minimum
hydrogen mass fraction (%). The decision variables for problem (14) and their
units are listed in Table I.

Variables | Set | Units (SI)
---|---|---
$s_{m}^{H_{2}}$ | $m\in\mathcal{G}_{s}^{H_{2}}$ | kg$\cdot$ s-1
$s_{m}^{NG}$ | $m\in\mathcal{G}_{s}^{NG}$ | kg$\cdot$ s-1
$d_{m}$ | $m\in\mathcal{G}_{d}$ | kg$\cdot$ s-1
$\alpha_{ij}$ | $(i,j)\in\mathcal{C}$ | -
$\phi_{ij}$ | $(i,j)\in\mathcal{E}$ | kg$\cdot$ s-1
$\gamma_{ij}$ | $(i,j)\in\mathcal{E}$ | -
$\gamma_{j}$ | $j\in\mathcal{V}$ | -
$P_{j}$ | $j\in\mathcal{V}$ | kg$\cdot$ m-1$\cdot$ s-2 (Pa)
TABLE I: Collection of optimization variables for Problem (14).

##  IV Case Studies

We examine solutions to the problem (14) using two test networks – a single
pipe with a compressor at the start (Fig. 1) and an 8-node network with a loop
and three compressors (Fig. 6. Problem 14 is solved in the Julia programming
language v1.7.2 using the JuMP package v1.1.0 [17], which is an optimization
modeling toolkit. We use IPOPT v1.0.2, a large scale optimization package
[18], as the nonlinear program solver. The case studies are computed on an AMD
EPYC 7742 64-core processor with 16 GB of RAM. The solve times for both
networks are mere seconds, and sensitivity analyses involving batches of
optimization problems require approximately five minutes or less, depending on
network size and other factors.

Figure 1: Schematic of a single pipe network with compressor C1, natural gas
and hydrogen injection at gNodes, S1 and S2, respectively, and a withdrawal
gNode D1. The pipe and compressor characteristics are shown, as well as offer
and bid prices at the supply and demand node, respectively.

###  IV-A Non-Dimensionalization and Rescaling

We non-dimensionalize the governing equations prior to solving problem (14) in
order to avoid numerical issues [15]. In addition, we re-scale (1a) because
the wave speed $V$ in blended gas is orders of magnitude larger than other
variables in the equation. Let $\bar{P}=P/P_{0}$, $\bar{L}=L/l_{0}$,
$\bar{D}=D/l_{0}$, $\bar{A}=A/A_{0}$, and
$\bar{\phi}=\phi/\phi_{0}=\phi/(\rho_{0}u_{0}A_{0})$. Equation (1a) then
becomes

$\displaystyle\bar{P}_{i}^{2}-\bar{P}_{j}^{2}=\frac{\lambda_{ij}\bar{L}_{ij}}{\bar{D}_{ij}\bar{A}_{ij}^{2}}\bar{V}_{ij}\bar{\phi_{ij}}\left|\bar{\phi}_{ij}\right|\cdot\frac{u_{0}^{2}}{a_{0}^{2}}\cdot\quad\forall(i,j)\in\mathcal{E},$
(15a)
$\displaystyle\bar{V}_{ij}\triangleq\frac{V_{ij}(\gamma_{ij})}{a_{0}^{2}}\quad\forall(i,j)\in\mathcal{E}.$
(15b)

We use a wave speed re-scaling factor of $a_{0}=635.06$ m/s. For the single
pipe and 8-node cases, the nominal pressures are $P_{0}=5$ MPa and
$P_{0}=3.04$ MPa, respectively. In addition, the nominal length, area, density
and velocity for both cases are as follows: $l_{0}=5000$ m, $A_{0}=1$ m2,
$\rho_{0}=P_{0}/a_{0}^{2}$, and $u_{0}=\lceil a_{0}\rceil/300$, where $a_{0}$
is the speed geometric mean, used in the re-scaling factor. We compute $a_{0}$
as $a_{0}=\sqrt{a_{NG}\cdot a_{H_{2}}}$ where $a_{NG}=370$ m/s and
$a_{H_{2}}=1090$ m/s is the wave speeds of NG and H2, respectively.
$a_{NG}=\sqrt{RT/M_{NG}}$ and $a_{H_{2}}=\sqrt{RT/M_{H_{2}}}$. Here, $R=8.314$
J/mol/K is the universal gas constant. $M_{NG}=0.01737$ kg/mol and
$M_{H_{2}}=0.002016$ kg/mol are the molecular mass of NG and H2, respectively.

###  IV-B Single Pipe

For the single pipe test network, we examine how the solution changes due to
1) variation of the constraint bound $g_{m}^{max}$ for maximum energy demand
(13a); 2) a similar analysis with fixed nodal values of carbon emissions
offsets $c_{m}^{CO_{2}}E_{m}$ added to the objective; 3) variation of the
constraint bound $\gamma_{j}^{min}$ for minimum hydrogen concentration (10b);
and 4) variation of the carbon emissions offset price $c_{m}^{CO_{2}}$. The
structure and characteristics of the single pipe network are described in
Figure 1, and the four sensitivity analyses and their results are described in
Sections IV-B1 through IV-B4.

####  IV-B1 Sensitivity analysis with respect to the maximum energy demand at
the withdrawal gNode, $g^{\max}_{D1}$, without carbon emissions offset values

Figure 2: Results of sensitivity analysis for the single pipe network in
Section IV-B1 are shown as the change in physical variables at the withdrawal
gNode D1, the shadow price at the physical node J3, and the compressor ratio
at C1. The maximum energy demand $g_{D1}^{max}$ is varied from 100 to 160
MJ/s. There is no carbon emissions offset at the withdrawal node
$(c_{D1}^{CO_{2}}\equiv 0)$.

Here we vary the constraint bound $g^{\max}_{D1}$ from 100 MJ/s to 160 MJ/s in
increments of 1 MJ/s, and solve problem (14) for each instance. The variations
in the solutions are shown in Figure 2, and can be divided into three regions
where the transitions arise from activation of new binding constraints. In the
first region, where $100\leq g^{\max}_{D1}\leq 115$ MJ/s, the energy demand
constraint binds because the system is able to meet the demand. The pressure
at each node is within the defined bounds and there is no work done by the
compressor. At 115 MJ/s, the pressure at the receiving node D1 hits the lower
bound of 3 MPa, and therefore the compressor must operate. In the region where
$115\leq g^{\max}_{D1}\leq 150$ MJ/s, the compressor ratio increases to meet
the rising energy demand. In the final region, where $150\leq g^{\max}_{D1}$
MJ/s, energy and mass flow delivered taper off to become constant regardless
of further increases in energy demand, with the compressor ratio binding at
its limit of 1.4.

The shadow prices of H2 and NG are shown to change with $g^{\max}_{D1}$ as
well in Figure 2. These quantities are computed by the solver, and are given
as the duals of the respective mass balance constraints at the physical
withdrawal node. The shadow price of the mixture is the linear combination of
the two shadow prices as weighted by mass concentration ratios. When a
constraint becomes active and a phase transition occurs, the shadow prices
adjust accordingly. As more power must be used for gas compression to deliver
energy, there is a steady increase in the shadow prices. Likewise, when the
energy demand can no longer be met, there is a hike in the shadow prices.

####  IV-B2 Sensitivity analysis with respect to the maximum energy demand at
the withdrawal gNode, $g^{\max}_{D1}$, with carbon emissions offset values

Figure 3: Results of sensitivity analysis for the single pipe network in
Section IV-B2 are shown as the change in physical variables at the withdrawal
gNode D1, the shadow price at the physical node J3, and the compressor ratio
at C1. The maximum energy demand $g_{D1}^{max}$ is varied from 100 to 160 MJ/s
when $c_{D1}^{CO_{2}}=\$0.055/kg$ at the withdrawal node.

Here we repeat the analysis in Section IV-B1, and include a carbon emissions
offset value where $c_{D1}^{CO_{2}}=\$0.055$/kg CO2 avoided. The changes in
the solution in this case are shown in Figure 3. The bound value
$g^{\max}_{D1}$ increases from 100 MJ/s to 160 MJ/s in increments of 1 MJ/s as
with the previous case. The solutions can be categorized into four regions.
Providing a value of reducing carbon emissions incentivizes end-use of
hydrogen, so the solution initially sees hydrogen injected at the maximum
allowable concentration of 0.1. In this first region of four, where
$g^{\max}_{D1}\leq 105$ MJ/s, the maximum energy demand is the binding
constraint. A transition occurs at 105 MJ/s where the pressure at junction J3
binds at the lower limit of 3 MPa. This pressure constraint is binding at a
lower maximum energy demand compared to the previous case, because hydrogen
blending reduces the pipeline pressure. For $105\leq g^{\max}_{D1}\leq 137$
MJ/s, there is compressor work done to ensure that the maximum energy delivery
constraint is being met while pressure does not fall below the lower limit. In
this region, the slope of the withdrawal flow rate becomes steeper because
hydrogen in the gas mixture requires more natural gas by mass to substitute
for the same energy. In region three, where $137\leq g^{\max}_{D1}\leq 150$
MJ/s, the solution switches to injecting pure NG to maintain allowable
pressure while meeting increasing energy demands. This behavior is contrary to
the previous case where the energy delivered tapers off. In the final region,
for $g^{\max}_{D1}\geq 150$ MJ/s, the energy delivered and withdrawal flow
become constant as the system becomes congested cannot supply energy at the
rate $g^{\max}_{D1}$.

Figure 4: Results of sensitivity analysis in Section IV-B3 with respect to the
minimum allowable hydrogen concentration at the withdrawal gNode, D1. Here,
the maximum energy demand is constant at 140 MJ/s.

####  IV-B3 Sensitivity analysis with respect to the minimum allowable H2
concentration, $\gamma_{J3}^{\min}$, at the withdrawal gNode with a carbon
emissions offset

This section explores the sensitivity of the solution to adjustment of the
minimum hydrogen concentration $\gamma_{J3}^{\min}$ at the physical withdrawal
node J3 corresponding to the withdrawal gNode D1. This scenario reflects
instances where a consumer is required to mitigate their carbon emissions by
bidding for a minimum H2 concentration at the delivery node. We incrementally
increase $\gamma_{J3}^{\min}$ from 0.0 to 0.1 at an increment of 0.025 while
$g_{D1}^{\max}$ is constant at 140 MJ/s. Figure 4 provides a summary of the
results. The solutions can be categorized in three regions. For
$0\leq\gamma_{J3}^{\min}\leq 0.053$, the energy demand $g_{D1}^{\max}$ can be
met, and the actual delivered concentration increases with the minimum bound
value. The compressor ratio must also increase to keep the line pressure high
enough to meet the minimum pressure bound at J3. For
$0.053\leq\gamma_{J3}^{\min}\leq 0.068$, The maximum pressure at J2 binds so
that the pipe is congested, and increasing the minimum hydrogen concentration
results in less energy delivered. At $\gamma_{J3}^{\min}\approx 0.068$, the
concentration binds at the maximum value. Observe that the primal and dual
solutions change in a non-monotone manner as concentration parameters change,
which is consistent with prior observations on transport of heterogeneous gas
blends.

####  IV-B4 Sensitivity analysis with respect to the carbon emissions offset
price $c_{m}^{CO_{2}}$ at the withdrawal gNode, D1

Figure 5: Results of sensitivity analysis in Section IV-B4 with respect to the
carbon emissions offset price $c_{m}^{CO_{2}}$ at the withdrawal gNode, D1.
Here, the maximum energy demand is constant at 140 MJ/s.

For this analysis, we vary the carbon emissions offset price $c_{m}^{CO_{2}}$
at the withdrawal gNode D1 from 0 to $0.09 $/kgCO2. The solutions can be
categorized into three regions of interest. The energy demand $g_{D1}^{\max}$
is set to 140 MJ/s, and the minimum concentration $\gamma_{J3}^{\min}$ is set
to 0. Figure 5 shows a summary of the results. When $c_{m}^{CO_{2}}\leq
0.055$, no H2 is injected into the system, and $\gamma_{J3}\equiv 0$. The
energy delivered is binding at the maximum, and the flow is pure NG. Though
the shadow price of H2 gradually increases to reflect its increased value, it
is still not sufficient for injection of H2 to add overall economic value for
the system. In the region $0.055\leq c_{m}^{CO_{2}}\leq 0.065$, several
transition occur. First, the value of carbon mitigation surpasses the value of
delivered energy, so that H2 starts to be injected into the system, and the
maximum pressure at node J2 binds so that the pipeline becomes constrained. As
a result, less energy is delivered than the request at gNode D1, the
withdrawal flow drops, and more hydrogen is injected into the system since its
carbon offset value continues to be greater than the value of delivered energy
at the offer price. Then in the region $c_{m}^{CO_{2}}\geq 0.065$, the maximum
H2 concentration constraint at the physical node J3 binds, and the only
subsequent change is an increase in the shadow price of H2 to reflect its
increasing carbon mitigation value.

###  IV-C 8-Node Network

Figure 6: Schematic of an 8-node network with one limited NG supplier at gNode
S1, a supplier of unlimited H2 at gNode S2, and a limited supplier of H2 at
S3. The limited supplier at S3 provides H2 at a relatively cheaper price than
the H2 supplier at S2. At D1 and D2 $c_{D1}^{CO_{2}}$ = $0.055/kg CO2 avoided.
D3 is a fixed demand node with $\bar{g}_{m}^{\max}$ = 100 MJ/s and
$c_{D1}^{CO_{2}}$ = 0. gNodes S1 and S2 share the same physical node J1.
gNodes D2 and D3 share J5 as the physical node. Pipe length, diameter,
friction factor varies.

We use the 8-node network shown in Figure 6 to demonstrate the behavior of
solutions to optimization problem (14) in the case of a more complex network
topology. First, we see that the computational implementation readily scales
to looped topologies and multiple edges and junctions. We also perform a
sensitivity analysis by varying the maximum energy demands $g_{D1}^{\max}$ and
$g_{D2}^{\max}$ at both withdrawal gNodes with a carbon emissions offset fixed
at $c_{D1}^{CO_{2}}=\$0.055$/kg and the minimum H2 concentrations fixed at
$\gamma_{J3}^{\min}=\gamma_{J5}^{\min}=0.05$. We increase $g_{D1}^{\max}$ and
$g_{D2}^{\max}$ from 120 MJ/s to 180 MJ/s in increments of 1 MJ/s, and present
the results seen at gNode D1 in Figure 7 and results for gNode D2 in Figure 8.
We observe two regions in the solution. In the first region where
$g_{D1}^{\max},g_{D2}^{\max}\leq 140$ MJ/s, the solution can continually meet
both $g_{D1}^{\max}$ and $g_{D2}^{\max}$. In the second region,
$g_{D1}^{\max},g_{D2}^{\max}\geq 140$ MJ/s, the energy and the flow delivered
at D2 gradually declines because the system becomes congested due to a
constraint on how much gas can be delivered to physical node J5. Because D2
shares a physical node with D3, a gas consumer with a fixed demand, the
solution has to reduce supply to D2. Meanwhile, D1 receives uninterrupted
supply of the blended gas because of its proximity to the injection site.

Figure 7: Plots of the changing physical characteristics at the withdrawal
gNode D1, the shadow price at the physical node J3, and the compressor ratio
at C2 from a sensitivity analysis $g_{D1}^{max}$ on the 8-node pipe network.
Here, $c_{D1}^{CO_{2}}.$ = $0.055/kg CO2 avoided at the withdrawal gNode.
Figure 8: Plots of the changing physical characteristics at the withdrawal
gNode D2, the shadow price at the physical node J5, and the compressor ratio
at C3 from a sensitivity analysis $g_{D1}^{max}$ on the 8-node pipe network.
Here, $c_{D1}^{CO_{2}}.$ = $0.055/kg CO2 avoided at the withdrawal gNode.

While the H2 concentration of the gas supplied to D1 is 0.1, that of the blend
delivered to D2 is in the neighborhood of 0.06 as a result of the mixing of
gases from pipes P3 and P4. Moreover, D3 is a fixed demand consumer with no
carbon emissions offset, thus the solution does not seek to deliver H2 to at
J5. Moreover, we observe that the shadow prices computed at the two physical
nodes differ, which reflects the difference in the locational values of
energy, the constituent gases, and carbon emissions mitigation.

##  V Discussion

We have performed extensive sensitivity analyses on the problem (14) applied
to a single pipe in order to demonstrate the intuitive changes to the solution
with respect to key parameters. In analysis IV-B1, given the pricing
structure, no hydrogen injection occurs because hydrogen has lower energy
content at the mixture’s common pressure, and thus there is no incentive for
H2 injection to appear in the solution. This scenario demonstrates the
behavior of the system as demand is increased, constraints bind, the pipeline
is conjested, and the shadow prices exhibit jumps and ramps to reflect how
energy value increases with the demand.

In analysis IV-B2, the first analysis is repeated but now with an incentive to
replace NG with H2 in order to mitigate carbon emissions. We see that when the
system is not congested, i.e., the pipeline is not constrained by pressure
and/or compressor constraints at both sending and receiving ends, H2 is
supplied to the maximum possible extent. As energy demand increases, the
pipeline becomes congested, and because the value of energy is greater than
the value of carbon mitigation, there is a transition to pure NG utilization.

In analysis IV-B3, the demand is fixed, and the setting reflects an increasing
H2 utilization mandate. Energy transport remains in the form of NG to the
extent possible. Eventually, delivery of energy at the desired level cannot be
maintained. Counter-intuitively, the shadow price of H2 jumps and then
gradually decreases between changes in the binding constraint set. This is
evidence of the nonlinearity and non-monotonicity of heterogeneous gas
transport. This also indicates that mandating H2 blending could lead to
complex market structures with unexpected outcomes.

Finally, the analysis IV-B4 provides an insightful result. There the energy
demand upper bound is given, and the carbon offset price is increased. We
observe increasing shadow prices of both NG and H2 with an increasing carbon
offset price, and the solutions show monotone change with respect to this
parameter. When the value of offsetting emissions becomes greater than the
value of delivered energy, the gas blend quickly transitions to H2 injection
at the maximum possible level subject to engineering constraints. Future
studies could examine this break point in H2 concentration injection to
quantify the appropriate “green premium” for H2 blending, and examine whether
this quantity can be analytically derived from the Karush-Kuhn-Tucker
conditions. Potentially, a locational carbon offset value could be derived,
similarly to the notion of locational marginal carbon intensity that was
proposed for power grid generation [19]. The analysis of the 8-node case in
Section IV-C shows that the valuations are very much location-dependent. The
problem (14) could be used as an optimization-based market mechanism for
locational valuation of the emissions reduction achieved with hydrogen
injection into natural gas pipelines. In such a market, suppliers could
provide commodity prices of natural gas and hydrogen, and consumers could
provide price and quantity bids for energy, as well as the amount they are
willing to pay for carbon emissions reduction.

##  VI Conclusions

We present an economic optimization problem for allocating the flow of natural
gas and hydrogen blends through transportation pipeline networks, accounting
for delivered energy in withdrawn flows, the cost of natural gas and hydrogen
injections, and avoided carbon emissions. We examine the sensitivity of the
physical and dual solution to several parameters to examine carbon mitigation
mechanisms. Future work could compare locational and global carbon pricing,
and analytically derive dual variables for pricing solutions.

## References

  * [1] Andrzej Witkowski, Andrzej Rusin, Mirosław Majkut, and Katarzyna Stolecka.  Analysis of compression and transport of the methane/hydrogen mixture in existing natural gas pipelines.  International Journal of Pressure Vessels and Piping, 166:24–34, 2018.
  * [2] Arun S. K. Raju, Alfredo Martinez-Moralez, Oren Lever, and Love Asiedu-Akrofi.  Hydrogen blending impacts study, 2022.  California Public Utilities Commission.
  * [3] Dennis Anderson and Matthew Leach.  Harvesting and redistributing renewable energy: on the role of gas and electricity grids to overcome intermittency through the generation and storage of hydrogen.  Energy policy, 32(14):1603–1614, 2004.
  * [4] Maciej Chaczykowski, Filip Sund, Paweł Zarodkiewicz, and Sigmund Mongstad Hope.  Gas composition tracking in transient pipeline flow.  Journal of Natural Gas Science and Engineering, 55:321–330, 2018\.
  * [5] Devinder Mahajan, Kun Tan, T. Venkatesh, Pradheep Kileti, and Clive R. Clayton.  Hydrogen blending in gas pipeline networks—a review.  Energies, 15(10):3582, 2022.
  * [6] Marc W. Melaina, Olga Antonia, and Michael Penev.  Blending hydrogen into natural gas pipeline networks: A review of key issues.  NREL Tech. Rep. No., NREL/TP-5600-51995 6112, 2013.
  * [7] Francis Bainier and Rainer Kurz.  Impacts of h2 blending on capacity and efficiency on a gas transport network.  In Turbo Expo: Power for Land, Sea, and Air, volume 58721, page V009T27A014. American Society of Mechanical Engineers, 2019.
  * [8] Sebastian Schuster, Hans Josef Dohmen, and Dieter Brillert.  Challenges of compressing hydrogen for pipeline transportation with centrifugal compressors.  pages 2504–4400, 2020.
  * [9] Thorsten Koch, Benjamin Hiller, Marc E. Pfetsch, and Lars Schewe.  Evaluating gas network capacities.  SIAM, 2015.
  * [10] A. Rudkevich and A. Zlotnik.  Locational Marginal Pricing of Natural Gas subject to Engineering Constraints.  In Proc. of the 50th Hawaii International Conference on System Sciences, pages 3092–3101, 2017.
  * [11] Falk M. Hante and Martin Schmidt.  Complementarity-based nonlinear programming techniques for optimal mixing in gas networks.  EURO journal on computational optimization, 7(3):299–323, 2019\.
  * [12] Dale D. Erickson, James Holbeach, Danny Golcznski, and Shane A. Morrissy.  The Importance of Tracking Hydrogen H2 in Complex Natural Gas Networks.  In Offshore Technology Conference, 2022.
  * [13] Boreum Lee, Hyun-Seok Cho, Heehyang Kim, Dongjun Lim, Wonchul Cho, Chang-Hee Kim, and Hankwon Lim.  Integrative techno-economic and environmental assessment for green h2 production by alkaline water electrolysis based on experimental data.  Journal of Environmental Chemical Engineering, 9(6):106349, 2021\.
  * [14] Roger Z Ríos-Mercado and Conrado Borraz-Sánchez.  Optimization problems in natural gas transportation systems: A state-of-the-art review.  Applied Energy, 147:536–555, 2015.
  * [15] Shriram Srinivasan, Kaarthik Sundar, Vitaliy Gyrya, and Anatoly Zlotnik.  Numerical solution of the steady-state network flow equations for a non-ideal gas.  arXiv preprint arXiv:2204.00071, 2022.
  * [16] E. Shashi Menon.  Gas pipeline hydraulics.  Crc Press, 2005.
  * [17] Iain Dunning, Joey Huchette, and Miles Lubin.  Jump: A modeling language for mathematical optimization.  SIAM review, 59(2):295–320, 2017.
  * [18] Andreas Wächter and Lorenz T. Biegler.  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.  Mathematical programming, 106(1):25–57, 2006.
  * [19] Pablo A. Ruiz and Aleksandr Rudkevich.  Analysis of marginal carbon intensities in constrained power networks.  In 2010 43rd Hawaii International Conference on System Sciences, pages 1–9. IEEE, 2010.