A Common Approach to Singular Perturbation and Homogenization I: Quasilinear
ODE Systems

Nikolay N. Nefedov (Moscow) and Lutz Recke (Berlin)

###### Abstract

We consider perodic homogenization of boundary value problems for quasilinear
second-order ODE systems in divergence form of the type

$a(x,x/\varepsilon,u(x),u^{\prime}(x))^{\prime}=f(x,x/\varepsilon,u(x),u^{\prime}(x))\mbox{
for }x\in[0,1].$

For small $\varepsilon>0$ we show existence of weak solutions
$u=u_{\varepsilon}$ as well as their local uniqueness for
$\|u-u_{0}\|_{\infty}\approx 0$, where $u_{0}$ is a given non-degenerate
solution to the homogenized boundary value problem, and we describe the rate
of convergence to zero for $\varepsilon\to 0$ of the homogenization error
$\|u_{\varepsilon}-u_{0}\|_{\infty}$. In particular, we show that this rate
depends on the smoothness of the maps $a(\cdot,y,u,u^{\prime})$ and
$f(\cdot,y,u,u^{\prime})$.

Our assumptions are, roughly speaking, as follows: The maps
$a,f:[0,1]\times\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$
are continuous, the maps $a(x,y,\cdot,\cdot)$ and $f(x,y,\cdot,\cdot)$ are
$C^{1}$-smooth, the maps $a(x,\cdot,u,u^{\prime})$ and
$f(x,\cdot,u,u^{\prime})$ are 1-periodic, and the maps $a(x,y,u,\cdot)$ are
strongly monotone and Lipschitz continuous uniformly with respect to $x$, $y$
and bounded $u$. Neither global solution uniqueness is supposed nor
$W^{2,2}$-regularity of $u_{0}$.

The main tool of the proofs is an abstract result of implicit function theorem
type which in the past has been applied to singularly perturbed nonlinear ODEs
and elliptic and parabolic PDEs and, hence, which permits a common approach to
existence and local uniqueness results for singularly perturbed problems and
and for homogenization problems.

##  1 Introduction

In this paper we present an abstract result of implicit function theorem type
(see Section 2), which in the past has been applied in [5, 6, 7, 16, 20, 22,
23] to singularly perturbed nonlinear ODEs and PDEs and, in Part II [17], to
periodic homogenization of semilinear elliptic PDE systems. In the present
paper we apply it to describe periodic homogenization for systems of
quasilinear second-order ODEs in divergence form of the type

$a(x,x/\varepsilon,u(x),u^{\prime}(x))^{\prime}=f(x,x/\varepsilon,u(x),u^{\prime}(x))\mbox{
for }x\in[0,1]$ (1.1)

with one Dirichlet and one natural boundary condition

$u(0)=a(1,1/\varepsilon,u(1),u^{\prime}(1))=0$ (1.2)

as well as with other boundary conditions (see Section 4). Here
$\varepsilon>0$ is the small homogenization parameter, and we look for vector
valued solutions $u:[0,1]\to\mathbb{R}^{n}$. The coefficient functions
$a,f:[0,1]\times\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$
are supposed to be 1-periodic with respect to the second argument, i.e.

$a(x,y+1,u,u^{\prime})=a(x,y,u,u^{\prime})\mbox{ and
}f(x,y+1,u,u^{\prime})=f(x,y,u,u^{\prime})$ (1.3)

for all $x\in[0,1]$, $y\in\mathbb{R}$ and $u,u^{\prime}\in\mathbb{R}^{n}$.
Further, we suppose that the maps $a$ and $f$ are continuous and that their
first partial derivatives with respect to the third and fourth arguments exist
and are continuous, i.e.

$a,\partial_{u}a,\partial_{u^{\prime}}a,f,\partial_{u}f,\partial_{u^{\prime}}f\mbox{
are continuous on
}[0,1]\times\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}^{n}.$ (1.4)

And finally, we supose that the maps $a(x,y,u,\cdot)$ are strongly monotone
and Lipschitz continuous uniformly with respect to $x$, $y$ and bounded $u$,
i.e. that there exist constants $M\geq m>0$ such that for all $x\in[0,1]$,
$y\in\mathbb{R}$ and $u,u_{1}^{\prime},u_{2}^{\prime}\in\mathbb{R}^{n}$ with
$\|u\|\leq 1$ we have

$\left.\begin{array}[]{l}\big{(}a(x,y,u,u_{1}^{\prime})-a(x,y,u,u_{2}^{\prime})\big{)}\cdot(u_{1}^{\prime}-u_{2}^{\prime})\geq
m\|u_{1}^{\prime}-u_{2}^{\prime}\|^{2},\\\
\|a(x,y,u,u_{1}^{\prime})-a(x,y,u,u_{2}^{\prime})\|\|\leq
M\|u^{\prime}_{1}-u^{\prime}_{2}\|.\end{array}\right\\}$ (1.5)

Here and in what follows we denote by $v\cdot w$ the Euclidean scalar product
of vectors $v,w\in\mathbb{R}^{n}$, and $\|v\|:=\sqrt{v\cdot v}$ is the
Euclidean norm of the vector $v\in\mathbb{R}^{n}$.

From assumption (1.5) it follows that for all $x\in[0,1]$, $y\in\mathbb{R}$
and $u\in\mathbb{R}^{n}$ with $\|u\|\leq 1$ the maps $a(x,y,u,\cdot)$ are
bijective from $\mathbb{R}^{n}$ onto $\mathbb{R}^{n}$. We denote

$b(x,y,u,\cdot):=a(x,y,u,\cdot)^{-1},\mbox{ i.e.
}a(x,y,u,b(x,y,u,u^{\prime}))=b(x,y,u,a(x,y,u,u^{\prime}))=u^{\prime}$ (1.6)

and

$b_{0}(x,u,u^{\prime}):=\int_{0}^{1}b(x,y,u,u^{\prime})dy$ (1.7)

for all $x\in[0,1]$, $y\in\mathbb{R}$ and $u,u^{\prime}\in\mathbb{R}^{n}$ with
$\|u\|\leq 1$. Also the maps $b_{0}(x,u,\cdot)$ are strongly monotone and
Lipschitz continuous and, hence, bijective from $\mathbb{R}^{n}$ onto
$\mathbb{R}^{n}$, and we denote

$a_{0}(x,u,\cdot):=b_{0}(x,u,\cdot)^{-1},\mbox{ i.e.
}a_{0}(x,u,b_{0}(x,u,u^{\prime}))=b_{0}(x,u,a_{0}(x,u,u^{\prime}))=u^{\prime}$
(1.8)

and

$f_{0}(x,u,u^{\prime}):=\int_{0}^{1}f(x,y,u,b(x,y,u,u^{\prime}))dy$ (1.9)

for $x\in[0,1]$ and $u,u^{\prime}\in\mathbb{R}^{n}$ with $\|u\|\leq 1$, and
the boundary value problem

$a_{0}(x,u(x),u^{\prime}(x))^{\prime}=f_{0}(x,u(x),u^{\prime}(x))\mbox{ for
}x\in[0,1],\;u(0)=a_{0}(1,u(1),u^{\prime}(1))=0$ (1.10)

is the homogenized version of the boundary value problem (1.1)-(1.2).

A vector function $u\in C^{1}([0,1];\mathbb{R}^{n})$ is called weak solution
to (1.1)-(1.2) if it satisfies the Dirichlet boundary condition $u(0)=0$ and
the variational equation

$\left.\begin{array}[]{r}\displaystyle\int_{0}^{1}\Big{(}a(x,x/\varepsilon,u(x),u^{\prime}(x))\cdot\varphi^{\prime}(x)+f(x,x/\varepsilon,u(x),u^{\prime}(x))\cdot\varphi(x)\Big{)}dx=0\\\
\mbox{ for all }\varphi\in C^{1}([0,1];\mathbb{R}^{n})\mbox{ with
}\varphi(0)=0,\end{array}\right\\}$ (1.11)

and similar for the homogenized boundary value problem (1.10) and its
linearization (1.12). Weak solutions to (1.1)-(1.2) are not classical
solutions, i.e. they are not $C^{2}$-smooth, in general, because the maps
$a(\cdot,\cdot,u,u^{\prime})$ are not $C^{1}$-smooth, in general.

Now we formulate our result about existence and local uniqueness of weak
solutions $u=u_{\varepsilon}$ to (1.1)-(1.2) with $\varepsilon\approx 0$,
which are close to a given non-degenerate solution $u=u_{0}$ to (1.10), and
about the rate of convergence to zero for $\varepsilon\to 0$ of the
homogenization error $\|u_{\varepsilon}-u_{0}\|_{\infty}$. Here and in what
follows we denote by

$\|u\|_{\infty}:=\max_{x\in[0,1]}\|u(x)\|$

the maximum norm in the function space $C([0,1];\mathbb{R}^{n})$.

###### Theorem 1.1

Suppose (1.3)-(1.5), and let $u=u_{0}$ be a weak solution to (1.10) such that
$\|u_{0}\|_{\infty}<1$ and that the linearized boundary value problem

$\left.\begin{array}[]{l}\Big{(}\partial_{u}a_{0}(x,u_{0}(x),u_{0}^{\prime}(x))u(x)+\partial_{u}^{\prime}a_{0}(x,u_{0}(x),u_{0}^{\prime}(x))u^{\prime}(x)\Big{)}^{\prime}\\\
=\partial_{u}f_{0}(x,u_{0}(x),u_{0}^{\prime}(x))u(x)+\partial_{u^{\prime}}f_{0}(x,u_{0}(x),u_{0}^{\prime}(x))u^{\prime}(x)\mbox{
for }x\in[0,1],\\\
u(0)=\partial_{u}a_{0}(1,u_{0}(1),u_{0}^{\prime}(1))u(1)+\partial_{u^{\prime}}a_{0}(1,u_{0}(1),u_{0}^{\prime}(1))u^{\prime}(1)=0\end{array}\right\\}$
(1.12)

does not have weak solutions $u\not=0$. Then the following is true:

(i) There exist $\varepsilon_{0}>0$ and $\delta>0$ such that for all
$\varepsilon\in(0,\varepsilon_{0}]$ there exists exactly one weak solution
$u=u_{\varepsilon}$ to (1.1)-(1.2) with $\|u-u_{0}\|_{\infty}\leq\delta$.
Moreover,

$\|u_{\varepsilon}-u_{0}\|_{\infty}\to 0\mbox{ for }\varepsilon\to 0.$ (1.13)

(ii) If also

$\partial_{x}a,\partial_{x}f\mbox{ exist and are continuous on
}[0,1]\times\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}^{n},$ (1.14)

then

$\|u_{\varepsilon}-u_{0}\|_{\infty}=O(\varepsilon)\mbox{ for }\varepsilon\to
0.$ (1.15)

###### Remark 1.2

Our notion of weak solutions $u$ to (1.1)-(1.2) does not include the
requirement $u\in W^{1,2}((0,1);\mathbb{R}^{n})$, as usual, but the stronger
requirement $u\in C^{1}([0,1];\mathbb{R}^{n})$. We do this in order to avoid
to suppose growth restrictions for the reaction functions $f(x,y,u,\cdot)$. If
we would suppose that the functions $f(x,y,u,\cdot)$ have linear growth, then
any solution $u\in W^{1,2}((0,1);\mathbb{R}^{n})$ to (1.11) with $u(0)=0$
would be $C^{1}$-smooth and, hence, a weak solution to (1.1)-(1.2) in the
sense introduced above.

###### Remark 1.3

The homogenized version $a_{0}$ of the map $a$ depends on $a$ only (cf. (1.6)
and (1.8)), but the homogenized version $f_{0}$ of the map $f$ depends not
only on $f$, but also on $a$ (cf. (1.9)), i.e. the homogenization of $f$ is
”relative to $a$”. For linear problems this effect is well-known, cf. [2,
Remark 1.13.1], [28] and [29, formula (3.9)].

###### Remark 1.4

It is easy to verify that the map $b$, which is defined in (1.6), is
continuous and its first partial derivatives with respect to the third and
fourth arguments exist and are continuous, that the maps
$b(x,\cdot,u,u^{\prime})$ are 1-periodic, and that

$\displaystyle\big{(}b(x,y,u,u_{1}^{\prime})-b(x,y,u,u_{2}^{\prime})\big{)}\cdot(u_{1}^{\prime}-u_{2}^{\prime})\geq\frac{m}{M^{2}}\|u_{1}^{\prime}-u_{2}^{\prime}\|^{2},$
(1.16)
$\displaystyle\|b(x,y,u,u_{1}^{\prime})-b(x,y,u,u_{2}^{\prime})\|\leq\frac{1}{m}\|u^{\prime}_{1}-u^{\prime}_{2}\|$
(1.17)

for all $x\in[0,1]$, $y\in\mathbb{R}$ and
$u,u_{1}^{\prime},u_{2}^{\prime}\in\mathbb{R}^{n}$ with $\|u\|\leq 1$.
Similarly, also the map $a_{0}$, which is defined in (1.8), is continuous and
its first partial derivatives with respect to the third and fourth arguments
exist and are continuous, and

$\displaystyle(a_{0}(x,u,u^{\prime}_{1})-a_{0}(x,u,u^{\prime}_{2}))\cdot(u^{\prime}_{1}-u^{\prime}_{1})$
$\displaystyle\geq$
$\displaystyle\frac{m^{3}}{M^{2}}\|u^{\prime}_{1}-u^{\prime}_{2}\|^{2},$
$\displaystyle\|a_{0}(x,u,u^{\prime}_{1})-a_{0}(x,u,u^{\prime}_{2})\|$
$\displaystyle\leq$
$\displaystyle\frac{M^{2}}{m}\|u^{\prime}_{1}-u^{\prime}_{2}\|$

for all $x\in[0,1]$, $u,u^{\prime}_{1},u^{\prime}_{2}\in\mathbb{R}^{n}$ with
$\|u\|\leq 1$.

###### Remark 1.5

In many applications the maps $a(x,y,u,\cdot)$ and $f(x,y,u,\cdot)$ are
affine, i.e.

$a(x,y,u,u^{\prime})=A(x,y,u)u^{\prime}+\bar{a}(x,y,u)\mbox{ and
}f(x,y,u,u^{\prime})=F(x,y,u)u^{\prime}+\bar{f}(x,y,u)$

with $n\times n$-matrices $A(x,y,u)$ and $F(x,y,u)$ and vectors
$\bar{a}(x,y,u)$ and $\bar{f}(x,y,u)$. Then also the maps $a_{0}(x,u,\cdot)$
and $f_{0}(x,u,\cdot)$ are affine, i.e.

$a_{0}(x,u,u^{\prime})=A_{0}(x,u)u^{\prime}+\bar{a}_{0}(x,u)\mbox{ and
}f_{0}(x,u,u^{\prime})=F_{0}(x,u)u^{\prime}+\bar{f}_{0}(x,u)$

with

$\displaystyle A_{0}(x,u)$ $\displaystyle:=$
$\displaystyle\left(\int_{0}^{1}A(x,y,u)^{-1}dy\right)^{-1},$
$\displaystyle\bar{a}_{0}(x,u)$ $\displaystyle:=$
$\displaystyle\left(\int_{0}^{1}A(x,y,u)^{-1}dy\right)^{-1}\int_{0}^{1}A(x,y,u)^{-1}\bar{a}(x,y,u)dy,$
$\displaystyle F_{0}(x,u)$ $\displaystyle:=$
$\displaystyle\int_{0}^{1}F(x,y,u)A(x,y,u)^{-1}dy,$
$\displaystyle\bar{f}_{0}(x,u)$ $\displaystyle:=$
$\displaystyle\int_{0}^{1}\left(\bar{f}(x,y,u)-F(x,y,u)A(x,y,u)^{-1}\bar{a}(x,y,u)\right)dy.$

Hence, $A_{0}$ depends on $A$ only, but $\bar{a}_{0}$ depends on $\bar{a}$ and
$A$, $F_{0}$ depends on $F$ and $A$, and $\bar{f}_{0}$ depends on $\bar{f}$,
$A$, $\bar{a}$ and $F$.

###### Remark 1.6

The assumption of Theorem 1.1, that there do not exist nontrivial weak
solutions to (1.10), is rather implicit. But there exist simple explicit
sufficient contitions for it. For example, if not only the matrices
$\partial_{u^{\prime}}a_{0}(x,u_{0}(x),u^{\prime}_{0}(x))$ are positive
definit (this follows from assumption (1.5)), but also the matrices
$\partial_{u^{\prime}}f_{0}(x,u_{0}(x),u^{\prime}_{0}(x))$, and if the
corresponding definitness coefficients are sufficiently large in comparison
with the matrix norms of $\partial_{u}a_{0}(x,u_{0}(x),u^{\prime}_{0}(x))$ and
$\partial_{u}f_{0}(x,u_{0}(x),u^{\prime}_{0}(x))$, then there do not exist
nontrivial weak solutions to (1.10). In order to verify this one can use the
formulas

$\displaystyle\partial_{u^{\prime}}a_{0}(x,u,u^{\prime})=\left(\int_{0}^{1}\partial_{u^{\prime}}a(x,y,u,u^{\prime})^{-1}dy\right)^{-1},$
$\displaystyle\partial_{u}a_{0}(x,u,u^{\prime})=\left(\int_{0}^{1}\partial_{u^{\prime}}a(x,y,u,u^{\prime})^{-1}dy\right)^{-1}\int_{0}^{1}\partial_{u^{\prime}}a(x,y,u,u^{\prime})^{-1}\partial_{u}a(x,y,u,u^{\prime})dy,$
$\displaystyle\partial_{u^{\prime}}f_{0}(x,u,u^{\prime})=\int_{0}^{1}\partial_{u^{\prime}}f(x,y,u,u^{\prime})\partial_{u^{\prime}}a(x,y,u,u^{\prime})^{-1}dy,$
$\displaystyle\partial_{u}f_{0}(x,u,u^{\prime})=\int_{0}^{1}\left(\partial_{u}f(x,y,u,u^{\prime})-\partial_{u^{\prime}}f(x,y,u,u^{\prime})\partial_{u^{\prime}}a(x,y,u,u^{\prime})^{-1}\partial_{u}a(x,y,u,u^{\prime})\right)dy.$

###### Remark 1.7

The assertions of Theorem 1.1 remain true also in cases where the maps
$a(x,\cdot,u,u^{\prime})$ or $f(x,\cdot,u,u^{\prime})$ are allowed to be
discontinous, for example, if

$a(x,y,u,u^{\prime})=a_{1}(x,u,u^{\prime})a_{2}(y)\mbox{ and
}f(x,y,u,u^{\prime})=f_{1}(x,u,u^{\prime})f_{2}(y)$

with vector functions $a_{1},f_{1}\in
C^{1}([0,1]\times\mathbb{R}^{n}\times\mathbb{R}^{n};\mathbb{R}^{n})$ and
1-periodic functions $a_{2},f_{2}\in L^{\infty}(\mathbb{R})$, or, more
general, if the maps $(x,u,u^{\prime})\mapsto a(x,\cdot,u,u^{\prime})$ and
$(x,u,u^{\prime})\mapsto f(x,\cdot,u,u^{\prime})$ are continuous from
$[0,1]\times\mathbb{R}^{n}\times\mathbb{R}^{n}$ into $L^{\infty}(\mathbb{R})$
(cf. also Remark 3.3). For the case of linear scalar equations see [2,
Theorems 6.1 and 6.3] and [30, Theorem 1.2].

###### Remark 1.8

$L^{\infty}$-estimates of the homogenization error $u_{\varepsilon}-u_{0}$
exist, to the best of our knowledge, for linear homogenization problems only:
For scalar ODEs of the type
$\big{(}a(x/\varepsilon)u^{\prime}(x)\big{)}^{\prime}=f(x)$ (with a smooth
1-periodic function $a:\mathbb{R}\to\mathbb{R}$ and a smooth function
$f:[0,1]\to\mathbb{R}$) in [18, Section 1], for scalar ODEs with stratified
structure of the type
$\big{(}a(x,\rho(x)/\varepsilon)u^{\prime}(x)\big{)}^{\prime}=f(x)$ in [30,
Theorem 1.2]. For $L^{\infty}$ homogenization error estimates for scalar
linear elliptic PDEs of the type $\mbox{\rm div}\,a(x/\varepsilon)\nabla
u(x)=f(x)$ see, e.g. [2, Chapter 2.4] and [15] and for linear elliptic systems
[24, Theorem 7.5.1].

What concerns existence and local uniqueness for nonlinear homogenization
problems (without assumption of global uniqueness) we know only the result [4]
for scalar semilinear elliptic PDEs of the type $\mbox{\rm
div}\,a(x/\varepsilon)\nabla u(x)=f(x)g(u(x)),$ where the nonlinearity $g$ is
supposed to have a sufficiently small local Lipschitz constant (on an
appropriate bounded interval). Let us mention also [13, 14], where existence
and local uniqueness for a homogenization problem for the linear Poisson
equation with periodic nonlinear Robin boundary conditions is shown. There the
specific structure of the problem (no highly oscillating diffusion
coefficients) allows to apply the classical implicit function theorem.

###### Remark 1.9

Consider the quasilinear elliptic PDE $\mbox{\rm
div}\,a(x,x/\varepsilon,u(x),\nabla u(x))=f(x)$ with flux function
$a:\Omega\times\mathbb{R}^{d}\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}^{d}$
(with $\Omega\subseteq\mathbb{R}^{d}$) such that $a(x,y+e_{j},u,v)=a(x,y,u,v)$
for all $j=1,\ldots,d$
($e_{1}:=(1,0,\ldots,0,0),\ldots,e_{d}:=(0,0,\ldots,0,1)$ is the standard
basis in $\mathbb{R}^{d}$) and that $a(x,y,u,\cdot)$ is strongly monotone and
Lipschitz continuous uniformly with respect to $x$, $y$ and bounded $u$. The
usual formula for the homogenized flux function is (cf., e.g. [8, 11, 21, 27])

$a_{0}(x,u,v):=\int_{[0,1]^{d}}a(x,y,u,v+\nabla_{y}w(x,y,u,v))dy,$ (1.18)

where $w(x,\cdot,u,v)$ is the solution (which depends parametrically on $x$,
$u$ and $v$) of the cell problem
${div}_{y}\;a(x,y,u,v+\nabla_{y}w(x,y,u,v))=0$, $w(x,y+e_{j},u,v)=w(x,y,u,v)$
and $\int_{[0,1]^{d}}w(x,y,u,v)dy=0$. In space dimension one, i.e. $d=1$, this
looks as follows:

$a_{0}(x,u,v):=\int_{0}^{1}a(x,y,u,v+\partial_{y}w(x,y,u,v))dy$ (1.19)

and

$\left.\begin{array}[]{l}\displaystyle\frac{d}{dy}a(x,y,u,v+\partial_{y}w(x,y,u,v))=0,\\\
w(x,y+1,u,v)=w(x,y,u,v),\;\displaystyle\int_{0}^{1}w(x,y,u,v)dy=0.\end{array}\right\\}$
(1.20)

From (1.20) it follows that $a(x,y,u,v+\partial_{y}w(x,y,u,v))$ is constant
with respect to $y$. Therefore (1.19) yields that
$a(x,y,u,v+\partial_{y}w(x,y,u,v))=a_{0}(x,u,v)$ and, hence,
$b(x,y,u,a_{0}(x,u,v))=v+\partial_{y}w(x,y,u,v),$ i.e.

$v=\int_{0}^{1}(v+\partial_{y}w(x,y,u,v))dy=\int_{0}^{1}b(x,y,u,a_{0}(x,u,v))dy=b_{0}(x,u,a_{0}(x,u,v)).$

On the other hand, the solution to (1.20) with $v=b_{0}(x,u,\bar{v})$ (with
arbitrary $\bar{v}\in\mathbb{R}$) is

$w(x,y,u,b_{0}(x,u,\bar{v}))=\int_{0}^{y}(b(x,z,u,\bar{v})-b_{0}(x,u,\bar{v}))dz-\int_{0}^{1}\int_{0}^{z_{1}}(b(x,z_{2},u,\bar{v})-b_{0}(x,u,\bar{v}))dz_{2}dz_{1}.$

Therefore

$\displaystyle a_{0}(x,u,b_{0}(x,u,\bar{v}))$ $\displaystyle=$
$\displaystyle\int_{0}^{1}a(x,y,u,b_{0}(x,u,\bar{v})-\partial_{y}w(x,y,u,b_{0}(x,u,\bar{v})))dy$
$\displaystyle=$
$\displaystyle\int_{0}^{1}a(x,y,u,b(x,y,u,\bar{v})dy=\bar{v}.$

It follows that $a_{0}(x,u,\cdot)=b_{0}(x,u,\cdot)^{-1}$. In other words: Our
definition (1.8) of the homogenized flux function $a_{0}$ is the same as the
usual one for PDEs, i.e. for (1.18), considered in the case $d=1$. In the
linear case this has been shown in [2, Remark 2.3] and [10, Proposition 6.16].
In [2, Remark 5.9] the formulas (1.8) and (1.19) are called dual formulas
(there the linear case is considered, but with multidimensional space variable
$x$).

Our paper is organized as follows: In Section 2 we consider abstract nonlinear
parameter depending equations of the type

${\cal F}_{\varepsilon}(w)=0.$ (1.21)

Here $\varepsilon>0$ is the parameter. We prove a result on existence and
local uniqueness of a family of solutions $w=w_{\varepsilon}\approx w_{0}$ to
(1.21) with $\varepsilon\approx 0$, where $w_{0}$ is an approximate solution
to (1.21), i.e. an element with ${\cal F}_{\varepsilon}(w_{0})\to 0$ for
$\varepsilon\to 0$, and we estimate the norm of the error
$w_{\varepsilon}-w_{0}$ by the norm of the discrepancy ${\cal
F}_{\varepsilon}(w_{0})$. This type of generalized implicit function theorems
has been successfully applied to singularly perturbed ODEs and PDEs in [5, 6,
7, 16, 20, 22, 23]). Contrary to the classical implicit function theorem it is
not supposed that the linearized operators ${\cal
F}^{\prime}_{\varepsilon}(u)$ converge for $\varepsilon\to 0$ in the uniform
operator norm. And, indeed, in the applications to singularly perturbed
problems as well as to periodic homogenization problems they do not converge
for $\varepsilon\to 0$ in the uniform operator norm (cf. Remark 3.6 below).
Hence, the present paper is a first step (on the ODE level) to create a common
approach to existence, local uniqueness and error estimates for singularly
perturbed problems and for homogenization problems. In Part II [17] we apply
this approach to periodic homogenization for semilinear elliptic PDE systems.

In Section 3 we prove Theorem 1.1 by means of the results of Section 2. For
that reason we transform the boundary value problem (1.1)-(1.2) into the
system (3.3)-(3.4) of integral equations, and for that system of integral
equations we introduce an abstract setting of the type (1.21). For that
abstract setting we have to verify the key assumptions (2.1) and (2.4) of
Theorem 2.1, and we do this in the Subsections 3.1 and 3.2, respectively.

Finally, in Section 4 we show that Theorem 1.1 remains true also for
inhomogeneous natural boundary conditions, but not, in general, for
inhomogeneous Neumann boundary conditions. Remark that the difficulties with
inhomogeneous Neumann boundary conditions are well-known already for scalar
linear problems (see, e.g. [2, Remark 1.2.10 and Section 1.7.1] and [26]).
Further, we show how to prove that the assertions of Theorem 1.1 are true also
for two Dirichlet boundary conditions.

##  2 An abstract result of implicit function theorem type

In this section we formulate and prove Theorem 2.1 below.

###### Theorem 2.1

Let be given a Banach space $W$ with norm $\|\cdot\|_{W}$, an open set
$W_{0}\subseteq W$, an element $w_{0}\in W_{0}$ and a family of $C^{1}$-maps
${\cal F}_{\varepsilon}:W_{0}\to W$ with $\varepsilon>0$ as family parameter.
Suppose that

$\|{\cal F}_{\varepsilon}(w_{0})\|_{W}\to 0\mbox{ for }\varepsilon\to 0.$
(2.1)

Further, suppose that there exists $\varepsilon_{0}>0$ such that

$\displaystyle{\cal F}^{\prime}_{\varepsilon}(w_{0})\mbox{ is Fredholm of
index zero from $W$ into $W$ for all $\varepsilon\in(0,\varepsilon_{0}]$},$
(2.4) $\displaystyle\inf\\{\|{\cal
F}^{\prime}_{\varepsilon}(w_{0})w\|_{W}:\;\varepsilon\in(0,\varepsilon_{0}],\;\|w\|_{W}=1\\}=:\alpha>0,$
$\displaystyle\sup_{\|w\|_{W}\leq 1}\|({\cal
F}^{\prime}_{\varepsilon}(w_{0}+w_{1})-{\cal
F}_{\varepsilon}^{\prime}(w_{0}))w\|_{W}\to 0\mbox{ for
}\varepsilon+\|w_{1}\|_{W}\to 0.$

Then there exist $\varepsilon_{1}\in(0,\varepsilon_{0}]$ and $\delta>0$ such
that for all $\varepsilon\in(0,\varepsilon_{1}]$ there exists exactly one
$w=w_{\varepsilon}\in W_{0}$ with ${\cal F}_{\varepsilon}(w)=0$ and
$\|w-w_{0}\|_{W}\leq\delta$. Moreover,

$\|w_{\varepsilon}-w_{0}\|_{W}<\frac{2}{\alpha}\|{\cal
F}_{\varepsilon}(w_{0})\|_{W}.$ (2.5)

Proof  Assumptions (2.4) and (2.4) imply, that for all
$\varepsilon\in(0,\varepsilon_{0})$ the operator ${\cal
F}_{\varepsilon}^{\prime}(w_{0})$ is an isomorphism from $W$ onto $W$ and

$\left\|{\cal
F}_{\varepsilon}^{\prime}(w_{0})^{-1}w\right\|_{W}\leq\frac{1}{\alpha}\|w\|_{W}\mbox{
for all }w\in W.$ (2.6)

Hence, the map ${\cal G}_{\varepsilon}:W_{0}\to W$,

${\cal G}_{\varepsilon}(w):=w-{\cal F}_{\varepsilon}^{\prime}(w_{0})^{-1}{\cal
F}_{\varepsilon}(w)$

is well-defined. Obviously, $w$ is a fixed point of ${\cal G}_{\varepsilon}$
if and only if ${\cal F}_{\varepsilon}(w)=0$.

For $r>0$ denote $\mathbb{B}_{r}:=\\{w\in W:\;\|w-w_{0}\|_{W}\leq r\\}.$ We
are going to show that for sufficiently small $\varepsilon>0$ and $r>0$ the
map ${\cal G}_{\varepsilon}$ is strictly contractive from the closed ball
$\mathbb{B}_{r}$ into itself.

In order to verify the strict contractivity of ${\cal G}_{\varepsilon}$ we
take $\varepsilon\in(0,\varepsilon_{0}]$ and $v,w\in W_{0}$ and estimate as
follows:

$\displaystyle\|{\cal G}_{\varepsilon}(v)-{\cal
G}_{\varepsilon}(w)\|_{W}=\left\|v-w-{\cal F}_{\varepsilon}(w_{0})^{-1}({\cal
F}_{\varepsilon}(v)-{\cal F}_{\varepsilon}(w))\right\|_{W}$
$\displaystyle=\left\|{\cal
F}_{\varepsilon}^{\prime}(w_{0})^{-1}\int_{0}^{1}\left({\cal
F}_{\varepsilon}^{\prime}(w_{0})-{\cal
F}^{\prime}_{\varepsilon}(sv+(1-s)w)\right)ds(v-w)\right\|_{W}$
$\displaystyle\leq\frac{1}{\alpha}\int_{0}^{1}\|\left({\cal
F}_{\varepsilon}^{\prime}(w_{0})-{\cal
F}^{\prime}_{\varepsilon}(sv+(1-s)w)\right)(v-w)\|_{W}ds.$

Here we used (2.6). Because of assumption (2.4) there exist
$\varepsilon_{1}\in(0,\varepsilon_{0}]$ and $r_{0}>0$ such that
$\mathbb{B}_{r_{0}}\subset W_{0}$ and $\|\left({\cal
F}_{\varepsilon}^{\prime}(w_{0})-{\cal
F}^{\prime}_{\varepsilon}(sv+(1-s)w)\right)(v-w)\|_{W}\leq\frac{\alpha}{2}\|v-w\|_{W}$
for all $\varepsilon\in(0,\varepsilon_{1}]$, $s\in[0,1]$ and
$v,w\in\mathbb{B}_{r_{0}}$. Hence,

$\|{\cal G}_{\varepsilon}(v)-{\cal
G}_{\varepsilon}(w)\|_{W}\leq\frac{1}{2}\|v-w\|_{W}\mbox{ for all
}\varepsilon\in(0,\varepsilon_{1}]\mbox{ and }v,w\in\mathbb{B}_{r_{0}}.$ (2.7)

Now, let us show that ${\cal G}_{\varepsilon}$ maps $\mathbb{B}_{r_{0}}$ into
$\mathbb{B}_{r_{0}}$ for all sufficiently small $\varepsilon>0$. Take
$\varepsilon\in(0,\varepsilon_{1}]$ and $w\in\mathbb{B}_{r_{0}}$. Then (2.6)
and (2.7) imply

$\displaystyle\left\|{\cal
G}_{\varepsilon}(w)-w_{0}\right\|_{W}\leq\left\|{\cal
G}_{\varepsilon}(w)-{\cal G}_{\varepsilon}(w_{0})\right\|_{W}+\left\|{\cal
G}_{\varepsilon}(w_{0})-w_{0}\right\|_{W}$
$\displaystyle\leq\frac{1}{2}\left\|w-w_{0}\right\|_{W}+\left\|{\cal
F}_{\varepsilon}^{\prime}(w_{0})^{-1}{\cal
F}_{\varepsilon}(w_{0})\right\|_{W}\leq\frac{r_{0}}{2}+\frac{1}{\alpha}\left\|{\cal
F}_{\varepsilon}(w_{0})\right\|_{W}.$

But assumption (2.1) yields that, if $\varepsilon_{1}$ is taken sufficiently
small, for all $\varepsilon\in(0,\varepsilon_{1}]$ we have $\|{\cal
F}_{\varepsilon}(w_{0})\|_{W}\leq\alpha r_{0}/2$. Hence, for those
$\varepsilon$ we get $\left\|{\cal G}_{\varepsilon}(w)-w_{0}\right\|_{W}\leq
r_{0}$.

Therefore, Banach’s fixed point principle yields the following: For all
$\varepsilon\in(0,\varepsilon_{1}]$ there exists exactly one
$w=w_{\varepsilon}\in\mathbb{B}_{r_{0}}$ with ${\cal F}_{\varepsilon}(w)=0$.

Finally, let us prove (2.5). We take $\varepsilon\in(0,\varepsilon_{1}]$ and
estimate as above:

$\|w_{\varepsilon}-w_{0}\|_{W}\leq\|{\cal
G}_{\varepsilon}(w_{\varepsilon})-{\cal G}_{\varepsilon}(w_{0})\|_{W}+\|{\cal
G}_{\varepsilon}(w)-w_{0}\|_{W}\leq\frac{1}{2}\|w_{\varepsilon}-w_{0}\|_{W}+\frac{1}{\alpha}\|{\cal
F}_{\varepsilon}(w_{0})\|_{W}.$

Hence, (2.5) is true.

###### Remark 2.2

In [5, 6, 7, 16, 20, 22, 23]) slightly more general versions of Theorem 2.1
are used, i.e. those with ${\cal F}_{\varepsilon}$ mapping one Banach space
into another one, both with $\varepsilon$-depending norms. Moreover, there the
approximate solutions are allowed to be $\varepsilon$-depending, i.e. to be a
family of approximate solutions. Hence, these versions of Theorem 2.1 seem to
be appropriate for applications to homogenization problems with approximate
solutions defined by using correctors of first or higher-order (see, e.g. [1,
9, 12]).

For another result of the type of Theorem 2.1 and its applications to
semilinear elliptic PDE systems with numerically determined approximate
solutions see [3, Theorem 2.1].

##  3 Proof of Theorem 1.1

In this section we will prove Theorem 1.1 by means of Theorem 2.1. Hence, all
assumptions of Theorem 1.1 (i.e. (1.3)-(1.5), existence of the weak solution
$u=u_{0}$ to (1.10), non-existence of weak solutions $u\not=0$ to (1.12)) will
be supposed to be satisfied (without mentioning their use). At places, where
we use the additional assumption (1.14) of Theorem 1.1(ii), we will mention
this.

In order to transform the problem of weak solutions $u\approx u_{0}$ to the
boundary value problem (1.1)-(1.2) into the problem of solutions $w\approx
w_{0}$ to an appropriate operator equation ${\cal F}_{\varepsilon}(w)=0$, we
use the notation

$\mathbb{B}:=\\{u\in C([0,1];\mathbb{R}^{n}):\;\|u\|_{\infty}<1\\}$ (3.1)

and Lemmas 3.1, 3.2 and 3.4 below.

###### Lemma 3.1

For all $\varepsilon>0$ the following is true:

(i) If $u\in\mathbb{B}$ is a weak solution to (1.1)-(1.2) and if $v\in
C([0,1];\mathbb{R}^{n})$ is defined by

$v(x):=a(x,x/\varepsilon,u(x),u^{\prime}(x))\mbox{ for }x\in[0,1],$ (3.2)

then

$\displaystyle u(x)=\int_{0}^{x}b(y,y/\varepsilon,u(y),v(y))dy\mbox{ for
}x\in[0,1],$ (3.3) $\displaystyle
v(x)=-\int_{x}^{1}f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v(y)))dy\mbox{
for }x\in[0,1].$ (3.4)

(ii) If $(u,v)\in\mathbb{B}\times C([0,1];\mathbb{R}^{n})$ is a solution to
(3.3)-(3.4), then $u$ is a weak solution to (1.1)-(1.2).

Proof  (i) Let $u\in\mathbb{B}$ be a weak solution to (1.1)-(1.2). Take an
arbitrary test function $\varphi\in C^{1}([0,1];\mathbb{R}^{n})$ with
$\varphi(0)=0$. Then (1.11) yields that

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\left(a(x,x/\varepsilon,u(x),u^{\prime}(x))\cdot\varphi^{\prime}(x)+f(x,x/\varepsilon,u(x),u^{\prime}(x))\cdot\int_{0}^{x}\varphi^{\prime}(y)dy\right)dx.$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\left(a(x,x/\varepsilon,u(x),u^{\prime}(x))+\int_{x}^{1}f(y,y/\varepsilon,u(y),u^{\prime}(y))dy\right)\cdot\varphi^{\prime}(x)dx.$

Therefore
$a(x,x/\varepsilon,u(x),u^{\prime}(x))+\int_{x}^{1}f(y,y/\varepsilon,u(y),u^{\prime}(y))dy$
is constant with respect to $x$. In particular, the function $x\mapsto
a(x,x/\varepsilon,u(x),u^{\prime}(x))$ is $C^{1}$-smooth, and

$a(x,x/\varepsilon,u(x),u^{\prime}(x))^{\prime}=f(x,x/\varepsilon,u(x),u^{\prime}(x)).$
(3.5)

If $v\in C([0,1];\mathbb{R}^{n})$ is defined by (3.2), then

$u^{\prime}(x)=b(x,x/\varepsilon,u(x),v(x)).$ (3.6)

Inserting (3.2) and (3.6) into (3.5) we get (3.4). Moreover, the boundary
condition $u(0)=0$ and (3.6) yield (3.3).

(ii) Let $(u,v)\in\mathbb{B}\times C([0,1];\mathbb{R}^{n})$ be a solution to
(3.3)-(3.4). From (3.3) follows $u(0)=0$ and (3.6), and from (3.6) follows
$v(x)=a(x,x/\varepsilon,u(x),u^{\prime}(x))$. Therefore (3.4) implies
$v(1)=a(1,1/\varepsilon,u(1),u^{\prime}(1))=0$ and
$v^{\prime}(x)=a(x,x/\varepsilon,u(x),u^{\prime}(x))^{\prime}=f(x,x/\varepsilon,u(x),b(x,x/\varepsilon,u(x),v(x))).$
If we multiply this scalarly by an arbitrary test function $\varphi\in
C^{1}([0,1];\mathbb{R}^{n})$ with $\varphi(0)=0$, integrate with respect to
$x$ and use the boundary condition $v(1)=0$, then we get (1.11).

The following lemma is the only tool from classical homogenization theory
which we are going to use. For related results see, e.g. [18, Lemma 1.1], [24,
Proposition 2.2.2], [29, Lemma 3.1]. In Lemma 3.2 below we use the following
notation for maps $g\in C([0,1]\times\mathbb{R};\mathbb{R}^{n})$:

$\displaystyle\omega_{g}(\varepsilon)$ $\displaystyle:=$
$\displaystyle\sup\\{\|g(x_{1},y)-g(x_{2},y)\|:\;x_{1},x_{2}\in[0,1],\;y\in\mathbb{R},\;|x_{1}-x_{2}|\leq\varepsilon\\}\mbox{
for }\varepsilon>0,$ $\displaystyle\|g\|_{*}$ $\displaystyle:=$
$\displaystyle\sup\\{\|g(x,y)\|:\;(x,y)\in[0,1]\times\mathbb{R}\\}.$

###### Lemma 3.2

Let be given $g\in C([0,1]\times\mathbb{R};\mathbb{R}^{n})$ such that
$g(x,y+1)=g(x,y)$ for all $x\in[0,1]$ and $y\in\mathbb{R}$. Then for all
$x\in[0,1]$ and $\varepsilon>0$ we have

$\left\|\int_{0}^{x}\left(g(y,y/\varepsilon)-\int_{0}^{1}g(y,z)dz\right)dy\right\|\leq
2\left(\omega_{g}(\varepsilon)+\varepsilon\|g\|_{*}\right)$ (3.7)

and, if the partial derivative $\partial_{x}g$ exists and is continuous,

$\left\|\int_{0}^{x}\left(g(y,y/\varepsilon)-\int_{0}^{1}g(y,z)dz\right)dy\right\|\leq
2\varepsilon\left(\|g\|_{*}+\|\partial_{x}g\|_{*}\right).$ (3.8)

Proof  Define

$h(x,y):=g(x,y)-\int_{0}^{1}g(x,z)dz.$

Then $h(x,y+1)=h(x,y)$ and $\int_{y}^{y+1}h(x,z)dz=0$ and
$\omega_{h}(\varepsilon)\leq 2\omega_{g}(\varepsilon)$. Therefore for
$x\in[0,1]$ and $\varepsilon>0$ we have

$\displaystyle\int_{0}^{x}\left(g(y,y/\varepsilon)-\int_{0}^{1}g(y,z)dz\right)dy=\int_{0}^{x}h(y,y/\varepsilon)dy=\varepsilon\int_{0}^{x/\varepsilon}h(\varepsilon
y,y)dy$
$\displaystyle=\varepsilon\left(\sum_{j=1}^{[x/\varepsilon]}\int_{j-1}^{j}h(\varepsilon
y,y)dy+\int_{[x/\varepsilon]}^{x/\varepsilon}h(\varepsilon y,y)dy\right)$
$\displaystyle=\varepsilon\left(\sum_{j=1}^{[x/\varepsilon]}\int_{j-1}^{j}\left(h(\varepsilon
y,y)dy-h(\varepsilon
j,y)\right)dy+\int_{[x/\varepsilon]}^{x/\varepsilon}h(\varepsilon
y,y)dy\right),$

where $[x/\varepsilon]$ is the integer part of $x/\varepsilon$, i.e. the
largest integer which is not larger than $x/\varepsilon$. In particular,
$\varepsilon[x/\varepsilon]\leq x\leq 1$.

For $y\in[j-1,j]$ we have that $0\leq\varepsilon(j-y)\leq\varepsilon$ and,
hence, that $|h(\varepsilon y,y)-h(\varepsilon j,y)|\leq w_{h}(\varepsilon)$.
Therefore

$\left\|\int_{0}^{x}\left(g(y,y/\varepsilon)-\int_{0}^{1}g(y,z)dz\right)dy\right\|\leq\varepsilon\left([x/\varepsilon]w_{h}(\varepsilon)+\|h\|_{*}\right)\leq
2\left(w_{g}(\varepsilon)+\varepsilon\|g\|_{*}\right),$

i.e. (3.7) is proved.

If $\partial_{x}g$ exists and is continuous, then
$g(x_{1},y)-g(x_{2},y)=(x_{1}-x_{2})\int_{0}^{1}\partial_{x}g(sx_{1}+(1-s)x_{2},y)ds$,
i.e. $\omega_{g}(\varepsilon)\leq\varepsilon\|\partial_{x}g\|_{*}$. Hence, in
that case (3.7) implies (3.8).

###### Remark 3.3

If $g(x,y)=g_{1}(x)g_{2}(y)$ with $g_{1}\in C([0,1];\mathbb{R}^{n})$ and
$g_{2}\in L^{\infty}(\mathbb{R})$ or, more general, if the map $x\mapsto
g(x,\cdot)$ is continuous from $[0,1]$ into $L^{\infty}(\mathbb{R})$, then the
assertions of Lemma 3.2 remain true.

Similarly to Lemma 3.1 we get, that the function $u_{0}$, which is by
assumtion of Theorem 1.1 a weak solution to (1.12), and the function $v_{0}\in
C([0,1];\mathbb{R}^{n})$, which is defined by

$v_{0}(x):=a_{0}(x,u_{0}(x),u_{0}^{\prime}(x))\mbox{ for }x\in[0,1],$ (3.9)

satisfy

$\left.\begin{array}[]{l}\displaystyle
u_{0}(x)=\int_{0}^{x}b_{0}(y,u_{0}(y),v_{0}(y))dy,\\\ \displaystyle
v_{0}(x)=-\int_{x}^{1}f_{0}(y,u_{0}(y),b_{0}(y,u_{0}(y),v_{0}(y)))dy\end{array}\right\\}\mbox{
for }x\in[0,1].$ (3.10)

###### Lemma 3.4

For all $\gamma>0$ there exists $\delta>0$ such that for all $\varepsilon>0$,
$u\in\mathbb{B}$ and $v\in C([0,1];\mathbb{R}^{n})$ with (3.4) and
$\varepsilon+\|u-u_{0}\|_{\infty}\leq\delta$ we have
$\|v-v_{0}\|_{\infty}\leq\gamma$.

Proof  Take arbitrary $\varepsilon>0$, $u\in\mathbb{B}$ and $v\in
C([0,1];\mathbb{R}^{n})$. Because of (3.4) and (3.10) we have

$\displaystyle v(x)-v_{0}(x)$ $\displaystyle=$
$\displaystyle\int_{x}^{1}\left(f_{0}(y,u_{0}(y),b_{0}(y,u_{0}(y),v_{0}(y)))-f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v(y))\right)dy$
$\displaystyle=$
$\displaystyle\alpha_{\varepsilon}(x)+\alpha_{\varepsilon,u}(x)+\alpha_{\varepsilon,u,v}(x)$

with

$\displaystyle\alpha_{\varepsilon}(x)$
$\displaystyle:=\int_{x}^{1}\left(f_{0}(y,u(y),b(y,y/\varepsilon,u(y),v(y)))-f(y,y/\varepsilon,u_{0}(y),b(y,y/\varepsilon,u_{0}(y),v_{0}(y))\right)dy,$
$\displaystyle\alpha_{\varepsilon,u}(x)$
$\displaystyle:=\int_{x}^{1}\left(f(y,y/\varepsilon,u_{0}(y),b(y,y/\varepsilon,u_{0}(y),v_{0}(y)))-f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v_{0}(y))\right)dy,$
$\displaystyle\alpha_{\varepsilon,u,v}(x)$
$\displaystyle:=\int_{x}^{1}\left(f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v_{0}(y)))-f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v(y))\right)dy.$

From (1.9) it follows that $\alpha_{\varepsilon}(x)$ equals to

$\int_{x}^{1}\left(\int_{0}^{1}f(y,z,u_{0}(y),b(y,z,u_{0}(y),v_{0}(y)))dz-f(y,y/\varepsilon,u_{0}(y),b(y,y/\varepsilon,u_{0}(y),v_{0}(y))\right)dy,$

and Lemma 3.2 with $g(x,y)=f(x,y,u_{0}(x),b(x,y,u_{0}(x),v_{0}(x)))$ yields
that $\|\alpha_{\varepsilon}\|_{\infty}\to 0$ for $\varepsilon\to 0$. Further,
we have

$\displaystyle\alpha_{\varepsilon,u}(x)=\int_{x}^{1}\int_{0}^{1}\Big{(}\partial_{u}f(y,y/\varepsilon,u_{s}(y),b(y,y/\varepsilon,u_{s}(y),v_{0}(y)))$
$\displaystyle\;\;\;\;\;+\partial_{u^{\prime}}f(y,y/\varepsilon,u_{s}(y),b(y,z,u_{s}(y),v_{0}(y))\partial_{u}b(y,y/\varepsilon,u_{s}(y),v_{0}(y))\Big{)}ds\,(u_{0}(y)-u(y))dy$

with $u_{s}(y):=su_{0}(y)+(1-s)u(y)$. Hence,
$\|\alpha_{\varepsilon,u}\|_{\infty}\leq\mbox{const}\|u-u_{0}\|_{\infty},$
where the constant does not depend on $\varepsilon$ and $u$. Finally, we have

$\alpha_{\varepsilon,u,v}(x)=\int_{x}^{1}\int_{0}^{1}\partial_{u^{\prime}}f(y,y/\varepsilon,u(y),b_{s}(y))ds\,(b(y,y/\varepsilon,u(y),v_{0}(y))-b(y,y/\varepsilon,u(y),v(y)))dy$

with
$b_{s}(y):=sb(y,y/\varepsilon,u(y),v_{0}(y))+(1-s)b(y,y/\varepsilon,u(y),v(y))$.
Hence, because of (1.17) there exists a constant $c>0$, which does not depend
on $\varepsilon$, $x$, $u$ and $v$, such that

$\|\alpha_{\varepsilon,u,v}\|\leq c\int_{x}^{1}\|v(y)-v_{0}(y)\|dy.$

It follows that
$\|v(x)-v_{0}(x)\|\leq\|\alpha_{\varepsilon}\|_{\infty}+\|\alpha_{\varepsilon,u}\|_{\infty}+c\int_{x}^{1}\|v(y))-v_{0}(y))\|dy,$
and Gronwall’s inequality yields that
$\|v(x)-v_{0}(x)\|\leq\left(\|\alpha_{\varepsilon}\|_{\infty}+\|\alpha_{\varepsilon,u}\|_{\infty}\right)\exp\left(c(1-x)\right),$
i.e.

$\|v-v_{0}\|_{\infty}\to 0\mbox{ for }\varepsilon+\|u-u_{0}\|_{\infty}\to 0.$

Now we are going to apply Theorem 2.1 in order to solve the boundary value
problem (1.1)-(1.2) with $\varepsilon\approx 0$ and
$\|u-u_{0}\|_{\infty}\approx 0$. We introduce the setting of Theorem 2.1 as
follows:

$\displaystyle
W:=C([0,1];\mathbb{R}^{n})^{2},\;\|(u,v)\|_{W}:=\|u\|_{\infty}+\|v\|_{\infty},\;W_{0}:=\mathbb{B}\times
C([0,1];\mathbb{R}^{n}),\;w_{0}:=(u_{0},v_{0}),$ $\displaystyle{\cal
F}_{\varepsilon}(u,v)=({\cal U}_{\varepsilon}(u,v),{\cal
V}_{\varepsilon}(u,v))$

with

$\displaystyle[{\cal
U}_{\varepsilon}(u,v)](x):=u(x)-\int_{0}^{x}b(y,y/\varepsilon,u(y),v(y))dy$
$\displaystyle[{\cal
V}_{\varepsilon}(u,v)](x):=v(x)+\int_{x}^{1}f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v(y)))dy.$

Here $\mathbb{B}$ is the open ball in $C([0,1];\mathbb{R}^{n})$, defined in
(3.1), $u_{0}$ is the solution to the linearized boundary value problen
(1.12), which is given by assumption of Theorem 1.1, and $v_{0}$ is defined in
(3.9)

Because of Lemmas 3.1, 3.2 and 3.4 we have the following: If $u\in\mathbb{B}$
is a weak solution to (1.1)-(1.2) with $\varepsilon\approx 0$ and
$\|u-u_{0}\|_{\infty}\approx 0$, then there exists $v\in
C([0,1];\mathbb{R}^{n})$ with $\|v-v_{0}\|_{\infty}\approx 0$ such that ${\cal
F}_{\varepsilon}(u,v)=0$. And if $(u,v)\in W_{0}$ satisfies ${\cal
F}_{\varepsilon}(u,v)=0$ with $\varepsilon\approx 0$ and
$\|u-u_{0}\|_{\infty}+\|v-v_{0}\|_{\infty}\approx 0$, then $u$ is a weak
solution to (1.1)-(1.2). Moreover, if all the assumptions of Theorem 2.1, i.e.
(2.1)-(2.4), are satisfied in the setting introduced above, then Theorem 2.1
yields the assertions of Theorem 1.1(i), in particular (2.5) yields (1.13).
If, moreover,

$\mbox{assumption (\ref{diffass1}) implies that }\|{\cal
F}_{\varepsilon}(u_{0},v_{0})\|_{W}=O(\varepsilon)\mbox{ for }\varepsilon\to
0$ (3.11)

in the setting introduced above, then the assertion (2.5) of Theorem 2.1
yields also assertion (1.15) of Theorem 1.1.

Hence, it remains to verify the assumptions (2.1)-(2.4) of Theorem 2.1 and the
assertion (3.11) in the setting introduced above.

###  3.1 Verification of (2.1) and (3.11)

Because of (1.7), (1.9) and (3.10) we have

$[{\cal
U}_{\varepsilon}(u_{0},v_{0})](x)=\int_{0}^{x}\left(\int_{0}^{1}b(y,z,u_{0}(y),v_{0}(y))dz-b(y,y/\varepsilon,u_{0}(y),v_{0}(y))\right)dy$

and

$\displaystyle[{\cal V}_{\varepsilon}(u_{0},v_{0})](x)$
$\displaystyle=\int_{x}^{1}\Big{(}f(y,y/\varepsilon,u_{0}(y),b(y,y/\varepsilon,u_{0}(y),v_{0}(y)))-\int_{0}^{1}f(y,z,u_{0}(y),b(y,z,u_{0}(y),v_{0}(y)))dzdz\Big{)}dy.$

Hence, Lemma 3.2 with

$g(x,y)=b(x,x/\varepsilon,u_{0}(x),v_{0}(x))\mbox{ and
}g(x,y)=f(x,x/\varepsilon,u_{0}(x),b(x,x/\varepsilon,u_{0}(x),v_{0}(x))),$

respectively, yields for $\varepsilon\to 0$ that

$\|{\cal U}_{\varepsilon}(u_{0},v_{0})\|_{\infty}+\|{\cal
V}_{\varepsilon}(u_{0},v_{0})\|_{\infty}=\left\\{\begin{array}[]{l}o(1),\\\
O(\varepsilon),\mbox{ if (\ref{diffass1}) is satified.}\end{array}\right.$

###  3.2 Verification of (2.4)-(2.4)

We have

$[{\cal F}^{\prime}_{\varepsilon}(u,v)](\bar{u},\bar{v})=(\partial_{u}{\cal
U}_{\varepsilon}(u,v)\bar{u}+\partial_{v}{\cal
U}_{\varepsilon}(u,v)\bar{v},\partial_{u}{\cal
V}_{\varepsilon}(u,v)\bar{u}+\partial_{v}{\cal V}_{\varepsilon}(u,v)\bar{v})$

with

$\displaystyle[\partial_{u}{\cal
U}_{\varepsilon}(u,v)\bar{u}](x)=\bar{u}(x)-\int_{0}^{x}\partial_{u}b(y,y/\varepsilon,u(y),v(y))\bar{u}(y)dy,$
$\displaystyle[\partial_{v}{\cal
U}_{\varepsilon}(u,v)\bar{v}](x)=-\int_{0}^{x}\partial_{u^{\prime}}b(y,y/\varepsilon,u(y),v(y))\bar{v}(y)dy,$
$\displaystyle[\partial_{u}{\cal
V}_{\varepsilon}(u,v)\bar{u}](x)=\int_{x}^{1}\Big{(}\partial_{u}f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v(y))$
$\displaystyle\;\;\;\;+\partial_{u^{\prime}}f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v(y))\partial_{u}b(y,y/\varepsilon,u(y),v(y))\Big{)}\bar{u}(y)dy,$
$\displaystyle[\partial_{v}{\cal V}_{\varepsilon}(u,v)\bar{v}](x)=\bar{v}(x)$
$\displaystyle\;\;\;\;+\int_{x}^{1}\partial_{u^{\prime}}f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v(y))\partial_{u^{\prime}}b(y,y/\varepsilon,u(y),v(y))\bar{v}(y)dy.$

In order to verify assumption (2.4) of Theorem 2.1 we calculate as follows:

$\displaystyle[(\partial_{u}{\cal
U}_{\varepsilon}(u_{0},v_{0})-\partial_{u}{\cal
U}_{\varepsilon}(u_{1},v_{1}))\bar{u}](x)$
$\displaystyle=\int_{0}^{x}\partial_{u}b(y,y/\varepsilon,u_{1}(y),v_{1}(y))-\partial_{u}b(y,y/\varepsilon,u_{0}(y),v_{0}(y))\bar{u}(y)dy.$

But $\partial_{u}b$ is uniformly continuous on
$\\{(x,y,u,v)\in[0,1]\times\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}^{n}:\;\|u\|\leq
1,\;\|v\|\leq R\\}$ for all $R>0$. Hence,

$\|\partial_{u}b(y,y/\varepsilon,u_{0}(y),v_{0}(y))-\partial_{u}b(y,y/\varepsilon,u_{1}(y),v_{1}(y))\bar{u}(y)\|\to
0\mbox{ for }\|u_{0}-u_{1}\|_{\infty}+\|v_{0}-v_{1}\|_{\infty}\to 0$

uniformly with respect to $\varepsilon>0$, $y\in[0,1]$ and
$\|\bar{u}\|_{\infty}\leq 1$. Similarly one can estimate the terms
$(\partial_{v}{\cal U}_{\varepsilon}(u_{0},v_{0})-\partial_{u}{\cal
U}_{\varepsilon}(u_{1},v_{1}))\bar{v}$, $(\partial_{u}{\cal
V}_{\varepsilon}(u_{0},v_{0})-\partial_{u}{\cal
V}_{\varepsilon}(u_{1},v_{1}))\bar{u}$ and $(\partial_{v}{\cal
V}_{\varepsilon}(u_{0},v_{0})-\partial_{u}{\cal
V}_{\varepsilon}(u_{1},v_{1}))\bar{v}$.

Further, for any $(u,v)\in{\cal W}$ we have that ${\cal
F}^{\prime}_{\varepsilon}(u,v)-I$ ($I$ is the identity in $W$) is a linear
bounded operator from $W$ into $C^{1}([0,1];\mathbb{R}^{n})^{2}$, where
$C^{1}([0,1];\mathbb{R}^{n})$ is equipped with its usual norm
$\|u\|_{\infty}+\|u^{\prime}\|_{\infty}$. Hence, the Arcela-Ascoli Theorem
yields that for any $(u,v)\in W$ the operator ${\cal
F}^{\prime}_{\varepsilon}(u,v)-I$ is compact from $W$ into $W$, and,
therefore, the operator ${\cal F}^{\prime}_{\varepsilon}(u,v)$ is Fredholm of
index zero from $W$ into $W$, i.e. (2.4) is satisfied.

Now, let us verify (2.4).

Suppose that (2.4) is not true, i.e. that it is not true that there exists
$\varepsilon_{0}>0$ such that $\inf\\{\|{\cal
F}_{\varepsilon}^{\prime}(w_{0})w\|_{W}:\;\varepsilon\in(0,\varepsilon_{0}],w\in
W,\|w\|_{W}=1\\}>0$. Then there exist sequences
$\varepsilon_{1},\varepsilon_{2},\ldots>0$ and $u_{1},u_{2},\ldots\in
C([0,1];\mathbb{R}^{n})$ and $v_{1},v_{2},\ldots\in C([0,1];\mathbb{R}^{n})$
such that

$\left.\begin{array}[]{r}\displaystyle\lim_{k\to\infty}\varepsilon_{k}=0,\\\
\displaystyle\lim_{k\to\infty}\|\partial_{u}{\cal
U}_{\varepsilon_{k}}(u_{0},v_{0})u_{k}+\partial_{v}{\cal
U}_{\varepsilon_{k}}(u_{0},v_{0})v_{k}\|_{\infty}=0,\\\
\displaystyle\lim_{k\to\infty}\|\partial_{u}{\cal
V}_{\varepsilon_{k}}(u_{0},v_{0})u_{k}+\partial_{v}{\cal
V}_{\varepsilon_{k}}(u_{0},v_{0})v_{k}\|_{\infty}=0,\end{array}\right\\}$
(3.12)

but

$\|u_{k}\|_{\infty}+\|v_{k}\|_{\infty}=1\mbox{ for all }k\in\mathbb{N}.$
(3.13)

Denote

$\displaystyle\bar{u}_{k}(x):=\int_{0}^{x}\Big{(}\partial_{u}b(y,y/\varepsilon_{n},u_{0}(y),v_{0}(y))u_{k}(y)+\partial_{u^{\prime}}b(y,y/\varepsilon_{n},u_{0}(y),v_{0}(y))v_{k}(y)\Big{)}dy,$
(3.14)
$\displaystyle\bar{v}_{k}(x):=\int_{x}^{1}\Big{(}\big{(}\partial_{u}f(y,y/\varepsilon_{n},u_{0}(y),b(y,y/\varepsilon_{n},u_{0}(y),v_{0}(y))$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;+\partial_{u^{\prime}}f(y,y/\varepsilon_{n},u_{0}(y),b(y,y/\varepsilon_{n},u_{0}(y),v_{0}(y))\partial_{u}b(y,y/\varepsilon_{n},u_{0}(y),v_{0}(y))\big{)}u_{k}(y)$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;+\partial_{u^{\prime}}f(y,y/\varepsilon_{n},u_{0}(y),b(y,y/\varepsilon_{n},u_{0}(y),v_{0}(y))\partial_{u^{\prime}}b(y,y/\varepsilon_{n},u_{0}(y),v_{0}(y))v_{k}(y)\Big{)}dy.$

Then

$\sup_{k\in\mathbb{N}}\big{(}\|\bar{u}^{\prime}_{k}\|_{\infty}+\|\bar{v}^{\prime}_{k}\|_{\infty}\big{)}<\infty.$
(3.15)

Hence, because of the Arzela-Ascoli Theorem without loss of generality we may
assume that there exist $\bar{u}_{0},\bar{v}_{0}\in C([0,1];\mathbb{R}^{n})$
such that

$\lim_{k\to\infty}\big{(}\|\bar{u}_{k}-\bar{u}_{0}\|_{\infty}+\|\bar{v}_{k}-\bar{v}_{0}\|_{\infty}\big{)}=0.$
(3.16)

But we have that

$\displaystyle u_{k}=\bar{u}_{k}+\partial_{u}{\cal
U}_{\varepsilon_{k}}(u_{0},v_{0})u_{k}+\partial_{v}{\cal
U}_{\varepsilon_{k}}(u_{0},v_{0})v_{k},$ $\displaystyle
v_{k}=\bar{v}_{n}-\partial_{u}{\cal
V}_{\varepsilon_{k}}(u_{0},v_{0})u_{k}+\partial_{v}{\cal
V}_{\varepsilon_{k}}(u_{0},v_{0})v_{k}.$

Hence, (3.12) yields that
$\lim_{k\to\infty}\big{(}\|u_{k}-\bar{u}_{0}\|_{\infty}+\|v_{k}-\bar{v}_{0}\|_{\infty}\big{)}=0,$
and (3.13) implies that

$\|\bar{u}_{0}\|_{\infty}+\|\bar{v}_{0}\|_{\infty}=1.$ (3.17)

We are going to show that (3.16) and (3.17) lead to a contradiction. Because
of (1.7) we have for all $x\in[0,1]$ that

$\displaystyle\int_{0}^{x}\left(\partial_{u}b_{0}(y,u_{0}(y),v_{0}(y))\bar{u}_{0}(y)+\partial_{u^{\prime}}b_{0}(y,u_{0}(y),v_{0}(y))\bar{v}_{0}(y)\right)dy$
$\displaystyle=\int_{0}^{x}\left(\partial_{u}\int_{0}^{1}b(y,z,u_{0}(y),v_{0}(y))dz\;\bar{u}_{0}(y)+\partial_{u^{\prime}}\int_{0}^{1}b(y,z,u_{0}(y),v_{0}(y))dz\;\bar{v}_{0}(y)\right)dy$
$\displaystyle=\int_{0}^{x}\int_{0}^{1}\left(\partial_{u}b(y,z,u_{0}(y),v_{0}(y))dz\bar{u}_{0}(y)+\partial_{u^{\prime}}b(y,z,u_{0}(y),v_{0}(y))dz\bar{v}_{0}(y)\right)dzdy,$

and because of Lemma 3.2 with
$g(x,y)=\partial_{u}b(x,y,u_{0}(x),v_{0}(x))\bar{u}_{0}(x)+\partial_{u^{\prime}}b(x,y,u_{0}(x),v_{0}(x))\bar{v}_{0}(x)$
this equals to

$\lim_{k\to\infty}\int_{0}^{x}\left(\partial_{u}b(y,y/\varepsilon_{k},u_{0}(y),v_{0}(y))dz\bar{u}_{0}(y)+\partial_{u^{\prime}}b(y,y/\varepsilon_{k},u_{0}(y),v_{0}(y))dz\bar{v}_{0}(y)\right)dy,$

and because of (3.16) this is

$\lim_{k\to\infty}\int_{0}^{x}\left(\partial_{u}b(y,y/\varepsilon_{k},u_{0}(y),v_{0}(y))\bar{u}_{k}(y)+\partial_{u^{\prime}}b(y,y/\varepsilon_{k},u_{0}(y),v_{0}(y))\bar{v}_{k}(y)\right)dy,$

and, finally, because of (3.14) and (3.16) this equals to
$\lim_{k\to\infty}\bar{u}_{k}=\bar{u}_{0}.$ We end up with

$\bar{u}_{0}(x)=\int_{0}^{x}\Big{(}\partial_{u}b_{0}(y,u_{0}(y),v_{0}(y))\bar{u}_{0}(y)+\partial_{u^{\prime}}b_{0}(y,u_{0}(y),v_{0}(y))\bar{v}_{0}(y)\Big{)}dy.$
(3.18)

Similarly one shows that

$\displaystyle\bar{v}_{0}(x)=-\int_{0}^{x}\Big{(}\big{(}\partial_{u}f_{0}(y,u_{0}(y),b_{0}(y,u_{0}(y),v_{0}(y))$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;+\partial_{u^{\prime}}f_{0}(y,u_{0}(y),b_{0}(y,u_{0}(y),v_{0}(y))\partial_{u}b_{0}(y,u_{0}(y),v_{0}(y))\big{)}\bar{u}_{0}(y)$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;+\partial_{u^{\prime}}f_{0}(y,u_{0}(y),b_{0}(y,u_{0}(y),v_{0}(y))\partial_{u^{\prime}}b_{0}(y,u_{0}(y),v_{0}(y))\bar{v}_{0}(y)\Big{)}dy$

and similarly to Lemma 3.1 it follows that $\bar{u}_{0}$ is a weak solution to
the linearized boundary value problem (1.12). Therefore, the assumption of
Theorem 1.1, that (1.12) does not have nontrivial solution, implies that
$u_{0}=0$. Then (3.18) implies that
$\partial_{u^{\prime}}b_{0}(x,u_{0}(x),v_{0}(x))\bar{v}_{0}(x)=0$ for all
$x\in[0,1]$. But (1.16) yields that

$\partial_{u^{\prime}}b_{0}(x,u_{0}(x),v_{0}(x))v\cdot
v\geq\frac{m}{M^{2}}\|v\|^{2}$ (3.19)

for all $x\in[0,1]$ and $v\in\mathbb{R}^{n}$. We get that $\bar{v}_{0}=0$,
what contradicts to (3.17).

###### Remark 3.5

In the proof above we used the well-known fact (cf., e.g. [19]) that the
operations linearization and homogenization commute.

###### Remark 3.6

It is easy to verify that the linear operators ${\cal
F}^{\prime}_{\varepsilon}(u,v)$ do not converge for $\varepsilon\to 0$ in the
uniform operator norm in ${\cal L}(W)$, in general. For example,
$\int_{0}^{x}\partial_{u}b(y,y/\varepsilon,u(y),v(y))\bar{u}(y)dy$ (with fixed
$u,v\in C([0,1];\mathbb{R}^{n})$) does not converge for $\varepsilon\to 0$
uniformly with respect to $x\in[0,1]$ and $\bar{u}\in C([0,1];\mathbb{R}^{n})$
with $\|\bar{u}\|_{\infty}\leq 1$, in general. But in the proof above a
subsequence of the sequence
$\int_{0}^{x}\partial_{u}b(y,y/\varepsilon_{k},u_{0}(y),v_{0}(y))\bar{u}_{k}(y)dy$
converges for $k\to\infty$ uniformly with respect to $x\in[0,1]$, and this is
because of (3.15).

##  4 Other boundary conditions

###  4.1 Inhomogeneous boundary conditions

If the homogeneous natural boundary condition in (1.2) is replaced by an
corresponding inhomogeneous one, i.e.

$u(0)=0,\;a(1,1/\varepsilon,u(1),u^{\prime}(1))=u^{1},$ (4.1)

then the weak formulation of the boundary value problem (1.1),(4.1) is as
follows: Find $u\in W^{1,2}((0,1);\mathbb{R}^{n})$ such that $u(0)=0$ and

$\begin{array}[]{r}\displaystyle\int_{0}^{1}\Big{(}a(x,x/\varepsilon,u(x),u^{\prime}(x))\cdot\varphi^{\prime}(x)+f(x,x/\varepsilon,u(x),u^{\prime}(x))\cdot\varphi(x)\Big{)}dx=u^{1}\cdot\varphi(1)\\\
\mbox{ for all }\varphi\in C^{1}([0,1];\mathbb{R}^{n})\mbox{ with
}\varphi(0)=0.\end{array}$

The system (3.3)-(3.4) of integral equations has to be replaced by

$\displaystyle u(x)=\int_{0}^{x}b(y,y/\varepsilon,u(y),v(y))dy,$
$\displaystyle
v(x)=u^{1}-\int_{x}^{1}f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v(y)))dy,$

and the results of Theorem 1.1 remain unchanged.

But if the inhomogeneous natural boundary condition in $x=1$ in (4.1) is
replaced by an inhomogeneous Neumann condition, i.e.
$u(0)=u^{0},\;u^{\prime}(1)=u^{1}\not=0,$ then, in general, there do not exist
solution families $u=u_{\varepsilon}$ which converge pointwise for
$\varepsilon\to 0$. For example, the solution to the scalar linear boundary
value problem

$\left(\frac{u^{\prime}(x)}{2+\sin(2\pi
x/\varepsilon)}\right)^{\prime}=1\mbox{ for
}x\in[0,1],\;u(0)=0,u^{\prime}(1)=u^{1}$ (4.2)

is

$u_{\varepsilon}(x)=\int_{0}^{x}(2+\sin(2\pi
y/\varepsilon))\left(\frac{u^{1}}{2+\sin(2\pi/\varepsilon)}+y-1\right)dy.$

Hence,

$u_{\varepsilon}(1)=-1+\frac{2u^{1}}{2+\sin(2\pi/\varepsilon)}+O(\varepsilon)\mbox{
for }\varepsilon\to 0,$

i.e. $u_{\varepsilon}(1)$ does not converge for $\varepsilon\to 0$ if
$u^{1}\not=0$.

###### Remark 4.1

Consider the scalar linear boundary value problem (4.2) with $u^{1}=0$. Then

$u_{\varepsilon}(1)=\int_{0}^{1}(2+\sin(2\pi
y/\varepsilon))(y-1)dy=-1-\frac{\varepsilon}{2\pi}+O(\varepsilon^{2})\mbox{
for }\varepsilon\to 0,$

and the solution to the corresponding homogenized problem
$\frac{1}{2}u^{\prime\prime}(x)=1$ for $x\in[0,1]$, $u(0)=u^{\prime}(1)=0$ is
$u_{0}(x)=x(x-2)$. Hence,

$u_{\varepsilon}(1)-u_{0}(1)=-\frac{\varepsilon}{2\pi}+O(\varepsilon^{2})\mbox{
for }\varepsilon\to 0,$

i.e. the asymptotic error estimate (1.15) of Theorem 1.1 is sharp.

###  4.2 Two Dirichlet boundary conditions

If we consider (1.1) with two homogeneous Dirichlet boundary conditions, i.e.
$u(0)=u(1)=0,$ then the weak formulation is as follows: Find $u\in
W^{1,2}((0,1);\mathbb{R}^{n})$ such that $u(0)=u(1)=0$ and

$\left.\begin{array}[]{r}\displaystyle\int_{0}^{1}\Big{(}a(x,x/\varepsilon,u(x),u^{\prime}(x))\cdot\varphi^{\prime}(x)+f(x,x/\varepsilon,u(x),u^{\prime}(x))\cdot\varphi(x)\Big{)}dx=0\\\
\mbox{ for all }\varphi\in C^{1}([0,1];\mathbb{R}^{n})\mbox{ with
}\varphi(0)=\varphi(1)=0,\end{array}\right\\}$ (4.3)

and the system (3.3)-(3.4) has to be changed to the system

$\displaystyle u(x)$ $\displaystyle=$
$\displaystyle\int_{0}^{x}b(y,y/\varepsilon,u(y),v(y))dy,$ $\displaystyle
v(x)$ $\displaystyle=$ $\displaystyle
w-\int_{x}^{1}f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v(y)))dy,$
$\displaystyle 0$ $\displaystyle=$
$\displaystyle\int_{0}^{1}b(x,x/\varepsilon,u(x),v(x))dx$

with unknowns $(u,v,w)\in C([0,1];\mathbb{R}^{n})\times
C([0,1];\mathbb{R}^{n})\times\mathbb{R}^{n}$. This system of integral
equations can be treated by Theorem 2.1 in the following setting:

$W:=C([0,1];\mathbb{R}^{n})^{2}\times\mathbb{R}^{n},\;\|(u,v,w)\|_{W}:=\|u\|_{\infty}+\|v\|_{\infty}+\|w\|,\;W_{0}:=\mathbb{B}\times
C([0,1];\mathbb{R}^{n})\times\mathbb{R}^{n}.$

The approximate solution $w_{0}$ of Theorem 2.1 is the triple
$(u_{0},v_{0},w_{0})\in C([0,1];\mathbb{R}^{n})^{2}\times\mathbb{R}^{n},$
where $u_{0}$ is the given weak solution of the homogenized problem, i.e.
$u_{0}\in W^{1,2}((0,1);\mathbb{R}^{n})$ with $\|u_{0}\|_{\infty}<1$ and
$u_{0}(0)=u_{0}(1)=0$ and

$\begin{array}[]{r}\displaystyle\int_{0}^{1}\Big{(}a_{0}(x,u_{0}(x),u_{0}^{\prime}(x))\cdot\varphi^{\prime}(x)+f_{0}(x,u_{0}(x),u_{0}^{\prime}(x))\cdot\varphi(x)\Big{)}dx=0\\\
\mbox{ for all }\varphi\in C^{1}([0,1];\mathbb{R}^{n})\mbox{ with
}\varphi(0)=\varphi(1)=0,\end{array}$

and $v_{0}(x):=a_{0}(x,u_{0}(x),u_{0}^{\prime}(x))$, $w_{0}:=v_{0}(1)$. The
homogenized vector functions $a_{0}$ and $f_{0}$ are defined as in (1.8) and
(1.9). Finally, the maps ${\cal F}_{\varepsilon}:W\to W$ are defined by

${\cal F}_{\varepsilon}(u,v,w):=({\cal U}_{\varepsilon}(u,v),{\cal
V}_{\varepsilon}(u,v,w),{\cal W}_{\varepsilon}(u,v))$

with

$\displaystyle[{\cal
U}_{\varepsilon}(u,v)](x):=u(x)-\int_{0}^{x}b(y,y/\varepsilon,u(y),v(y))dy,$
$\displaystyle[{\cal
V}_{\varepsilon}(u,v,w)](x):=v(x)-w+\int_{x}^{1}f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v(y)))dy,$
$\displaystyle{\cal
W}_{\varepsilon}(u,v):=\int_{0}^{1}b(x,x/\varepsilon,u(x),v(x))dx.$

As in Subsection 3.1 it follows that

$\displaystyle[{\cal
U}_{\varepsilon}(u_{0},v_{0})](x)=\int_{0}^{x}\Big{(}b_{0}(y,u_{0}(y),v_{0}(y))-b(y,y/\varepsilon,u_{0}(y),v_{0}(y))\Big{)}dy,$
$\displaystyle[{\cal V}_{\varepsilon}(u_{0},v_{0},w_{0})](x)$
$\displaystyle=\int_{x}^{1}\Big{(}f(y,y/\varepsilon,u(y),b(y,y/\varepsilon,u(y),v(y)))-f_{0}(y,u_{0}(y),b_{0}(y,u_{0}(y),v_{0}(y)))\Big{)}dy,$

and, hence, for $\varepsilon\to 0$ follows that

$\|{\cal U}_{\varepsilon}(u_{0},v_{0})\|_{\infty}+\|{\cal
V}_{\varepsilon}(u_{0},v_{0},w_{0})\|_{\infty}=\left\\{\begin{array}[]{l}o(1),\\\
O(\varepsilon),\mbox{ if (\ref{diffass1}) is satified.}\end{array}\right.$

Further, we have

$\displaystyle{\cal W}_{\varepsilon}(u_{0},v_{0})$ $\displaystyle=$
$\displaystyle\int_{0}^{1}b(x,x/\varepsilon,u_{0}(x),v_{0}(x))dx$
$\displaystyle=$
$\displaystyle\int_{0}^{1}(b(x,x/\varepsilon,u_{0}(x),v_{0}(x))-b_{0}(x,u_{0}(x),v_{0}(x))dx$

and, hence, for $\varepsilon\to 0$ follows

$\|{\cal
W}_{\varepsilon}(u_{0},v_{0})\|_{\infty}=\left\\{\begin{array}[]{l}o(1),\\\
O(\varepsilon),\mbox{ if (\ref{diffass1}) is satified.}\end{array}\right.$

Therefore, the assumptions (2.1) and (3.11) are satisfied in the setting
introduced above.

Finally, let us verify the assumption (2.4) of Theorem 2.1 in the setting
introduced above. Suppose that (2.4) is not true. Then there exist sequences
$\varepsilon_{1},\varepsilon_{2},\ldots>0$ and $u_{1},u_{2},\ldots\in
C([0,1];\mathbb{R}^{n})$ and $v_{1},v_{2},\ldots\in C([0,1];\mathbb{R}^{n})$
and $w_{1},w_{2},\ldots\in\mathbb{R}^{n}$ such that $\varepsilon_{n}\to 0$ for
$n\to\infty$ and

$\displaystyle\lim_{k\to\infty}\|\partial_{u}{\cal
U}_{\varepsilon_{k}}(u_{0},v_{0})u_{k}+\partial_{v}{\cal
U}_{\varepsilon_{k}}(u_{0},v_{0})v_{k}\|_{\infty}$ $\displaystyle=$
$\displaystyle 0,$ (4.4) $\displaystyle\lim_{k\to\infty}\|\partial_{u}{\cal
V}_{\varepsilon_{k}}(u_{0},v_{0})u_{k}+\partial_{v}{\cal
V}_{\varepsilon_{k}}(u_{0},v_{0})v_{k}+\partial_{w}{\cal
V}_{\varepsilon_{k}}(u_{0},v_{0})w_{k}\|_{\infty}$ $\displaystyle=$
$\displaystyle 0,$ (4.5) $\displaystyle\lim_{k\to\infty}\|\partial_{u}{\cal
W}_{\varepsilon_{k}}(u_{0},v_{0})u_{k}+\partial_{v}{\cal
W}_{\varepsilon_{k}}(u_{0},v_{0})v_{k}\|$ $\displaystyle=$ $\displaystyle 0,$

but

$\|u_{k}\|_{\infty}+\|v_{k}\|_{\infty}+\|w_{k}\|=1\mbox{ for all
}k\in\mathbb{N}.$ (4.6)

As in Subsection 3.2 one can show that without loss of generality we can
assume that there exist $\bar{u}_{0},\bar{v}_{0}\in C([0,1];\mathbb{R}^{n})$
and $\bar{w}_{0}\in\mathbb{R}^{n}$ such that

$\lim_{k\to\infty}\left(\|u_{k}-\bar{u}_{0}\|_{\infty}+\|v_{k}-\bar{v}_{0}\|_{\infty}+\|w_{k}-\bar{w}_{0}\|\right)=0$

and that $\bar{u}_{0}$ is a solution to the linearization (in $u_{0}$) of
(4.3), i.e. that $\bar{u}_{0}=0$. Hence, from (4.4) follows that for all
$x\in[0,1]$ we have

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\lim_{k\to\infty}[\partial_{v}{\cal
U}_{\varepsilon}(u_{0},v_{0})v_{k}](x)=\lim_{k\to\infty}\int_{0}^{x}\partial_{u^{\prime}}b(y,y/\varepsilon,u_{0}(y),v_{0}(y))v_{k}(y)dy$
$\displaystyle=$
$\displaystyle\int_{0}^{x}\partial_{u^{\prime}}b_{0}(y,u_{0}(y),v_{0}(y))\bar{v}_{0}(y)dy.$

It follows that
$\partial_{u^{\prime}}b_{0}(x,u_{0}(x),v_{0}(x))\bar{v}_{0}(x)=0$ for all
$x\in[0,1]$ and, hence, (3.19) yields that $\bar{v}_{0}=0$. Therefore, (4.5)
implies that

$0=\lim_{k\to\infty}\|\partial_{w}{\cal
V}(u_{0},v_{0},w_{0})\bar{w}_{k}\|_{\infty}=\lim_{k\to\infty}\|w_{k}\|=\|\bar{w}_{0}\|,$

and we get a contradiction to (4.6).

## References

  * [1] G. Allaire, M. Briane, M. Vanninathan, A comparison between two-scale asymptotic expansions and Bloch wave expansions for the homogenization of periodic structures. SeMA 73 (2016), 237-259.
  * [2] A. Bensoussan, J.L. Lions, G. Papanicolaou, Asymptotic Analysis for Periodic Structures. Studies in Mathematics and its Applications vol. 3, North-Holland, 1978.
  * [3] M. Breden, R. Castelli, Existence and instability of steady states for a triangular cross-diffusion system: a computer-assisted proof. J. Differ. Equations 264 (2018), 6418–6458.
  * [4] R. Bunoiu, R. Precup, Localization and multiplicity in the homogenization of nonlinear problems. Adv. Nonlinear Anal. 9 (2020), 292–304.
  * [5] V.F. Butuzov, N.N. Nefedov, O.E. Omel’chenko, L. Recke, Time-periodic boundary layer solutions to singularly perturbed parabolic problems. J. Differ. Equations 262 (2017), 4823–4862.
  * [6] V.F. Butuzov, N.N. Nefedov, O.E. Omel’chenko, L. Recke, Boundary layer solutions to singularly perturbed quasilinear systems. Discrete Cont. Dyn. Syst., Series B 27 (2022), 4255–4283.
  * [7] V.F. Butuzov, N.N. Nefedov, O.E. Omel’chenko, L. Recke, K.R. Schneider, An implicit function theorem and applications to nonsmooth boundary layers. In: Patterns of Dynamics, ed. by P. Gurevich, J. Hell, B. Sandstede, A. Scheel, Springer Proc. in Mathematics & Statistics vol. 205, Springer, 2017, 111–127.
  * [8] G. Cardone, S. E. Pastukhova, V. V. Zhikov, Some estimates for non-linear homogenization. Rend. Accad. Naz. Sci. XL, Mem. Mat. Appl. 24 (2005) 101-110 (2005).
  * [9] K. D. Cherednichenko, V. P. Smyshlyaev, On full two-scale expansion of the solutions of nonlinear periodic rapidly oscillating problems and higher-order homogenised variational problems. Arch. Rational Mech. Anal.174 (2004), 385–442.
  * [10] D. Cioranescu, P. Donato, An Introduction to Homogenization. Oxford Lecture Series in Mathematics and its Applications vol. 17, Oxford University Press, 1999.
  * [11] N. Fusco, G. Moscariello, On the homogenization of quasilinear divergence structure operators. Ann. Mat. Pura Appl., IV. Ser. 146 (1987) 1-13. ,
  * [12] Tebib Hawa, Chacha Djamal Ahmed, Third-order corrections in periodic homogenization for elliptic problem. Mediterr. J. Math. 18 (2021), Paper No. 135.
  * [13] M. Lanza de Cristoforis, P. Musolino, Two-parameter homogenization for a nonlinear periodic Robin problem for a Poisson equation: a functional analytic approach. Rev. Mat. Complut. 31 (2018), 63–110.
  * [14] M. Lanza de Cristoforis, P. Musolino, Asymptotic behaviour of the energy integral of a two-parameter homogenization problem with nonlinear periodic Robin boundary conditions. Proc. Edinb. Math. Soc. II. Ser. 62 (2019), 985–1016.
  * [15] Wen Ming He, Jun Zhi Cui, Error estimate of the homogenization solution for elliptic problems with small periodic coefficients on $L^{\infty}(\Omega)$. Science China Mathematics 53 (2010), 1231–1252.
  * [16] N.N. Nefedov, A.O. Orlov, L. Recke, K.R. Schneider, Nonsmooth regular perturbations of singularly perturbed problems. J. Differ. Equations 375 (2023), 206–236.
  * [17] N.N. Nefedov, L. Recke, A common approach to singular perturbation and homogenization II: Semilinear elliptic systems. arXiv:2309.14108.
  * [18] S. Neukamm, An introduction to the qualitative and quantitative theory of homogenization. Interdisciplinary Information Sciences 22 (2016),147–186.
  * [19] S. Neukamm, M. Schäffner, Quantitative homogenization in nonlinear elasticity for small loads. Arch. Ration. Mech. Anal. 230 (2018), 343–396.
  * [20] O.E. Omel’chenko, L. Recke, Existence, local uniqueness and asymptotic approximation of spike solutions to singularly perturbed elliptic problems. Hiroshima Math. J. 45 (2015), 35–89.
  * [21] S. E. Pastukhova, Operator estimates in nonlinear problems of reiterated homogenization. Proc. Steklov Inst. Math. 261 (2008), 214-228.
  * [22] L. Recke, Use of very weak approximate boundary layer solutions to spatially nonsmooth singularly perturbed problems. J. Math. Anal. Appl. 506 (2022), Article ID 125552.
  * [23] L. Recke, O.E. Omel’chenko, Boundary layer solutions to problems with infinite dimensional singular and regular perturbations. J. Differ. Equations 245 (2008), 3806–3822.
  * [24] Zongwei Shen, Periodic Homogenization of Elliptic Systems. Operator Theory: Advances and Applications vol. 269, Birkhäuser, 2018.
  * [25] Zongwei Shen, Jinping Zhuge, Convergence rates in periodic homogenization of systems of elasticity. Proc. Am. Math. Soc. 145 (2017), 1187–1202.
  * [26] T. Suslina, Homogenization of the Neumann problem for elliptic systems with periodic coefficients. SIAM J. Math. Anal. 45 (2013), 3453–3498.
  * [27] Li Wang, Qiang Xu, Peihao Zhao, Quantitative estimates on periodic homogenization of nonliear elliptic operators. arXiv:1807.10865v1 (2018).
  * [28] Qiang Xu, Convergence rates for general elliptic homogenization problems in Lipschitz domains. SIAM J. Math. Anal. 48 (2016), 3742–3788.
  * [29] Shixin Xu, Xingye Yue, Changrong Zhang, Homogenization: in mathematics or physics? Discrete Contin. Dyn. Syst. Series S 9 (2016), 1575–1590.
  * [30] Yao Xu, Weisheng Niu, Homogenization of elliptic systems with stratified structure revisited. Commun. Partial Differ. Equations 45 (2020), 655–689.


# Crystalline phases at finite winding densities in a quantum link ladder

Paolo Stornati<EMAIL_ADDRESS>ICFO-Institut de Ciencies Fotoniques,
The Barcelona Institute of Science and Technology, Av. Carl Friedrich Gauss 3,
08860 Castelldefels (Barcelona), Spain  Philipp Krah  TU Berlin, Institute of
Mathematics, Straße des 17. Juni 136, 10623 Berlin, Germany  Karl Jansen
Deutsches Elektronen-Synchrotron DESY,Platanenallee 6, 15738 Zeuthen, Germany
Debasish Banerjee  Theory Division, Saha Institute of Nuclear Physics, 1/AF
Bidhan Nagar, Kolkata 700064, India  Homi Bhabha National Institute, Training
School Complex, Anushaktinagar, Mumbai 400094,India

###### Abstract

Condensed matter physics of gauge theories coupled to fermions can exhibit a
rich phase structure, but are nevertheless very difficult to study in Monte
Carlo simulations when they are afflicted by a sign problem. As an alternate
approach, we use tensor network methods to explore the finite density physics
of Abelian gauge theories without dynamical matter. As a concrete example, we
consider the $U(1)$ gauge invariant quantum link ladder with
spin-$\frac{1}{2}$ gauge fields in an external electric field which cause the
winding electric fluxes to condense in the ground state. We demonstrate how
the electric flux tubes arrange themselves in the bulk giving rise to
crystalline patterns, whose period can be controlled by tuning the external
field. We propose observables to detect the transitions in ground state
properties not only in numerical experiments, but also in future cold-atom
realizations. A systematic procedure for reaching the thermodynamic limit, as
well as extending the studies from ladders to extended geometries is outlined.

## Introduction.–

Finite chemical potentials are expected to give rise to novel phases and
correlations otherwise absent in the ground state of quantum field theories or
quantum many-body systems. Two physically relevant examples are Quantum
Chromodynamics (QCD) and the Hubbard model. Markov Chain Monte Carlo (MCMC)
methods to solve QCD regulated on the lattice can explain properties of
hadrons, such as their masses, binding energies, and scattering cross-
sections. At finite baryon densities, $\mu_{B}$, relevant for e.g., the
description of the interior of neutron stars or the very early universe, the
MCMC methods suffer from the infamous sign problem. The Hubbard model, on the
other hand, is a pedagogical system to describe a variety of phases of
strongly correlated electrons. At finite doping, it is expected to host high-
temperature superconducting phases and provide a model for many physically
interesting materials. Once again, the regime of non-zero doping is difficult
to investigate numerically using Monte Carlo methods due to the sign problem.

Finite density physics of scalar and fermionic theories in various space-time
dimensions have been extensively investigated [1, 2, 3, 4, 5, 6, 7]. We extend
such studies which dealt with point particles to pure gauge theories without
dynamical matter fields containing loop operators. The simplest scenario is an
$U(1)$ Abelian lattice gauge theory in a finite volume and in 2+1 dimensions,
where gauge-invariant winding electric flux strings can be excited by coupling
a chemical potential to each of the global $U(1)$ centre-symmetry generators.
Each sector is labelled by a set of integers
$(\mathbb{Z}_{1},\mathbb{Z}_{2})$, indicating the number of windings in a
specified spatial direction. Moreover, these sectors are topological in
nature, and states in a given winding number sector cannot be smoothly
deformed to another sector. Further, the electric flux tubes are non-local
extended excitations, unlike the point-like bosonic or fermionic particles,
and their properties at finite densities could in principle be considerably
different.

Flux tubes have been played a prominent role in the description of various
physical phenomena. Nielsen and Olesen [8] introduced the field theory of a
vortex-line model, also identified with dual strings. These are flux tubes,
similar to the ones that occur in the theory of type-II superconductors, and
are responsible for most of the low-energy physics in the strong coupling
limit. Classical and semi-classical analysis involving electric fluxes
interacting with a gas of monopoles, giving rise to confinement have been
discussed in [9, 10]. Non-abelian generalizations of such operators, called
disorder operators, were introduced by ’t Hooft to analyse the phases of non-
Abelian gauge theories [11].

We consider the condensed matter physics of these flux tubes in
2+1-dimensional U(1) gauge theory. Previous studies have used the path
integral formulation by either exploiting the dual representation of Abelian
lattice gauge theories [12, 13], or by using the multilevel algorithm [14] and
explored properties such as the profile of the electric flux lines connecting
static charges, or the variation of the potential between two charges with
increasing the representation of the charges. Among other things, this
provides valuable insights about the attractive or repulsive nature of the
flux tubes.

In this article, we use the Hamiltonian formulation of a $U(1)$ quantum link
ladder (QLL) [15]. This theory is known to have novel crystalline confined
phases which carry fractional electric flux excitations [16], possess
anomalously localized excited states [17], and are the building blocks of
spin-ice compounds [18, 19]. While it is known how to simulate the theory with
an improved cluster algorithm at zero and finite temperature [20], this method
has not been extended to deal with the scenario at finite winding chemical
potential. Instead, we use tensor network methods (see for review [21]) to
perform an ab-initio study of the system at finite winding density. Thanks to
the rapid development of quantum simulators, the key elements for realizing
this microscopic model on digital and analogue quantum computers are already
available [22, 23, 24]. The finite density physics investigated in the article
is ideal to be observed in a quantum computing setup. The open boundaries and
gauge invariance realized with quantum spin operators are very natural for
quantum simulators.

$\hat{x}$$\hat{y}$$\circlearrowright$$\circlearrowleft$$\circlearrowright$$\ncirclearrowright$$W_{y}$$W_{x}$

Figure 1: Ladder geometry of the lattice. The periodicity in $\hat{y}$ is
indicated by the dashed lines. Two flippable plaquettes ($\circlearrowleft$,
$\circlearrowright$) and a non flippable plaquette ($\ncirclearrowright$) are
also shown. The dotted lines indicate the links which need to be summed to
obtain the $x$\- and $y$-windings.

## The $U(1)$ quantum link ladder.–

To illustrate our ideas, we consider the setup of the $U(1)$ QLL with the
gauge fields represented by quantum spins in the spin-$\frac{1}{2}$
representation on a rectangular lattice $L_{x}\times L_{y}$, with $L_{y}=2$
and $L_{x}=6,\dots,64$, illustrated in Figure 1. Each link degree of freedom
has a two-dimensional Hilbert space, and the gauge field operator raises (or
lowers) the electric flux basis state:
$U_{r,\hat{i}}=S_{r,\hat{i}}^{+},\leavevmode\nobreak\
U^{\dagger}_{r,\hat{i}}=S_{r,\hat{i}}^{-},\leavevmode\nobreak\
E_{r,\hat{i}}=S_{r,\hat{i}}^{z}$. The Hamiltonian consists of two types of
plaquette operators:

$\mathcal{H}_{\square}=-J\sum_{\square}(U_{\square}+U^{\dagger}_{\square})+\lambda\sum_{\square}(U_{\square}+U^{\dagger}_{\square})^{2}\,,$
(1)

where
$U_{\square}=U_{r,\hat{i}}U_{r+\hat{i},\hat{j}}U^{\dagger}_{r+\hat{j},\hat{j}}U^{\dagger}_{r,\hat{j}}$.
One could have added the square of the electric field energy
$\sum_{r,\hat{i}}E^{2}_{r,\hat{i}}$, but for the spin-$\frac{1}{2}$
representation, this is a trivial constant and can be neglected. As shown in
Figure 1, the first operator flips any flippable plaquette, while the second
operator counts the total number of flippable plaquettes. Only two of the 16
states on a plaquette are non-trivially acted upon by the plaquette operators.
The reduction in the number of physical states is due to a local $U(1)$
symmetry, generated by the Gauss law

$G_{r}=\sum_{\hat{i}=\hat{x},\hat{y}}\left(E_{r-\hat{i},\hat{i}}-E_{r,\hat{i}}\right)=\sum_{\hat{i}=\hat{x},\hat{y}}\left(S^{z}_{r-\hat{i},\hat{i}}-S^{z}_{r,\hat{i}}\right).$
(2)

Physical states satisfy $G_{r}\ket{\psi}=0$, which implies the absence of any
charge on the lattice. In addition, the model has several global symmetries:
the lattice translation symmetry (by one lattice spacing), the reflection and
the rotation symmetry. In addition, there is the $Z_{2}$ charge conjugation
symmetry: $U\rightarrow U^{\dagger},E\rightarrow-E$. However, the main object
of our interest are the $U(1)^{2}$ global winding number symmetries, generated
by the operators:

$W_{x}=\frac{1}{2L_{y}}\sum_{r}S^{z}_{r,\hat{y}}\quad\text{and}\quad
W_{y}=\frac{1}{2L_{x}}\sum_{r}S^{z}_{r,\hat{x}}.$ (3)

where the sum over $r$ runs over all lattice sites. These operators commute
with the Hamiltonian and thus classify the eigenstates in terms of the number
of times the flux loops wind the system either along the $x$\- or the
$y$-direction. Therefore, it is natural to couple chemical potentials with
strengths $\mu_{x},\mu_{y}$ to the Hamiltonian and extend the full Hamiltonian
as: $\mathcal{H}={\mathcal{H}}_{\square}-\mu_{x}W_{x}-\mu_{y}W_{y}$.

The windings $W_{x,y}$ are good quantum numbers for periodic boundary
conditions. However, for using open boundary conditions (as we impose in the
longer directions, since we use matrix product states (MPS) in our
calculations), one can show that an external field $(h_{x},h_{y})$ that
couples to the $x$-links and $y$-links respectively serves the same purpose,
keeping $W_{x,y}$ to be good quantum numbers. With the external field, there
is a non-trivial contribution from the kinetic energy term:
$\sum_{r,\hat{x}}(E_{r,\hat{x}}-h_{x})^{2}+\sum_{r,\hat{y}}(E_{r,\hat{y}}-h_{y})^{2}=-h_{x}\sum_{r,\hat{x}}E_{r,\hat{x}}-h_{y}\sum_{r,\hat{y}}E_{r,\hat{y}}+{\rm
const}=-2h_{x}W_{y}-2h_{y}W_{x}+{\rm const}$, which is equivalent to coupling
the system with $\mu_{x,y}$. We will use the latter notation for the rest of
the article.

$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$1.6$$1.8$$2$$0$$0.2$$0.4$$0.6$$0.8$$1$$\displaystyle\mu_{x}$$\displaystyle\langle
W_{x}\rangle/L_{x}$$L_{x}=4$$L_{x}=8$$L_{x}=16$$L_{x}=32$$L_{x}=64$
Figure 2: Staircase structure of the winding numbers $\langle W_{x}\rangle$
with increasing $\mu$. The plateaux correspond to ground states where the
winding flux remains fixed as $\mu$ is varied. In the thermodynamic limit, the
curve becomes continuous. For large $\mu$, the curve saturates just as for
fermion (or hard-core bosons).

## Numerical methods.–

We begin by noting that the model considered here a rich ground state phase
diagram [16] and realizes novel crystalline confined phases. The physics of
excited states have revealed the existence of quantum scar states, and
atypical real-time dynamics [17]. While the former used an efficient cluster
Monte-Carlo algorithm, the latter used large scale exact diagonalization (ED).
In this work, we aim to go for system sizes beyond the reach of ED, but
efficient algorithms at finite $\mu$ are non-trivial to construct. While the
existing cluster algorithms can update all sectors at finite temperatures, it
is unclear on how to extend this algorithm for finite $\mu$. Therefore, we use
density matrix renormalization group (DMRG) on MPS states to simulate the
ground state phases with increasing values of
$\mu=\sqrt{\mu_{x}^{2}+\mu_{y}^{2}}$. The $\mu_{y}=0$ is kept throughout the
calculations to ensure that there is no condensation of strings in the
$y$-direction.

## Condensation of strings.–

The effect of increasing $\mu_{x}$ on various system sizes is shown in Figure
2. The more familiar examples of condensation phenomena are known from bosons
and fermions, which are point particles. We notice that with flux strings,
too, one has the ”silver blaze” problem, in which the ground state is
unaffected by the chemical potential until a threshold value
$\mu_{\mathrm{c}}(L_{x})$ is reached, after which the vacuum becomes unstable
to the creation of net flux strings periodically winding around $L_{y}$. On
the smaller lattices, one can clearly observe the step-like structure that
results, with each step indicating the number of winding strings that have
condensed in the vacuum. Plotted in terms of the winding density, we notice
the smooth approach to the thermodynamic limit in the data for lattices when
reaching $L_{x}=64$ (see Figure 2). Note in particular that both the threshold
chemical potential, $\mu_{\mathrm{c}}(L_{x})$, at which condensation phenomena
starts and the saturation chemical potential $\mu_{s}(L_{x})$ have well-
defined thermodynamic limits. In Figure 3 we show the behaviour of
$\mu_{\mathrm{c}}(L_{x})$ with increasing volume. It is interesting to note
that the finite volume dependence is very well described with the same formula
that governs the dependence of a massive particle in finite volume [25]. We
note that the step behaviour of magnetization with an external magnetic field
at zero temperature is well known for frustrated spin systems [26]. Recently,
a similar behaviour has been reported for the ladder Rydberg systems [27].

48163264$0.3$$\displaystyle L_{x}$$\displaystyle\mu_{\mathrm{c}}$fitteddata
Figure 3: Finite size dependence of
$\mu_{\mathrm{c}}(L_{x})=a\,\mathrm{exp}(-bL_{x})+\mu_{\mathrm{c}}^{\infty}$
on $L_{x}$. From the fit we determine: $\mu_{\mathrm{c}}^{\infty}=0.269$. The
error bars are the magnitude of the finite step $\Delta_{\mu_{x}}=0.0025$
taken to identify the phase transition point.

While we have demonstrated the thermodynamic limit for $L_{y}=2$ ladders, more
work is essential to extend the results to other geometries. In particular,
the 2-d system can be thought of as a sequence of ladders with increasing
$L_{y}$ at each step. At each fixed $L_{y}$, we can first take the
$L_{x}\to\infty$ limit. Thus, strings that are in general non-local in $L_{y}$
can condense in an infinitely wide ladder. For a confining theory, increasing
$L_{y}\to\infty$ is expected to yield $\mu_{c}(L_{x}\rightarrow\infty,L_{y})$
that increases linearly with $L_{y}$. We postpone the demonstration of the
thermodynamic limit of larger ladders in a future study, and turn to
understanding the nature of the phases that are realized in the ground states
at finite density.

## Crystalline structures.–

As we demonstrate now, once the winding strings start condensing in the ground
state, they modulate existing crystalline properties. At $\mu_{x}=0$ and
$\lambda=-1$, the ground state breaks both translation invariance and charge
conjugation spontaneously [16]. The novel feature at finite $\mu_{x}$ is the
repulsion of condensed strings in the $x$-direction and their subsequent
arrangement in periodic intervals. This necessarily modulates the pattern of
electric fields, $E_{r,\hat{i}}$, from the zero density case.

$-0.25$$0$$0.25$$\displaystyle\langle
E_{x,\hat{i}}\rangle$$\displaystyle\mu_{x}=0.313$$-0.25$$0$$0.25$$\displaystyle\langle
E_{x,\hat{i}}\rangle$$\displaystyle\mu_{x}=1.3$$0$$10$$20$$30$$40$$50$$60$$-0.25$$0$$0.25$$\displaystyle
x$$\displaystyle\langle E_{x,\hat{i}}\rangle$$\displaystyle\mu_{x}=1.91$

Figure 4: Vertical electric field, $\braket{E_{x,\hat{i}}}$, for the
$L_{x}=64$ lattice for three different regimes of winding density. The dashed
lines correspond to the upper rung and solid lines to the electric field on
the lower rung of the ladder

In Figure 4, we show the spatial distribution of the fluxes in the
$y$-direction, $E_{x,\hat{i}}$, at three different $\mu_{x}$ values for the
largest lattice $L_{x}=64$, representative of three distinct regimes. We call
these three different winding regimes: dilute gas regime, half filled and
close to saturation regime.

The first regime occurs when the system has just started to condense isolated
strings, and the system can be treated as a dilute gas of strings. The top
panel of Figure 4 at $\mu_{x}=0.313$ illustrates this case. The three regions
where the $\braket{E_{x,\hat{i}}}\approx 0$ marks the location of the winding
strings wrapping along the $y$-direction. We infer that the preference of the
strings to stay as far away from each other as possible is indicative of their
repulsive interaction. Moreover, in between the location of the fluxes, the
$\braket{E_{x,\hat{i}}}$ displays a regular oscillatory pattern, as also
expected for $\mu_{x}=0$. This arrangement of the fluxes maximizes the total
number of flippable plaquettes, as preferred by the $\lambda=-1$ term in the
Hamiltonian.

On increasing the filling fraction of the winding density, we notice that the
long-wavelength modulations of the electric flux disappear. As shown in the
representative middle panel, for $\mu_{x}=1.3$, the long-range modulations of
$\braket{E_{x,\hat{i}}}$ disappear. The short range oscillations of the
horizontal fluxes are still present with twice the period than the previous
case: the dashed and the solid lines take their maximum positive and negative
values ($\approx\pm 0.25$) 16 times. This regime corresponds to the half-
filling of winding strings, now distributed evenly through the system,
removing traces of previous spatial modulations. Making the system denser
causes one to approach to the saturation regime, where the electric fields
further rearrange to produce a smooth coherent oscillation. The bottom panel
in Figure 4, for $\mu_{x}=1.91$ shows the coherent oscillations for $L_{x}=64$
in this regime, spread over 12-15 lattice spacing.

We can also understand the physical properties from the sum of the electric
fields on the vertical links, $E_{y,\hat{j}}$, which provides the analog of
the ”particle number density”. Following our previous discussions, we also
expect these profiles to show modulations, which are plotted in Figure 5(a)
for our biggest lattice $L_{x}=64$. Three distinct regimes are also visible in
this plot. The set of blue curves represents the regime where the system has
just started to condense isolated strings. It is clear that as $\mu_{x}$ is
slowly increased, the winding strings condense in such a way as to maintain
maximal separation from each other and the highly polarized boundaries. The
first such string excitation sits in the middle of the lattice, as shown by
the maximum in the density profile. The case with three peaks in
$E_{x,\hat{j}}$ (at $x=10a,30a,50a$) correspond to the profile of
$E_{x,\hat{i}}$ at $\mu_{x}=0.313$ shown in Figure 4. The presence of the
fluxes (wrapping vertically) makes the plaquettes non-flippable, which is
exactly the locations where the horizontal fields are minimum and the vertical
fields maximum, demonstrating that the strings affect all the local
properties.

On increasing $\mu_{x}$, the $E_{x,\hat{j}}$ looses the modulations that
identify individual fluxes, and a smooth distribution, modulated more at the
boundaries than in the bulk are visible. Closer to the saturation region upon
further increasing of $\mu_{x}$, again longer ranged smooth modulations of the
”particle density” appear, which now stretch over several lattice spacing.
Interestingly, this length scale seems to be dynamically generated in this
regime and rather sensitive to the external $\mu_{x}$. The wave number of the
oscillations can thus be controlled by tuning the $\mu_{x}$.

(a) $0.9$$1$$0.5$$0.55$particle number$0$$20$$40$$60$$0$$0.1$$\displaystyle x$

(b) $0.5$$1$$1.5$$2$$0$$2$$4$$\displaystyle\mu_{x}$wave
number$\mu_{x}<0.35$$1.2<\mu_{x}<1.3$$1.8<\mu_{x}$

Figure 5: Winding number regimes of the quantum link ladder. (a) The winding
number distribution as a function of the distance to one boundary for the
three different regimes. (b) The wave number as a function of the chemical
potential for the states where the particle number is non-zero or non-
saturated. The three different winding regimes are highlighted with colored
markers/lines: dilute gas regime ($\bullet$), half filled ($\bullet$) and
close to saturation regime ($\bullet$).

Figure 5(b) shows the wave number of the oscillations as a function of
$\mu_{x}$, obtained by identifying the dominant wave-number that contributes
in the Fourier transform of the vertical electric flux profiles,
$E_{x,\hat{j}}$. The information in this observable is thus the same as in the
structure factor up to a global factor, which is given as a Fourier transform
of the electric flux correlation function at a particular momentum $k$ (the
wave number is $k/2\pi$ in our context). Even in this plot, the aforementioned
three regimes in $\mu_{x}$ are clearly visible. The first non-trivial
excitations for small $\mu_{x}$, present long range oscillations whose wave
numbers keep decreasing until they saturate to a small value. This is the
regime where the system is approximately half-filled, and for the $L_{x}=64$
spans from $\mu_{x}=1.2,\dots,1.3$. In this region, the translational
invariance is approximately recovered. When the chemical potential is
increased, the oscillations rise again with a much faster rate, as already
apparent from the earlier observables.

## Conclusions and Outlook.–

In this Letter, we have explored the phenomenon of string condensation in an
$U(1)$ Abelian lattice gauge theory realized as a spin-$1/2$ QLM. We have
demonstrated that our ladder system posses a smooth thermodynamic limit for a
fixed $L_{y}$. The system starts to condense strings with the increase in
$\mu_{x}$, and the system exhibits at least three different regimes before
saturation is reached. Through the profiles of the horizontal and vertical
electric fluxes, we have shown that the winding strings arrange themselves in
patterns which behave distinctly in each of the three regimes. In the dilute
regime, isolated string excitations can be identified, while the half-filled
regime is marked by an approximate restoration of translation invariance. In
the dense region, there is a dynamically generated length scale which changes
rapidly with $\mu_{x}$ before the system saturates. Our observables are
perfectly suited to be measured in cold atom experiments of lattice gauge
theory models [28, 29, 30, 31, 32, 33].

There are several directions in which the analysis can be extended. The most
obvious is to repeat the calculation for larger ladders and study the
different regimes that manifest themselves. Other observables, such as the
central charge, and finite size scaling of correlation functions could be
useful in attempting to understand if there is a phase transition between the
different regimes. The nature of the origin of the length scale in the dense
region is also an open question, which might be understood better from an
effective field theory approach. Another obvious question is if similar
phenomena can also be observed in QLMs in the spin-$1$ representation, which
are very similar to the lattice gauge theory formulation by Wilson.

## Acknowledgements.–

We would like to thank Luca Barbiero, Marcello Dalmonte, Adam Nahum, Arnab
Sen, and Uwe-Jens Wiese for helpful discussions.

PS acknowledges support from: ERC AdG NOQIA; Ministerio de Ciencia y
Innovation Agencia Estatal de Investigaciones
(PGC2018-097027-B-I00/10.13039/501100011033,
CEX2019-000910-S/10.13039/501100011033, Plan National FIDEUA
PID2019-106901GB-I00, FPI, QUANTERA MAQS PCI2019-111828-2, QUANTERA DYNAMITE
PCI2022-132919, Proyectos de I+D+I “Retos Colaboración” QUSPIN
RTC2019-007196-7); European Union NextGenerationEU (PRTR); Fundació Cellex;
Fundació Mir-Puig; Generalitat de Catalunya (European Social Fund FEDER and
CERCA program (AGAUR Grant No. 2017 SGR 134, QuantumCAT U16-011424, co-funded
by ERDF Operational Program of Catalonia 2014-2020); Barcelona Supercomputing
Center MareNostrum (FI-2022-1-0042); EU Horizon 2020 FET-OPEN OPTOlogic (Grant
No 899794); National Science Centre, Poland (Symfonia Grant
No.2016/20/W/ST4/00314); European Union’s Horizon 2020 research and innovation
programme under the Marie-Skłodowska-Curie grant agreement No 101029393
(STREDCH) and No 847648 (“La Caixa” Junior Leaders fellowships ID100010434:
LCF/BQ/PI19/11690013, LCF/BQ/PI20/11760031, LCF/BQ/PR20/11770012,
LCF/BQ/PR21/11840013). PK acknowledges support from the Research Training
Group ”Differential Equation- and Data-driven Models in Life Sciences and
Fluid Dynamics: An Interdisciplinary Research Training Group (DAEDALUS)” (GRK
2433) funded by the German Research Foundation (DFG).

## References

  * Banerjee and Chandrasekharan [2010] D. Banerjee and S. Chandrasekharan, Finite size effects in the presence of a chemical potential: A study in the classical nonlinear o(2) sigma model, Physical Review D 81, 10.1103/physrevd.81.125007 (2010).
  * Aarts and James [2010] G. Aarts and F. A. James, On the convergence of complex langevin dynamics: the three-dimensional XY model at finite chemical potential, Journal of High Energy Physics 2010, 10.1007/jhep08(2010)020 (2010).
  * Katz  _et al._ [2017] S. Katz, F. Niedermayer, D. Nogradi, and C. Torok, Comparison of algorithms for solving the sign problem in the O(3) model in 1+1 dimensions at finite chemical potential, Phys. Rev. D 95, 054506 (2017), arXiv:1611.03987 [hep-lat] .
  * Bloch  _et al._ [2021] J. Bloch, R. G. Jha, R. Lohmayer, and M. Meister, Tensor renormalization group study of the three-dimensional o(2) model, Physical Review D 104, 10.1103/physrevd.104.094517 (2021).
  * Gupta [2010] S. Gupta, QCD at finite density, PoS LATTICE2010, 007 (2010), arXiv:1101.0109 [hep-lat] .
  * Ayyar  _et al._ [2018] V. Ayyar, S. Chandrasekharan, and J. Rantaharju, Benchmark results in the 2D lattice Thirring model with a chemical potential, Phys. Rev. D 97, 054501 (2018), arXiv:1711.07898 [hep-lat] .
  * Bañuls  _et al._ [2017] M. C. Bañuls, K. Cichy, J. I. Cirac, K. Jansen, and S. Kühn, Density Induced Phase Transitions in the Schwinger Model: A Study with Matrix Product States, Phys. Rev. Lett. 118, 071601 (2017), arXiv:1611.00705 [hep-lat] .
  * Nielsen and Olesen [1973] H. Nielsen and P. Olesen, Vortex-line models for dual strings, Nuclear Physics B 61, 45 (1973).
  * Banks  _et al._ [1977] T. Banks, R. Myerson, and J. Kogut, Phase transitions in abelian lattice gauge theories, Nuclear Physics B 129, 493 (1977).
  * Polyakov [1977] A. M. Polyakov, Quark Confinement and Topology of Gauge Groups, Nucl. Phys. B 120, 429 (1977).
  * ’t Hooft [1978] G. ’t Hooft, On the Phase Transition Towards Permanent Quark Confinement, Nucl. Phys. B 138, 1 (1978).
  * Trottier and Woloshyn [1993] H. D. Trottier and R. M. Woloshyn, Flux tubes in three-dimensional lattice gauge theories, Physical Review D 48, 2290 (1993).
  * Zach  _et al._ [1998] M. Zach, M. Faber, and P. Skala, Flux tubes and their interaction in $U(1)$ lattice gauge theory, Nucl. Phys. B 529, 505 (1998), arXiv:hep-lat/9709017 .
  * Koma  _et al._ [2004] Y. Koma, M. Koma, and P. Majumdar, Static potential, force, and flux-tube profile in 4D compact $U(1)$ lattice gauge theory with the multi-level algorithm, Nuclear Physics B 692, 209 (2004).
  * Chandrasekharan and Wiese [1997] S. Chandrasekharan and U. J. Wiese, Quantum link models: A Discrete approach to gauge theories, Nucl. Phys. B 492, 455 (1997), arXiv:hep-lat/9609042 .
  * Banerjee  _et al._ [2013] D. Banerjee, F.-J. Jiang, P. Widmer, and U.-J. Wiese, The (2 + 1)-d $u(1)$ quantum link model masquerading as deconfined criticality, J. Stat. Mech. Theory Exp. 2013, P12010 (2013).
  * Banerjee and Sen [2021] D. Banerjee and A. Sen, Quantum scars from zero modes in an abelian lattice gauge theory on ladders, Phys. Rev. Lett. 126, 220601 (2021).
  * Shannon  _et al._ [2004] N. Shannon, G. Misguich, and K. Penc, Cyclic exchange, isolated states, and spinon deconfinement in an $xxz$ heisenberg model on the checkerboard lattice, Phys. Rev. B 69, 220403 (2004).
  * Benton  _et al._ [2012] O. Benton, O. Sikora, and N. Shannon, Seeing the light: Experimental signatures of emergent electromagnetism in a quantum spin ice, Physical Review B 86, 10.1103/physrevb.86.075154 (2012).
  * Banerjee [2021] D. Banerjee, Recent progress on cluster and meron algorithms for strongly correlated systems, Indian Journal of Physics 95, 1669 (2021).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of physics 326, 96 (2011).
  * Celi  _et al._ [2020] A. Celi, B. Vermersch, O. Viyuela, H. Pichler, M. D. Lukin, and P. Zoller, Emerging two-dimensional gauge theories in rydberg configurable arrays, Phys. Rev. X 10, 021057 (2020).
  * Paulson  _et al._ [2021] D. Paulson, L. Dellantonio, J. F. Haase, A. Celi, A. Kan, A. Jena, C. Kokail, R. van Bijnen, K. Jansen, P. Zoller, and C. A. Muschik, Simulating 2d effects in lattice gauge theories on a quantum computer, PRX Quantum 2, 030334 (2021).
  * Huffman  _et al._ [2021] E. Huffman, M. G. Vera, and D. Banerjee, Real-time dynamics of plaquette models using nisq hardware (2021).
  * Luscher [1986] M. Luscher, Volume Dependence of the Energy Spectrum in Massive Quantum Field Theories. 1. Stable Particle States, Commun. Math. Phys. 104, 177 (1986).
  * Honecker  _et al._ [2004] A. Honecker, J. Schulenburg, and J. Richter, Magnetization plateaus in frustrated antiferromagnetic quantum spin models, Journal of Physics: Condensed Matter 16, S749 (2004).
  * Sarkar  _et al._ [2022] M. Sarkar, M. Pal, A. Sen, and K. Sengupta, Quantum order-by-disorder induced phase transition in rydberg ladders with staggered detuning, arXiv preprint arXiv:2204.12515 (2022).
  * Aidelsburger  _et al._ [2022] M. Aidelsburger, L. Barbiero, A. Bermudez, T. Chanda, A. Dauphin, D. González-Cuadra, P. R. Grzybowski, S. Hands, F. Jendrzejewski, J. Jünemann, G. Juzeliūnas, V. Kasper, A. Piga, S.-J. Ran, M. Rizzi, G. Sierra, L. Tagliacozzo, E. Tirrito, T. V. Zache, J. Zakrzewski, E. Zohar, and M. Lewenstein, Cold atoms meet lattice gauge theory, Philos. Trans. R. Soc. A 380, 20210064 (2022).
  * Martinez  _et al._ [2016] E. A. Martinez _et al._ , Real-time dynamics of lattice gauge theories with a few-qubit quantum computer, Nature 534, 516 (2016), arXiv:1605.04570 [quant-ph] .
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, and et al., Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579–584 (2017).
  * Schweizer  _et al._ [2019] C. Schweizer, F. Grusdt, M. Berngruber, L. Barbiero, E. Demler, N. Goldman, I. Bloch, and M. Aidelsburger, Floquet approach to $\mathbb{Z}_{2}$ lattice gauge theories with ultracold atoms in optical lattices, Nature Physics 15, 1168–1173 (2019).
  * Mil  _et al._ [2020] A. Mil, T. V. Zache, A. Hegde, A. Xia, R. P. Bhatt, M. K. Oberthaler, P. Hauke, J. Berges, and F. Jendrzejewski, A scalable realization of local $U(1)$ gauge invariance in cold atomic mixtures, Science 367, 1128 (2020), arXiv:1909.07641 [cond-mat.quant-gas] .
  * Yang  _et al._ [2020] B. Yang, H. Sun, R. Ott, H.-Y. Wang, T. V. Zache, J. C. Halimeh, Z.-S. Yuan, P. Hauke, and J.-W. Pan, Observation of gauge invariance in a 71-site Bose–Hubbard quantum simulator, Nature 587, 392 (2020), arXiv:2003.08945 [cond-mat.quant-gas] .
  * Fishman  _et al._ [2020] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor software library for tensor network calculations (2020), arXiv:2007.14822 .

##  Appendix A Supplementary Material

## Numerical Implementation.–

For an efficient representation of the ground states in the quantum link
ladder, we use matrix product states implemented in the itensor software
package [34]. All the result presented in this work have been extrapolated to
infinite bond dimension. During the numerical simulations, the bond dimension
$D$ of the matrix product state has been increased up to 1000
($D_{\mathrm{max}}$). We have also checked that the final state can be
compressed with $D<D_{\mathrm{max}}$. The stopping criteria used in the
optimization procedure is that the difference of the energy after 5 sweeps
should be smaller than $10^{-10}$. For the larger volumes and close to the
saturation, in the regime of the $\mu_{x}\sim 1.8$, the numerical simulation
becomes unstable. If we run the optimization procedure different times, we
might find different the ground state in different topological sectors. This
phenomenon is caused because the energy gap between the different topological
sectors approaches zero in the thermodynamic limit and close to the saturation
regime. Since the gap in this region is smaller than the numerical precision,
we could have found inconsistent result for certain parameter values.
Nevertheless, we are still able to capture the important properties and
extract physical information.

## Observables at $\mu_{x}=\mu_{y}=0$.–

In this section, we briefly sketch the physics of the model at zero chemical
potential, but for varying $\lambda$ for completeness. This has already been
discussed at length in [16] for periodic boundary conditions and using an
efficient quantum Monte Carlo algorithm. It was demonstrated in that study
that for large negative values of $\lambda$, where the $J$ term is
insignificant, the ground state physics is dominated by the states with the
largest number of flippable plaquettes. This state spontaneously breaks both
the charge conjugation and the lattice translation symmetry. For decreasing
$\lambda$, however, the symmetry breaking pattern changes as the $J$ term
increases in strength. We encounter a phase where the charge conjugation
symmetry is restored by the lattice translation symmetry remains broken. The
two phases are connected by a weak first order phase transition.

For the quantum link ladder with $L_{y}=2$ and open boundary conditions in the
$x$-direction, the symmetries are different. In particular, the lattice
translation symmetry is not exact any more, while the charge conjugation
symmetry is still exact. Thus, we expect a phase transition from a phase which
breaks charge conjugation symmetry into a phase where the symmetry is
restored. The susceptibility is expected to be large in the symmetry broken
phase and vanish in the phase where the symmetry is restored. Figure 6 shows
the expected behaviour of the susceptibility as a function of $\lambda$ for
three different lattice sizes.

$-1.2$$-1.1$$-1$$-0.9$$-0.8$$-0.7$$-0.6$$-0.5$$-0.4$$-0.3$$-0.2$$-0.1$$0$$0$$5$$10$$15$$20$$25$$30$$\displaystyle\lambda$$\displaystyle\chi$$L_{x}=16$$L_{x}=24$$L_{x}=32$
Figure 6: Susceptibility as a function of $\lambda$ at $\mu_{x}=0$. The phase
for large negative lambda has a broken symmetry, indicated by the increase in
$\chi$ with $L_{x}$. In the symmetry restored phase, $\chi$ does not increase
with the system size as expected. $0$$10$$\displaystyle
W_{x}$$L_{x}=16$$-20$$0$$\displaystyle\partial
E_{0}/\partial\mu_{x}$0.340.460.60.740.8611.121.261.361.48$0$$0.5$$\displaystyle\mu_{x}$$\displaystyle
M$ Figure 7: (Top) Winding number sectors for $L_{x}=16$. (Middle) The
derivative of the ground state energy with respect to the chemical potential.
(Bottom) The magnetization, $M$, defined as the difference of the plaquettes
flippable in the clockwise and anti-clockwise fashion respectively.

## Observables at $\mu_{x}>0$.–

As the $\mu_{x}$ is cranked up (keeping $\mu_{y}=0$), we study the behaviour
of the expectation of the winding number $\braket{W_{x}}$, and the derivative
of the ground state energy $E_{0}$ with respect to $\mu_{x}$, $\frac{\partial
E_{0}}{\partial\mu_{x}}$. In the top and the middle panel in Figure 7, we plot
the two quantities. As expected from the Feynman-Hellman theorem, these two
quantities show the same qualitative behaviour upto a constant factor (with an
overall negative sign). Finally, we also plot the flippability of the
plaquettes, $M$. The flippability is defined as the difference between the
total number of plaquettes which are flippable in the clockwise fashion and
the ones which are flippable in the anti-clockwise fashion. Interestingly,
this observable shows a staggered behaviour between zero and non-zero values
for odd and even winding number of electric fluxes respectively. This is a
clear indication of a co-operative behaviour of the plaquettes across the
lattice.