# Analytic and Gevrey class regularity for parametric elliptic eigenvalue
problems

Alexey Chernov  Institut für Mathematik, Carl von Ossietzky Universität
Oldenburg Campus Oldenburg-Wechloy, Carl von Ossietzky Straße 9-11, 26129
Oldenburg<EMAIL_ADDRESS>and  Tung Le  Institut für
Mathematik, Carl von Ossietzky Universität Oldenburg<EMAIL_ADDRESS>Dedicated to Professor Christoph Schwab on the occasion of his 60th birthday

(Date: May 25, 2023 and, in revised form …)

###### Abstract.

We investigate a class of parametric elliptic eigenvalue problems with
homogeneous essential boundary conditions where the coefficients (and hence
the solution $u$) may depend on a parameter $\boldsymbol{y}$. For the
efficient approximate evaluation of parameter sensitivities of the first
eigenpairs on the entire parameter space we propose and analyse Gevrey class
and analytic regularity of the solution with respect to the parameters. This
is made possible by a novel proof technique which we introduce and demonstrate
in this paper. Our regularity result has immediate implications for
convergence of various numerical schemes for parametric elliptic eigenvalue
problems, in particular, for elliptic eigenvalue problems with infinitely many
parameters arising from elliptic differential operators with random
coefficients.

###### Key words and phrases:

Elliptic eigenvalue problems, parametric regularity analysis, numerical
integration

###### 2020 Mathematics Subject Classification:

Primary 65N25, 65C30, 65D30, 65N30

##  1\. Introduction

Multiparametric eigenvalue problems arise in numerous applications in
engineering and physics. For the perturbation theory we refer to the works by
Rellich and Kato [18, 22] and more recent works by Andreev and Schwab [1],
Gilbert et al. [12, 13], Grubiŝić et al. [15], and Dölz and Ebert [9] where
the particular case of stochastic parameter perturbation has been addressed.
Irrespectively from the nature of the perturbation the study of the regularity
of the solution is important since it allows to design appropriate
approximation schemes in the parameter domain, tailored to a specific
regularity class. For example, articles [12, 13] propose efficient Quasi-Monte
Carlo approximations for EVP with stochastic coefficients. The corresponding
regularity analysis is a part of the convergence proof.

In this paper, we address the regularity for the solution of real second-order
elliptic eigenvalue problems (EVP) in the general form

(1.1)
$\begin{split}-\nabla\cdot(\tilde{a}(\boldsymbol{x},\boldsymbol{y})\nabla
u(\boldsymbol{x},\boldsymbol{y}))+b(\boldsymbol{x},\boldsymbol{y})u(\boldsymbol{x},\boldsymbol{y})&=\lambda(\boldsymbol{y})c(\boldsymbol{x},\boldsymbol{y})u(\boldsymbol{x},\boldsymbol{y})\qquad(\boldsymbol{x},\boldsymbol{y})\in
D\times U,\\\ u(\boldsymbol{x},\boldsymbol{y})&=0\hskip
93.89409pt(\boldsymbol{x},\boldsymbol{y})\in\partial D\times U,\end{split}$

where the derivative operator $\nabla$ acts in the physical variable
$\boldsymbol{x}\in D$, where $D\subset\mathbb{R}^{d}$ is a bounded Lipschitz
domain. The vector of parameters $\boldsymbol{y}=(y_{1},y_{2},\dots)\in U$ has
either finitely many or countably many components. For example, if
$\boldsymbol{y}$ is a random parameter, the model with
$U:=[-\frac{1}{2},\frac{1}{2}]^{\mathbb{N}}$ and $\boldsymbol{y}\in U$ being a
countably-dimensional vector of independently and identically distributed
uniform random variables has been frequently used in the literature [5, 6, 20,
21]. We mention that for our analysis boundedness of the parameter domain is
necessary as it implies the positivity of the relative spectral gap introduced
in Section 4 below.

It turns out that it is more convenient to work with the rescaled version of
the diffusion coefficient

(1.2)
$a(\boldsymbol{x},\boldsymbol{y})=\chi_{1}\tilde{a}(\boldsymbol{x},\boldsymbol{y}),$

where $\chi_{1}$ is the smallest eigenvalue of the deterministic Dirichlet-
Laplace eigenvalue problem in the domain $D$

(1.3) $\chi_{1}=\inf_{v\in H^{1}_{0}(D)}\dfrac{\int_{D}|\nabla
u|^{2}\,dx}{\int_{D}|u|^{2}\,dx}.$

On the one hand, for a fixed domain $D$ the eigenvalue $\chi_{1}$ is a fixed
positive number and therefore the rescaling (1.2) is no restriction of
generality. On the other hand, this trivial rescaling helps to simplify the
expressions in the forthcoming proofs. Moreover, we may argue that $\tilde{a}$
is the only parameter in (1.1) which depends as $\mathcal{L}^{2}$ on the
length units $\mathcal{L}$, whereas by (1.3) $\chi_{1}$ depends on
$\mathcal{L}^{-2}$. In this sense, the rescaling (1.2) is natural, since it
renders $a(\boldsymbol{x},\boldsymbol{y})$ invariant w.r.t. the change of
length units.

In the following we assume that the coefficients $a$, $b$, $c$ admit the
uniform bounds

(1.4) $\frac{\overline{a}}{2}\geq
a(\boldsymbol{x},\boldsymbol{y})\geq\underline{a}>0,\qquad\frac{\overline{b}}{2}\geq
b(\boldsymbol{x},\boldsymbol{y})\geq 0,\qquad\frac{\overline{c}}{2}\geq
c(\boldsymbol{x},\boldsymbol{y})\geq\underline{c}>0$

for all $\boldsymbol{y}\in U$ and almost all $\boldsymbol{x}\in D$. In this
case for every fixed $\boldsymbol{y}\in U$ the eigenvalues of (1.1) are real,
positive, bounded from below, and the smallest eigenvalue $\lambda_{1}$ is
simple, see [11, Theorem 8.37] and [8, 16]. Since the coefficients depend on
the parameters $\boldsymbol{y}$, the eigenvalues $\lambda$ and corresponding
eigenfunctions $u$ will depend on $\boldsymbol{y}$ as well. Particularly, if
$\boldsymbol{y}$ is random, then $\lambda(\boldsymbol{y})$ and
$u(\boldsymbol{x},\boldsymbol{y})$ will be random too.

In this paper we present a rigorous regularity analysis for the smallest
eigenvalue $\lambda_{1}$ and the corresponding eigenfunction $u_{1}$ with
respect to the parameter $\boldsymbol{y}$ in the general case where the
coefficients $a,b,c$ are infinitely differentiable functions of
$\boldsymbol{y}$ belonging to the Gevrey class $G^{\delta}$ for some fixed
$\delta\geq 1$. The scale of Gevrey classes is a nested scale of the parameter
$\delta$ that fills the gap between analytic and $C^{\infty}$ functions

$G^{\delta}\subset G^{\delta^{\prime}}\subset C^{\infty},\qquad
1\leq\delta<\delta^{\prime}.$

The class of analytic functions corresponding to the special case $G^{\delta}$
with $\delta=1$, is the smallest and arguably most important of this scale (we
refer to Section 2 for precise definitions) and has been addressed for
parametric/stochastic EVP before. The work [1] proves the analyticity of the
smallest eigenpair using elegant complex analysis arguments. More recently
Gilbert et al. have used inductive real-variable arguments to derive explicit
upper bounds for all $\boldsymbol{y}$-derivatives of the eigenpair that are
close to bounds for real analytic functions, but are still include
$\epsilon$-suboptimal terms in the factorial term. Moreover, all above
mentioned works consider the special case of the affine parametrization of the
coefficients of the type

(1.5)
$a(\boldsymbol{x},\boldsymbol{y})=a_{0}(\boldsymbol{x})+\sum_{m=1}^{\infty}y_{m}a_{m}(\boldsymbol{x}).$

Let us also mention a specific periodic model

(1.6)
$a(\boldsymbol{x},\boldsymbol{y})=a_{0}(\boldsymbol{x})+\sum_{m=1}^{\infty}\sin(2\pi
y_{m})a_{m}(\boldsymbol{x})$

introduced and analysed for parametric elliptic source problems in a recent
work [17].

In this paper we go beyond these settings in various ways.

  1. (1)

We demonstrate that the reason for the suboptimality in the real-variable
inductive proofs is indeed an artefact of the proof strategy. We introduce a
modified argument and obtain optimal upper bounds for the
$\boldsymbol{y}$-derivatives of the smallest eigenvalue and the corresponding
eigenfunction.

  2. (2)

We go beyond affine coefficient expansions as in (1.5) or the periodic
expansion (1.6) and allow for general analytic ($\delta=1$) and Gevrey class
($\delta\geq 1$) dependency on $\boldsymbol{y}$ in all coefficients $a,b,c$ in
(1.1).

  3. (3)

The case of Gevrey class non-analytic coefficients ($\delta>1$) has not been
considered in the existing literature (notice that complex analysis arguments
cannot be be applied in this case). On the other hand, it greatly extends the
palette of parametrizations that can be used to model various phenomena,
covering e.g. partition-of-unity parametrizations. We found that our proof is
essentially the same as in the analytic case ($\delta=1$) and therefore cover
Gevrey-$\delta$ class parametrizations in our proofs.

The proofs of new, more general regularity results become possible due to the
novel proof technique that we call the _alternative-to-factorial technique_.
The aim of this paper is to introduce this approach and demonstrate it in the
example of the parametric eigenvalue problems. However, our new proof
technique is not limited to the class of EVPs considered here and can be
applied for more general nonliear parametric PDEs, and hence this is a subject
of our current research.

The structure of the paper is as follows. In Section 2 we demonstrate the
reasons for deficiency of the real-value inductive arguments from the existing
literature and explain how our novel _alternative-to-factorial technique_
helps to overcome the difficulties. In Section 3 we formulate the regularity
assumptions on the coefficients of the eigenvalue problem (1.1) and state our
main theorem claiming the same type regularity for the smallest eigenvalue and
the corresponding eigenfunction. In Section 4 we summarize the properties of
elliptic eigenvalue problems needed for the forthcoming regularity analysis.
In Section 5 we present the proof of the main result. The meaning and validity
of the main regularity result is illustrated by numerical experiments in
Section 6. The appendix contains two technical combinatorial results required
in the proofs.

In that follows we denote by $\mathbb{N}$ the set of positive and by
$\mathbb{N}_{0}$ of nonnegative integers. By $L^{2}(D)$ and $L^{\infty}(D)$ we
denote the spaces of square integrable and bounded functions equipped with
standard norms. $V:=H^{1}_{0}(D)$ is the subspace of $L^{2}(D)$ with square
integrable partial derivatives and vanishing Dirichlet trace equipped with the
norm

(1.7) $\|v\|_{V}:=\chi_{1}^{-1/2}\|\nabla v\|_{L^{2}(D)},$

where $\chi_{1}$ is the smallest eigenvalue of the deterministic Dirichlet-
Laplace problem (1.3). Similarly as in (1.2), the multiplier $\chi_{1}$ is
added here for scaling reasons. This will help us to simplify the expressions
throughout the paper. We will also occasionally use the short notation
$\left\|{\cdot}\right\|_{\infty}$ for the essential supremum norm
$\left\|{\cdot}\right\|_{L^{\infty}(D)}$ in the physical domain $D$.

##  2\. Preliminaries

###  2.1. Deficiency of the real-variable inductive argument

The deficiency of the real-variable inductive argument for nonlinear problems
can be seen already in the one-dimensional case. It is in fact a consequence
of the Leibniz product rule and the triangle inequality. To illustrate this we
collect some elementary results on real analytic functions on the real line,
see e.g. [19], Section 1.

######  Definition 2.1.

Let $U\subseteq\mathbb{R}$ be an open domain. A function $f$ is called real
analytic at $y_{0}\in U$, if it may be represented by a convergent power
series

(2.1) $f(y)=\sum_{j=0}^{\infty}a_{j}(y-y_{0})^{j}$

in some neighbourhood of $y_{0}$. The function $f$ is called real analytic in
an open subset $I\subseteq U$, if it is analytic in all $y_{0}\in I$.

The uniqueness of the power series implies the identity
$a_{n}=f^{(n)}(y_{0})/n!$ and (after some additional work) the following
alternative characterization of real analytic functions.

######  Theorem 2.2.

Let $f\in C^{\infty}(I)$ for some open interval $I$. The function $f$ is real
analytic if and only if for each $y_{0}\in I$, there are an open interval $J$,
with $y_{0}\in J\subseteq I$, and constants $R>0$ and $C>0$ such that the
derivatives of $f$ satisfy

(2.2) $\frac{|f^{(n)}(y)|}{n!}\leq\frac{C}{R^{n}},\qquad\forall y\in J,\quad
n\in\mathbb{N}_{0}.$

The radius of convergence $\rho$ of the power series (2.1) at some $y_{0}\in
I$ can be determined as the supremum of $R$ such that (2.2) holds with some
$C=C(R)$.

Let now $f$ be some analytic function normalized to $C=1$, centred at $y$ with
radius of convergence $\rho_{f}$ and consider $g=f^{2}$. Clearly, $g$ is
analytic, but what is its radius of convergence $\rho_{g}$? If we use the
Leibniz product rule, the triangle inequality and (2.2), we get

(2.3)
$\frac{|g^{(n)}(y)|}{n!}=\frac{1}{n!}\sum_{m=0}^{n}\binom{n}{m}|f^{(m)}(y)|\,|f^{(n-m)}(y)|\leq\frac{n+1}{R^{n}}.$

since

(2.4) $\sum_{m=0}^{n}\binom{n}{m}m!(n-m)!=(n+1)n!$

Observe that the upper bound has increased by the factor $(n+1)$ and therefore
does not fit the upper bound (2.2) with the same $R$. One way to proceed is to
use the crude estimate $n+1\leq 2^{n}$ which holds for all
$n\in\mathbb{N}_{0}$, but according to the upper bound this implies that the
radius of convergence $\rho_{g}$ of $g$ degrades to $\rho_{f}/2$. This kind of
argument can clearly not be used in inductive proofs, since when $n\to\infty$
the radius of convergence will degenerate to zero. A more successful option is
to replace the term $n!$ by $(n!)^{1+\epsilon}$ for $\epsilon>0$, see e.g. the
work by Gilbert et al. in [12] for further details. In this case
$\epsilon$-suboptimal terms pollute the upper bounds and analytic-type
estimates as in (2.2) cannot be obtained for $g$.

On the other hand, it is the classic result that the function $g=f^{2}$ is
analytic with _the same radius of convergence_ as the function $f$ itself.
This can either be seen from the root test applied to (2.3) or the following
explicit bound for the remainder of the power series of $g$. Let
$S_{N}(y)=\sum_{m=0}^{N}b_{m}(y-y_{0})^{m}$ be the partial sum of $g$. Then
$b_{m}=\sum_{j+k=m}a_{j}a_{k}$ and there holds

$|g(y)-S_{2N+1}(y)|=\left|\sum_{m=2N+2}^{\infty}\sum_{j+k=m}a_{j}a_{k}(y-y_{0})^{j+k}\right|\leq
2\sum_{j=N+1}^{\infty}|a_{j}|\,|y-y_{0}|^{j}\sum_{k=0}^{\infty}|a_{k}|\,|y-y_{0}|^{k}$

Since the power series of $f$ converges uniformly and absolutely for
$|y-y_{0}|<\rho_{f}$, the $k$-sum is bounded and the $j$-sum tends to zero for
$N\to\infty$ and all $|y-y_{0}|\leq\rho_{f}$. Hence $\rho_{g}=\rho_{f}$. This
demonstrates that inductive proofs based on the Leibniz product rule, the
triangle inequality in combination with (2.2) can hardly lead to optimal
regularity estimates for nonlinear parametric problems.

###  2.2. The falling factorial estimates as a remedy

Observe that the quadratic type nonlinearity $g=f^{2}$ discussed in the
previous section is the simplest, yet very important type of nonlinearity
appearing e.g. in Navier-Stokes equations and eigenvalue problems that serve
as a model for this paper. The use of the Leibniz product rule together with
the triangle inequality in the parametric regularity proofs is a desirable
tool that splits the nonlinear part into basic contributions and therefore
should be preserved. _Instead we seek for an alternative equivalent for ( 2.2)
that can be used in the course of induction_. To get a better insight, let us
understand whether there is a function $f$ defined in a neighbourhood $J$ of
$y_{0}=0$, such that

(2.5) $|g^{(n)}(y)|\leq|f^{(n)}(y)|\qquad\forall y\in J,\quad
n\in\mathbb{N}_{0}$

is valid for $g=f^{2}$? The answer is positive and an example of such a
function $f$ can be given in closed form, that is

(2.6) $f(y)=\tfrac{1}{2}(1-\sqrt{1-y}).$

Indeed, $g(y)=f(y)-\frac{y}{4}$ and hence (2.5) holds as equality for $n\geq
2$. Basic calculations also verify (2.5) for $n=0$ and $1$ and any
$J\subset[-3,1)$. The derivatives of $f$ are all positive for $y<1$ and
satisfy

(2.7)
$2f^{(n)}(y)=(-1)^{n-1}\left(\tfrac{1}{2}\right)\left(\tfrac{1}{2}-1\right)\dots\left(\tfrac{1}{2}-n+1\right)(1-y)^{\tfrac{1}{2}-n}=(-1)^{n-1}\left(\tfrac{1}{2}\right)_{n}(1-y)^{\tfrac{1}{2}-n}$

for $n\geq 1$ and the _falling factorial_ notation

(2.8) $(q)_{n}:=\left\\{\begin{array}[]{cl}1,&n=0,\\\
q(q-1)\dots(q-n+1),&n\geq 1,\end{array}\right.\qquad q\in\mathbb{R},\quad
n\in\mathbb{N}_{0}.$

For $q<1$ the falling factorial $(q)_{n}$ is a sign-alternating sequence of
$n$. As already indicated by (2.6) and (2.7), for our purposes the particular
value $q=\frac{1}{2}$ will be important. To further simplify the notation and
avoid keeping track of the sign alteration, we denote the _absolute value of
the falling factorial of $\frac{1}{2}$_ by

$\left[\tfrac{1}{2}\right]_{n}:=\left|\left(\tfrac{1}{2}\right)_{n}\right|.$

This notation appears somewhat non-standard, but quite convenient, as we will
see in the sequel. Notice that the defining property of the Gamma function
$\Gamma(z+1)=z\Gamma(z)$ and $\Gamma(\tfrac{1}{2})=\sqrt{\pi}$ imply the
representation for all nonnegative integers $n\in{\mathbb{N}}_{0}$

(2.9)
$\left[\tfrac{1}{2}\right]_{n}=\frac{1}{2\sqrt{\pi}}|\Gamma(n-\tfrac{1}{2})|\quad\text{or
equivalently}\quad
n!=\xi_{n}\left[\tfrac{1}{2}\right]_{n}\quad\text{with}\quad\xi_{n}=2\sqrt{\pi}\frac{\Gamma(n+1)}{|\Gamma(n-\tfrac{1}{2})|}$

and the trivial bound $1\leq\xi_{n}\leq 2\cdot 2^{n}$ implying

(2.10) $[\tfrac{1}{2}]_{n}\leq n!\leq 2\cdot 2^{n}[\tfrac{1}{2}]_{n}.$

Observe next that $g^{(n)}(0)=f^{(n)}(0)$ holds for $n\geq 2$ and therefore
this value can be computed in two different ways: On the one hand, it
satisfies (2.7), i.e. $g^{(n)}(0)=\tfrac{1}{2}\left[\tfrac{1}{2}\right]_{n}$.
On the other hand, by the Leibniz product rule

(2.11)
$g^{(n)}(0)=\sum_{i=0}^{n}\binom{n}{i}f^{(i)}(0)f^{(n-i)}(0)=\frac{1}{4}\sum_{i=1}^{n-1}\binom{n}{i}\left[\tfrac{1}{2}\right]_{i}\,\left[\tfrac{1}{2}\right]_{n-i},\qquad
n\geq 2$

where we used $f(0)=0$ in the last step. Since by definition
$\left[\tfrac{1}{2}\right]_{0}=1$, this implies the following combinatorial
identities.

######  Lemma 2.3.

For all integers $n\geq 2$ the following identities hold

$\sum_{i=1}^{n-1}\binom{n}{i}\left[\tfrac{1}{2}\right]_{i}\,\left[\tfrac{1}{2}\right]_{n-i}=2\left[\tfrac{1}{2}\right]_{n},\quad\sum_{i=1}^{n}\binom{n}{i}\left[\tfrac{1}{2}\right]_{i}\,\left[\tfrac{1}{2}\right]_{n-i}=3\left[\tfrac{1}{2}\right]_{n},\quad\sum_{i=0}^{n}\binom{n}{i}\left[\tfrac{1}{2}\right]_{i}\,\left[\tfrac{1}{2}\right]_{n-i}=4\left[\tfrac{1}{2}\right]_{n}.$

With the convention that the empty sum equals zero the statement extends to
all integers $n\geq 0$ as

$\sum_{i=1}^{n-1}\binom{n}{i}\left[\tfrac{1}{2}\right]_{i}\,\left[\tfrac{1}{2}\right]_{n-i}\leq
2\left[\tfrac{1}{2}\right]_{n},\quad\sum_{i=1}^{n}\binom{n}{i}\left[\tfrac{1}{2}\right]_{i}\,\left[\tfrac{1}{2}\right]_{n-i}\leq
3\left[\tfrac{1}{2}\right]_{n},\quad\sum_{i=0}^{n}\binom{n}{i}\left[\tfrac{1}{2}\right]_{i}\,\left[\tfrac{1}{2}\right]_{n-i}\leq
4\left[\tfrac{1}{2}\right]_{n}.$

Remarkably, the constants in the right hand sides of the above identities are
uniformly bounded, whereas the constant in the right hand side of (2.4) grows
as $(n+1)$. This suggests to replace the term $n!$ by the falling factorial
$[\frac{1}{2}]_{n}$ in the inductive proof. Moreover, there holds the
following alternative characterization of analytic functions.

######  Theorem 2.4.

Let $f\in C^{\infty}(I)$ for some open interval $I$. The function $f$ is real
analytic if and only if for each $y_{0}\in I$, there are an open interval $J$,
with $y_{0}\in J\subseteq I$, and constants $C>0$ and $R>0$ such that the
derivatives of $f$ satisfy

(2.12)
$|f^{(n)}(y)|\leq\frac{C}{R^{n}}\left[\tfrac{1}{2}\right]_{n},\qquad\forall
y\in J,\quad n\in\mathbb{N}_{0}.$

The radius of convergence $\rho$ of the power series (2.1) at some $y_{0}\in
I$ can be determined as the supremum of $R$ satisfying (2.12).

###### Proof.

Let us fix $y_{0}\in I$. By Theorem 2.2 the radius of convergence of the power
series (2.1) at $y_{0}$ can be written as

$\rho=\sup\\{R\leavevmode\nobreak\ :\leavevmode\nobreak\ R\text{ satisfies
\eqref{ana-bound-1d} with some constant }C=C(R)\\}.$

Thus, if we denote

$\tilde{\rho}=\sup\\{R\leavevmode\nobreak\ :\leavevmode\nobreak\ R\text{
satisfies \eqref{ff-bound-1d} with some constant }C=C(R)\\}$

our aim would be to show that $\rho=\tilde{\rho}$. On the one hand, there
holds $\left[\tfrac{1}{2}\right]_{n}\leq n!$. Together with (2.12) this
implies (2.2) with the same $R$ and $C$ and thereby $\rho\leq\tilde{\rho}$. On
the other hand, if (2.2) holds, i.e. $f$ is analytic at $y_{0}\in I$, the root
test shows that the second derivative $f^{\prime\prime}$ in analytic at
$y_{0}$ too and moreover has the same radius of convergence $\rho$ as the
original function $f$, see e.g. Proposition 1.1.13 in [19]. Now let $n\geq 2$.
By Theorem 2.2 this means

(2.13)
$|f^{(n)}(y)|=|(f^{\prime\prime}(y))^{(n-2)}(y)|\leq\frac{C}{R^{n-2}}(n-2)!$

From definition (2.9) of the falling factorial we find $(n-2)!\leq
4\left[\tfrac{1}{2}\right]_{n}$ and hence (2.13) implies (2.12) with the same
$R$ and $\widetilde{C}:=4CR^{2}$

$|f^{(n)}(y)|\leq\frac{\widetilde{C}}{R^{n}}\left[\tfrac{1}{2}\right]_{n},\qquad\forall
y\in J,\quad n\in\mathbb{N}_{0},\quad n\geq 2.$

By a suitable adjustment of $\widetilde{C}$ the range of $n$ can be extended
to cover the cases $n=0$ and $n=1$. In other words, (2.2) implies (2.12) with
the same $R$ and possibly adjusted $C$. This implies $\rho\geq\tilde{\rho}$
and thereby finishes the proof. ∎

###  2.3. Gevrey-class and analytic functions

The following standard multiindex notations, see e.g. [2, 5], will be used in
what follows. We denote the countable set of finitely supported sequences of
nonnegative integers by

(2.14)
$\mathcal{F}:=\left\\{\boldsymbol{\nu}=(\nu_{1},\nu_{2},\dots)\leavevmode\nobreak\
:\leavevmode\nobreak\ \nu_{j}\in\mathbb{N}_{0},\text{ and }\nu_{j}\neq 0\text{
for only a finite number of }j\right\\}$

where the summation $\boldsymbol{\alpha}+\boldsymbol{\beta}$ and the partial
order relations $\boldsymbol{\alpha}<\boldsymbol{\beta}$ and
$\boldsymbol{\alpha}\leq\boldsymbol{\beta}$ of elements in
$\boldsymbol{\alpha},\boldsymbol{\beta}\in\mathcal{F}$ are understood
componentwise. We write

$\displaystyle\left|\boldsymbol{\nu}\right|:=\sum_{j\geq
1}\nu_{j},\qquad\qquad\boldsymbol{\nu}!:=\prod_{j\geq
1}v_{j}!,\qquad\qquad\boldsymbol{R}^{\boldsymbol{\nu}}=\prod_{j\geq
1}R_{j}^{v_{j}}$

for the absolute value, the multifactorial and the power with the multi-index
$\boldsymbol{\nu}$ and a sequence $\boldsymbol{R}=\\{R_{j}\\}_{j\geq 1}$ of
positive real numbers. For two multiindices
$\boldsymbol{m},\boldsymbol{\nu}\in\mathcal{F}$ with
$\boldsymbol{m}\leq\boldsymbol{\nu}$ the binomial coefficient is defined by

$\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)=\frac{\boldsymbol{\nu}!}{(\boldsymbol{\nu}-\boldsymbol{m})!\boldsymbol{m}!}.$

Notice that $\left|\boldsymbol{\nu}\right|$ is finite if and only if
$\boldsymbol{\nu}\in\mathcal{F}$. For $\boldsymbol{\nu}\in\mathcal{F}$
supported in $\left\\{1,2,\dots,n\right\\}$, we define the partial derivative
with respect to the variables $\boldsymbol{y}$

$\displaystyle\partial^{\boldsymbol{\nu}}u=\frac{\partial^{\left|\boldsymbol{\nu}\right|}u}{\partial
y_{1}^{\nu_{1}}\partial y_{2}^{\nu_{2}}\dots\partial y_{n}^{\nu_{n}}}.$

The above multiindex notations are helpful for treatment of multiparametric
objects. In particular, our regularity considerations are based on the
following definition of Gevrey-$\delta$ functions with countably many
parameters.

######  Definition 2.5.

Let $\delta\geq 1$, $B$ be a Banach space, $I\subset\mathbb{R}^{\mathbb{N}}$
be an open domain and a function $f:I\to B$ be such that its
$\boldsymbol{y}$-derivatives $\partial^{\boldsymbol{\nu}}f:I\to B$ are
continuous for all $\boldsymbol{\nu}\in\mathcal{F}$. We call the function $f$
Gevrey-$\delta$ if for each $y_{0}\in I$ there exist an open set $J\subseteq
I$ and $y_{0}\in J$ and strictly positive constants
$\boldsymbol{R}=(R_{1},R_{2},\dots)\subset\mathbb{R}^{\mathbb{N}}$ and
$C\in\mathbb{R}$ that the derivatives of $f$ satisfy the bounds

(2.15)
$\|\partial^{\boldsymbol{\nu}}f(\boldsymbol{y})\|_{B}\leq\frac{C}{\boldsymbol{R}^{\boldsymbol{\nu}}}(|\boldsymbol{\nu}|!)^{\delta},\qquad\forall\boldsymbol{y}\in
J,\quad\forall\boldsymbol{\nu}\in\mathcal{F}.$

In this case we write $f\in G^{\delta}(U,B)$.

######  Remark 2.6.

For the case of finitely many parameters $\boldsymbol{y}=(y_{1},\dots,y_{M})$
the definition of Gevrey-$\delta$ functions with values in $B=\mathbb{R}$ or
$\mathbb{C}$ can e.g. be found in [23, Section 1.4] as well as some further
characterizations of this class in the case $M<\infty$. However, to the best
of authors’ knowledge, characterizations of the type (2.15) for $M=\infty$,
likewise the related characterization (3.2) below (neither for $M<\infty$ nor
$M=\infty$, or in the special case $\delta=1$) have not been discussed in the
literature before.

######  Remark 2.7.

Definition 2.5 is also suitable for the case of finitely many parameters. In
particular, when $\boldsymbol{y}=(y_{1},\dots,y_{M})$, $B=\mathbb{R}$ or
$\mathbb{C}$ and $\delta=1$, the bound (2.15) guarantees convergence of the
power series of $f$ and therefore characterizes the class of analytic
functions of $M$ variables, see e.g. in [19, Section 2.2]. This requires the
bound $|\boldsymbol{\nu}|!\leq M^{|\boldsymbol{\nu}|}\boldsymbol{\nu}!$ that
is valid for a multiindex $\boldsymbol{\nu}$ with $M$ nonzero components.

######  Remark 2.8.

In general, the estimate (2.15) does not guarantee convergence of the power
series. Indeed, in the case $\delta=1$ and countably many parameters
$\boldsymbol{y}=(y_{1},y_{2},\dots)$, we need an additional assumption, e.g.
that the sequence $\boldsymbol{R}=(R_{1},R_{2},\dots)$ grows such that
$\sum_{i=1}^{\infty}R_{i}^{-1}<\infty$, cf. [12]. To see this observe that for
all $\boldsymbol{y}\in U$ from the neighbourhood
$\|\boldsymbol{y}-\boldsymbol{y}_{0}\|_{\infty}<r$ (2.15) with $\delta=1$
yields

$\displaystyle|f(\boldsymbol{y})|=\left|\sum_{\boldsymbol{\nu}\in\mathcal{F}}\frac{\partial^{\boldsymbol{\nu}}f(\boldsymbol{y}_{0})}{\boldsymbol{\nu}!}(\boldsymbol{y}-\boldsymbol{y}_{0})^{\boldsymbol{\nu}}\right|\leq
C\sum_{\boldsymbol{\nu}\in\mathcal{F}}\frac{\left|\boldsymbol{\nu}\right|!}{\boldsymbol{\nu}!}\left(\frac{r}{\boldsymbol{R}}\right)^{\boldsymbol{\nu}}=\frac{C}{1-r\sum_{i=1}^{\infty}R_{i}^{-1}}.$

where we have used the identity
$\left(1-\sum_{i=1}^{\infty}z_{i}\right)^{-1}=\sum_{\boldsymbol{\nu}\in\mathcal{F}}\frac{\left|\boldsymbol{\nu}\right|!}{\boldsymbol{\nu}!}\boldsymbol{z}^{\boldsymbol{\nu}}$
in the last step.

##  3\. Regularity assumptions and the formulation of the main result

We are now ready to formulate the Gevrey-class regularity assumption on the
coefficients $a,b,c$ of the eigenvalue problem (1.1) which (as we prove in the
following) will imply the Gevrey-class regularity of the eigenpairs with the
same $\delta\geq 1$.

######  Assumption 3.1.

For all values $\boldsymbol{y}\in U$ and $\boldsymbol{\nu}\in\mathcal{F}$, the
coefficients
$a(\boldsymbol{x},\boldsymbol{y}),b(\boldsymbol{x},\boldsymbol{y}),c(\boldsymbol{x},\boldsymbol{y})$
are of Gevrey class $G^{\delta}(U,L^{\infty}(D))$, that is there exist a
sequence $\boldsymbol{R}=(R_{1},R_{2},\dots)$ with positive components
$R_{i}>0$ and constants $\overline{a},\overline{b},\overline{c}$ independent
of $\boldsymbol{\nu}$ such that

(3.1)
$\begin{split}\left\|{\partial^{\boldsymbol{\nu}}a(\cdot,\boldsymbol{y})}\right\|_{\infty}\leq\frac{\overline{a}}{2}\frac{(\left|\boldsymbol{\nu}\right|!)^{\delta}}{(2\boldsymbol{R})^{\boldsymbol{\nu}}},\quad\left\|{\partial^{\boldsymbol{\nu}}b(\cdot,\boldsymbol{y})}\right\|_{\infty}\leq\frac{\overline{b}}{2}\frac{(\left|\boldsymbol{\nu}\right|!)^{\delta}}{(2\boldsymbol{R})^{\boldsymbol{\nu}}},\quad\left\|{\partial^{\boldsymbol{\nu}}c(\cdot,\boldsymbol{y})}\right\|_{\infty}\leq\frac{\overline{c}}{2}\frac{(\left|\boldsymbol{\nu}\right|!)^{\delta}}{(2\boldsymbol{R})^{\boldsymbol{\nu}}}.\end{split}$

Notice that for $\boldsymbol{\nu}=0$ assumption (3.1) agrees with the upper
bounds in (1.4). Notice also that the components of $\boldsymbol{R}$ are
readily scaled by the factor of $2$. This leads to no loss of generality, but
allows to shorten the forthcoming expressions. For example, in view of (2.10)
assumption (3.1) immediately implies

(3.2)
$\max\bigg{(}\frac{\left\|{\partial^{\boldsymbol{\nu}}a(\cdot,\boldsymbol{y})}\right\|_{\infty}}{\overline{a}},\frac{\left\|{\partial^{\boldsymbol{\nu}}b(\cdot,\boldsymbol{y})}\right\|_{\infty}}{\overline{b}},\frac{\left\|{\partial^{\boldsymbol{\nu}}c(\cdot,\boldsymbol{y})}\right\|_{\infty}}{\overline{c}}\bigg{)}\leq\frac{\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}\right|}}{\boldsymbol{R}^{\boldsymbol{\nu}}(\left|\boldsymbol{\nu}\right|!)^{1-\delta}}$

for all $\boldsymbol{y}\in U$ and $\boldsymbol{\nu}\in\mathcal{F}$.

The main result of this paper is given by the following theorem.

######  Theorem 3.2.

Let $\delta\geq 1$, $\boldsymbol{\nu}\in\mathcal{F}$ and Assumption 3.1 be
valid. There exist constants $C_{\lambda},C_{u}>0$ and a sequence
$\widetilde{\boldsymbol{R}}=(\widetilde{R}_{1},\widetilde{R}_{2},\dots)$ with
positive components $\widetilde{R}_{i}>0$ such that for all $\boldsymbol{y}\in
U$ the derivative of the smallest eigenvalue of (4.1) is bounded by

(3.3)
$\displaystyle\left|\partial^{\boldsymbol{\nu}}\lambda_{1}(\boldsymbol{y})\right|\leq
C_{\lambda}\frac{(\left|\boldsymbol{\nu}\right|!)^{\delta}}{{\widetilde{\boldsymbol{R}}}^{\boldsymbol{\nu}}}$

and the norm of the derivative of the corresponding eigenfunction is bounded
by

(3.4) $\displaystyle\left\|{\partial^{\boldsymbol{\nu}}\nabla
u_{1}(\boldsymbol{y})}\right\|_{L^{2}(D)}\leq
C_{u}\frac{(\left|\boldsymbol{\nu}\right|!)^{\delta}}{{\widetilde{\boldsymbol{R}}}^{\boldsymbol{\nu}}}.$

###### Proof.

Recall that $\left\|{v}\right\|_{V}=\chi_{1}^{-1/2}\left\|{\nabla
v}\right\|_{L^{2}(D)}$. Then the statement follows directly from Theorem 5.1
from Section 5 below by using
$\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}\right|}\leq|\boldsymbol{\nu}|!$
and the constants $C_{\lambda}=\overline{\lambda_{1}}$,
$C_{u}=\chi_{1}^{\frac{1}{2}}\overline{u_{1}}$ and
$\widetilde{R}_{i}=R_{i}/(\sigma\rho)$ with the scaling factors
$\sigma,\rho\geq 1$ from Theorem 5.1. ∎

##  4\. Elliptic eigenvalue problems with countably many parameters

In this section we collect the required notations and facts from the theory of
variational eigenvalue problems. Throughout the paper, when it is unambiguous
we will drop the $\boldsymbol{x}$-dependence when referring to a function
defined on $D$ at a parameter value $\boldsymbol{y}$. For example, we will
write $a(\boldsymbol{y}):=a(\cdot,\boldsymbol{y})$ and analogously for the
other coefficients and eigenfunctions.

For a fixed $\boldsymbol{y}\in U$ the variational formulation of (1.1) reads

(4.1) $\chi_{1}^{-1}\int_{D}a(\boldsymbol{y})\nabla u(\boldsymbol{y})\,\nabla
v+\int_{D}b(\boldsymbol{y})u(\boldsymbol{y})\,v=\lambda(\boldsymbol{y})\int_{D}c(\boldsymbol{y})\,u(\boldsymbol{y})\,v$

where $\chi_{1}$ is the smallest Dirichlet-Laplace eigenvalue defined in (1.3)
being related to the Poincaré constant. For the rescaled norm (1.7) the
Poincaré inequality reads

(4.2) $\left\|{v}\right\|_{L^{2}(D)}\leq\chi_{1}^{-1/2}\left\|{\nabla
v}\right\|_{L^{2}(D)}=\left\|{v}\right\|_{V}.$

For a fixed parameter value $\boldsymbol{y}\in U$ we define the symmetric
bilinear forms ${\mathcal{A}}_{\boldsymbol{y}}:V\times V\rightarrow\mathbb{R}$
and $\mathcal{M}_{\boldsymbol{y}}:L^{2}(D)\times
L^{2}(D)\rightarrow\mathbb{R}$ as

(4.3)
${\mathcal{A}}_{\boldsymbol{y}}(w,v):=\chi_{1}^{-1}\int_{D}a(\boldsymbol{y})\nabla
w(\boldsymbol{y})\,\nabla
v+\int_{D}b(\boldsymbol{y})w\,v,\qquad\mathcal{M}_{\boldsymbol{y}}(w,v):=\int_{D}c(\boldsymbol{y})w\,v.$

In view of the bounds (1.4), ${\mathcal{A}}_{\boldsymbol{y}}$ and
$\mathcal{M}_{\boldsymbol{y}}$ are inner products on $V$ and $L^{2}(D)$,
satisfying

(4.4)
$\begin{split}\underline{a}\|\boldsymbol{v}\|_{V}^{2}\leq{\mathcal{A}}_{\boldsymbol{y}}(\boldsymbol{v},\boldsymbol{v})\leq\frac{\overline{a}+\overline{b}}{2}\|\boldsymbol{v}\|_{V}^{2},&\qquad\boldsymbol{v}\in
V,\\\
\underline{c}\|\boldsymbol{v}\|_{L^{2}(D)}^{2}\leq\mathcal{M}_{\boldsymbol{y}}(\boldsymbol{v},\boldsymbol{v})\leq\frac{\overline{c}}{2}\|\boldsymbol{v}\|_{L^{2}(D)}^{2},&\qquad\boldsymbol{v}\in
L^{2}(D).\end{split}$

Thus, for a fixed $\boldsymbol{y}\in U$, the variational equivalent of (1.1)
in reference domain is the problem of finding a nontrivial
$u(\boldsymbol{y})\in V$ and $\lambda(\boldsymbol{y})\in\mathbb{R}$ such that

(4.5)
$\begin{split}{\mathcal{A}}_{\boldsymbol{y}}(u(\boldsymbol{y}),v)&=\lambda(\boldsymbol{y})\mathcal{M}_{\boldsymbol{y}}(u(\boldsymbol{y}),v),\quad\forall
v\in V,\\\
\left\|{u(\boldsymbol{y})}\right\|_{\mathcal{M}_{\boldsymbol{y}}}&=1,\end{split}$

where the eigenvectors are normalized with respect to the norm
$\left\|{v}\right\|_{\mathcal{M}_{\boldsymbol{y}}}=\sqrt{\mathcal{M}_{\boldsymbol{y}}(v,v)}$.
Similarly as the Dirichlet-Laplacian, for a fixed $\boldsymbol{y}\in U$ the
formulation (4.5) has countably-many eigenvalues
$\\{\lambda_{k}(\boldsymbol{y})\\}_{k\in\mathbb{N}}$, which are positive, have
finite multiplicity and accumulate only at infinity, the smallest eigenvalue
is simple, cf. [11, 8, 16]. Counting multiplicities, we enumerate them in the
non-decreasing order as

$\displaystyle
0<\lambda_{1}(\boldsymbol{y})<\lambda_{2}(\boldsymbol{y})\leq\lambda_{3}(\boldsymbol{y})\leq\dots$

with $\lambda_{k}(\boldsymbol{y})\rightarrow\infty$ as $k\rightarrow\infty$.
For an eigenvalue $\lambda_{k}(\boldsymbol{y})$, we define its eigenspace to
be

$E(\lambda_{k}(\boldsymbol{y})):={\rm{span}}\left\\{u_{k}\,:\,u_{k}\text{ is
an eigenfunction corresponding to }\lambda_{k}(\boldsymbol{y})\right\\}.$

Since $\lambda_{1}$ is simple, $E(\lambda_{1}(\boldsymbol{y}))$ is one-
dimensional and $u_{k}(\boldsymbol{y})$ belong to the orthogonal complement

$E(\lambda_{1}(\boldsymbol{y}))^{\bot}:=\big{\\{}v\leavevmode\nobreak\
:\leavevmode\nobreak\
\mathcal{M}_{\boldsymbol{y}}(v,u_{1}(\boldsymbol{y}))=0\big{\\}}.$

The min-max principle and (4.4) imply uniform upper bounds for the smallest
eigenvalue

(4.6) $\lambda_{1}(\boldsymbol{y})=\inf_{v\in
V}\frac{{\mathcal{A}}_{\boldsymbol{y}}(v,v)}{\mathcal{M}_{\boldsymbol{y}}(v,v)}\leq\frac{\overline{a}+\overline{b}}{2\,\underline{c}}=:\overline{\lambda_{1}}$

and the corresponding eigenfunction

(4.7)
$\left\|{u_{1}(\boldsymbol{y})}\right\|_{V}\leq\sqrt{\frac{\lambda_{1}(\boldsymbol{y})}{\underline{a}}}\leq\sqrt{\frac{\overline{a}+\overline{b}}{2\,\underline{c}\,\underline{a}}}=:\overline{u_{1}}.$

######  Remark 4.1.

The assumption that $b(\boldsymbol{y})\geq 0$ in (1.4) is made without loss of
generality since any EVP with $b<0$, but satisfying the rest of Assumption
3.1, can be shifted to an equivalent problem with a coefficient “new $b$” non-
negative by adding $\sigma\cdot c(\boldsymbol{y})\cdot u(\boldsymbol{y})$ to
both sides of (1.1), where $\sigma$ is chosen such that
$-b(\boldsymbol{y})\leq\sigma\cdot c(\boldsymbol{y})$ for all
$\boldsymbol{x},\boldsymbol{y}$. Such $\sigma$ exists due to Assumption 3.1.
The eigenvalues of the original EVP will become the eigenvalues of the shifted
problem shifted by $-\sigma$, and the corresponding eigenspaces are unchanged.
This has been used previously e.g. in [10, 12].

We denote by $\mu$ the uniform bound of the _relative spectral gap_
$1-\frac{\lambda_{1}(\boldsymbol{y})}{\lambda_{2}(\boldsymbol{y})}$ over all
$\boldsymbol{y}\in U$

(4.8) $\mu:=\min_{\boldsymbol{y}\in
U}\left(1-\frac{\lambda_{1}(\boldsymbol{y})}{\lambda_{2}(\boldsymbol{y})}\right).$

Since $U$ is a bounded domain, the minimum is attained, cf. [7, Proposition
1], and therefore $\mu$ is strictly positive, see also [1] and [12]. Since the
eigenvalues are positive, this implies $0<\mu\leq 1$.

We begin with the following coercive-type estimate, which is required in order
to control the norm of the $\boldsymbol{y}$-derivatives of the eigenfunction
$u_{1}(\boldsymbol{y})$. The proof essentially follows the lines of [12, Lemma
3.1], but covers the case of a parameter-dependent coefficient $c$.

######  Lemma 4.2.

Let $\mu$ be the relative spectral gap defined in (4.8) and the coefficients
$a$, $b$, $c$ satisfy the bounds (1.4). The shifted bilinear form
$\widetilde{\mathcal{A}}_{\boldsymbol{y}}:={\mathcal{A}}_{\boldsymbol{y}}-\lambda_{1}(\boldsymbol{y})\mathcal{M}_{\boldsymbol{y}}$
is uniformly coercive on the orthogonal complement to $u_{1}(\boldsymbol{y})$,
that is

${\mathcal{A}}_{\boldsymbol{y}}(v,v)-\lambda_{1}(\boldsymbol{y})\mathcal{M}_{\boldsymbol{y}}(v,v)\geq
C_{\mu}\left\|{v}\right\|_{V}^{2},\qquad\forall v\in
E(\lambda_{1}(\boldsymbol{y}))^{\bot}$

with the coercivity constant $C_{\mu}=\underline{a}\,\mu$ being independent of
$\boldsymbol{y}$.

###### Proof.

Since the eigenfunctions $\\{u_{k}(\boldsymbol{y})\\}_{k\in\mathbb{N}}$ form a
basis in $V$ that is orthonormal with respect to the inner product
$\mathcal{M}_{\boldsymbol{y}}$, for $v\in
E(\lambda_{1}(\boldsymbol{y}))^{\bot}$, letting
$v_{k}(\boldsymbol{y}):=\mathcal{M}_{\boldsymbol{y}}(v,u_{k}(\boldsymbol{y}))u_{k}(\boldsymbol{y})$
for $k=1,2,\dots,$ we can write

$\displaystyle v(\boldsymbol{y})=\sum_{k=2}^{\infty}v_{k}(\boldsymbol{y})$

Notice that $v_{1}(\boldsymbol{y})=0$ since $v\in
E(\lambda_{1}(\boldsymbol{y}))^{\bot}$. We will suppress the dependence of the
$v$ and $v_{k}$ on $\boldsymbol{y}$ until the end of this proof for brevity.
For $v\in E(\lambda_{1}(\boldsymbol{y}))^{\bot}$, we have

$\displaystyle\widetilde{\mathcal{A}}_{\boldsymbol{y}}(v,v)=\sum_{k=2}^{\infty}\sum_{j=2}^{\infty}\big{(}{\mathcal{A}}_{\boldsymbol{y}}\left(v_{k},v_{j}\right)-\lambda_{1}(\boldsymbol{y})\mathcal{M}_{\boldsymbol{y}}(v_{k},v_{j})\big{)}$

Since all $v_{k}$ are just scaled versions of $u_{k}$, they also satisfy the
variational equation (4.1), so that
${\mathcal{A}}_{\boldsymbol{y}}(v_{k},v_{j})=\lambda_{k}\mathcal{M}_{\boldsymbol{y}}(v_{k},v_{j})$
and they are orthogonal with respect to
$\mathcal{M}_{\boldsymbol{y}}(\cdot,\cdot)$. This implies that
${\mathcal{A}}_{\boldsymbol{y}}(v_{k},v_{j})=0$ for $k\neq j$. Thus, the above
double sum is reduced to

$\displaystyle\widetilde{\mathcal{A}}_{\boldsymbol{y}}(v,v)$
$\displaystyle=\sum_{k=2}^{\infty}\left({\mathcal{A}}_{\boldsymbol{y}}\left(v_{k},v_{k}\right)-\frac{\lambda_{1}(\boldsymbol{y})}{\lambda_{k}(\boldsymbol{y})}{\mathcal{A}}_{\boldsymbol{y}}\left(v_{k},v_{k}\right)\right)$
$\displaystyle\geq\left(1-\frac{\lambda_{1}(\boldsymbol{y})}{\lambda_{2}(\boldsymbol{y})}\right)\sum_{k=2}^{\infty}{\mathcal{A}}_{\boldsymbol{y}}\left(v_{k},v_{k}\right)$
$\displaystyle=\left(1-\frac{\lambda_{1}(\boldsymbol{y})}{\lambda_{2}(\boldsymbol{y})}\right){\mathcal{A}}_{\boldsymbol{y}}\left(v,v\right)\geq\underline{a}\,\mu\left\|{v}\right\|_{V}^{2}=C_{\mu}\left\|{v}\right\|_{V}^{2}.$

This finishes the proof. ∎

##  5\. Parametric regularity: Proof of the main result

The aim of this section is the proof of the main regularity result. As we will
see, the scaling of the constants in the estimates for the
$\boldsymbol{\nu}$-th derivative of the eigenpairs will _solely depend on the
contrast_ of the coefficients of the eigenvalue problem (1.1) defined as
follows, cf. (4.4)

(5.1) $K_{a}:=\frac{\overline{a}+\overline{b}}{2\underline{a}},\qquad
K_{c}:=\frac{\overline{c}}{2\underline{c}}.$

Notice that $K_{a},K_{c}\geq 1$ and the following particular relations of
dimensionless quantities

(5.2)
$\frac{\overline{u_{1}}^{2}\overline{c}}{2}=K_{a}K_{c}=\frac{\overline{\lambda_{1}}\overline{c}}{2\underline{a}},\qquad\frac{\overline{a}+\overline{b}+\overline{\lambda_{1}}\,\overline{c}}{2\underline{a}}=K_{a}(1+K_{c}).$

######  Theorem 5.1.

Let $\delta\geq 1$, $\boldsymbol{\nu}\in\mathcal{F}$, $|\boldsymbol{\nu}|\geq
1$ and $\boldsymbol{R}=(R_{1},R_{2},\dots)$ be defined in Assumption 3.1. Then
there exist $\sigma,\rho\geq 1$ independent of $\boldsymbol{\nu}$ such that
the smallest eigenvalue $\lambda_{1}$ and the corresponding eigenfunction
$u_{1}$ of (4.5) admit the bounds

(5.3)
$\left|\partial^{\boldsymbol{\nu}}\lambda_{1}(\boldsymbol{y})\right|\leq\frac{\overline{\lambda_{1}}}{\rho}\left(\frac{\sigma\rho}{\boldsymbol{R}}\right)^{\boldsymbol{\nu}}\frac{\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}\right|}}{(\left|\boldsymbol{\nu}\right|!)^{1-\delta}}$

and

(5.4)
$\left\|{\partial^{\boldsymbol{\nu}}u_{1}(\boldsymbol{y})}\right\|_{V}\leq\frac{\overline{u_{1}}}{\rho}\left(\frac{\sigma\rho}{\boldsymbol{R}}\right)^{\boldsymbol{\nu}}\frac{\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}\right|}}{(\left|\boldsymbol{\nu}\right|!)^{1-\delta}}.$

The scaling factors $\sigma$, $\rho$ depend solely on the contrast of the
coefficients (5.1) and are defined as

(5.5)
$\begin{array}[]{ll}\sigma:=\mu^{-1}\sigma_{1}+K_{a}K_{c},&\quad\sigma_{1}:=2K_{a}(1+K_{c}),\\\\[8.61108pt]
\rho:=\mu^{-1}\rho_{1}+11K_{a}K_{c},&\quad\rho_{1}:=3\sigma_{1}+16K_{a}K_{c}.\end{array}$

The proof the above Theorem requires upper bounds for the derivatives of the
eigenvalue and the eigenfunction collected in the following Lemma 5.2, Lemma
5.3 and Lemma 5.4.

######  Lemma 5.2.

For sufficiently regular solutions of (4.1) and
$\boldsymbol{\nu}\in\mathcal{F}$ with $|\boldsymbol{\nu}|\geq 1$ there holds

(5.6)
$\begin{split}\left|\partial^{\boldsymbol{\nu}}\lambda_{1}\right|&\leq\bigg{(}\sum_{\boldsymbol{0}<\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\left\|{\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}}\right\|_{V}\left(\left\|{\partial^{\boldsymbol{m}}a}\right\|_{L^{\infty}(D)}+\left\|{\partial^{\boldsymbol{m}}b}\right\|_{L^{\infty}(D)}+\lambda_{1}\left\|{\partial^{\boldsymbol{m}}c}\right\|_{L^{\infty}(D)}\right)\\\
&\quad+\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\left|\partial^{\boldsymbol{\nu}-\boldsymbol{m}}\lambda_{1}\right|\sum_{0\leq\boldsymbol{\ell}\leq\boldsymbol{m}}\left(\boldsymbol{m}\atop\boldsymbol{\ell}\right)\left\|{\partial^{\boldsymbol{m}-\boldsymbol{\ell}}u_{1}}\right\|_{V}\left\|{\partial^{\boldsymbol{\ell}}c}\right\|_{L^{\infty}(D)}\bigg{)}\left\|{u_{1}}\right\|_{V}.\end{split}$

###### Proof.

The Leibniz general product rule for the $\boldsymbol{\nu}$-th derivative of
(4.1) with respect to $\boldsymbol{y}$ implies

(5.7) $\begin{split}\int_{D}&\chi_{1}^{-1}\partial^{\boldsymbol{\nu}}(a\nabla
u_{1})\cdot\nabla
v+\partial^{\boldsymbol{\nu}}(bu_{1})v-\lambda_{1}\partial^{\boldsymbol{\nu}}(cu_{1})v=\sum_{\boldsymbol{0}\leq\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\partial^{\boldsymbol{\nu}-\boldsymbol{m}}\lambda_{1}\int_{D}\,\partial^{\boldsymbol{m}}(cu_{1})v,\end{split}$

that is valid for all $v\in V$. Notice that the term with
$\boldsymbol{m}=\boldsymbol{\nu}$ is moved to the left-hand side. To obtain a
bound on the derivatives of the eigenvalue, we specify $v=u_{1}$ and split
$\partial^{\boldsymbol{\nu}}\lambda_{1}$ from the remaining terms. This
implies

$\begin{split}\int_{D}&\chi_{1}^{-1}\partial^{\boldsymbol{\nu}}(a\nabla
u_{1})\cdot\nabla
u_{1}+\partial^{\boldsymbol{\nu}}(bu_{1})u_{1}-\lambda_{1}\partial^{\boldsymbol{\nu}}(cu_{1})u_{1}=\partial^{\boldsymbol{\nu}}\lambda_{1}+\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\partial^{\boldsymbol{\nu}-\boldsymbol{m}}\lambda_{1}\int_{D}\,\partial^{\boldsymbol{m}}(cu_{1})u_{1},\end{split}$

where we have used the normalization
$\mathcal{M}_{\boldsymbol{y}}(u_{1},u_{1})=1$ in the right-hand side. Next, we
apply the Leibniz general product rule to expand the derivatives and arrive at

$\begin{split}\partial^{\boldsymbol{\nu}}\lambda_{1}=&\sum_{\boldsymbol{0}<\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\int_{D}\bigg{(}\chi_{1}^{-1}\partial^{\boldsymbol{m}}a\,\partial^{\boldsymbol{\nu}-\boldsymbol{m}}\nabla
u_{1}\cdot\nabla
u_{1}+\partial^{\boldsymbol{m}}b\,\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}u_{1}-\lambda_{1}\partial^{\boldsymbol{m}}c\,\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}u_{1}\bigg{)}\\\
&-\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\partial^{\boldsymbol{\nu}-\boldsymbol{m}}\lambda_{1}\sum_{0\leq\boldsymbol{\ell}\leq\boldsymbol{m}}\left(\boldsymbol{m}\atop\boldsymbol{\ell}\right)\int_{D}\,\partial^{\boldsymbol{\ell}}c\,\partial^{\boldsymbol{m}-\boldsymbol{\ell}}u_{1}u_{1}.\end{split}$

Observe that the terms with $\boldsymbol{m}=\boldsymbol{0}$ vanish in the
first sum. Indeed, if $\partial^{\boldsymbol{\nu}}u_{1}\in V$, equation (4.1)
implies
${\mathcal{A}}_{\boldsymbol{y}}(u_{1},\partial^{\boldsymbol{\nu}}u_{1})-\lambda_{1}\mathcal{M}_{\boldsymbol{y}}(u_{1},\partial^{\boldsymbol{\nu}}u_{1})=0$.
The statement of the Lemma follows now by the triangle inequality, the Cauchy-
Schwarz inequality and the bound
$\left\|{u_{1}}\right\|_{L^{2}(D)}\leq\left\|{u_{1}}\right\|_{V}$, cf. (4.2).
∎

To obtain the bound on the derivatives of the eigenfunction, consider the
unique bi-orthogonal decomposition

(5.8) $\partial^{\boldsymbol{\nu}}u_{1}=p_{1}+p^{\bot},\qquad p_{1}\in
E(\lambda_{1}(\boldsymbol{y})),\quad p^{\bot}\in
E(\lambda_{1}(\boldsymbol{y}))^{\bot}$

that is valid for every fixed $\boldsymbol{y}\in U$. The bi-orthogonality
means that the relations ${\mathcal{A}}_{\boldsymbol{y}}(p_{1},p^{\bot})=0$
and $\mathcal{M}_{\boldsymbol{y}}(p_{1},p^{\bot})=0$ are simultaneously
satisfied. The triangle inequality reveals

(5.9)
$\|\partial^{\boldsymbol{\nu}}u_{1}\|_{V}\leq\|p_{1}\|_{V}+\|p^{\bot}\|_{V}.$

In the following we will estimate the norms of the projections $p_{1}$ and
$p^{\bot}$ independently. The following Lemma shows that the norm of
$p^{\bot}$ admits a similar upper bound as the eigenvalue in Lemma 5.2 that
involves lower order derivatives of $\lambda_{1}$ and $u_{1}$ only.

######  Lemma 5.3.

For sufficiently regular solutions of (4.1) and
$\boldsymbol{\nu}\in\mathcal{F}$ with $|\boldsymbol{\nu}|\geq 1$ there holds

(5.10)
$\begin{split}C_{\mu}\|p^{\bot}\|_{V}&\leq\sum_{\boldsymbol{0}<\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\left\|{\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}}\right\|_{V}\left(\left\|{\partial^{\boldsymbol{m}}a}\right\|_{\infty}+\left\|{\partial^{\boldsymbol{m}}b}\right\|_{\infty}+\lambda_{1}\left\|{\partial^{\boldsymbol{m}}c}\right\|_{\infty}\right)\\\
&\quad+\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\left|\partial^{\boldsymbol{\nu}-\boldsymbol{m}}\lambda_{1}\right|\sum_{0\leq\boldsymbol{\ell}\leq\boldsymbol{m}}\left(\boldsymbol{m}\atop\boldsymbol{\ell}\right)\left\|{\partial^{\boldsymbol{m}-\boldsymbol{\ell}}u_{1}}\right\|_{V}\left\|{\partial^{\boldsymbol{\ell}}c}\right\|_{\infty}.\end{split}$

###### Proof.

Recall identity (5.7) that is valid for all $v\in V$. In the particular case
$v=p^{\bot}$ it reads

$\begin{split}\int_{D}&\chi_{1}^{-1}\partial^{\boldsymbol{\nu}}(a\nabla
u_{1})\cdot\nabla
p^{\bot}+\partial^{\boldsymbol{\nu}}(bu_{1})p^{\bot}-\lambda_{1}\partial^{\boldsymbol{\nu}}(cu_{1})p^{\bot}=\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\partial^{\boldsymbol{\nu}-\boldsymbol{m}}\lambda_{1}\int_{D}\,\partial^{\boldsymbol{m}}(cu_{1})p^{\bot}.\end{split}$

Notice that the term with $\boldsymbol{m}=\boldsymbol{0}$ on the right-hand
side has vanished, since $u_{1}$ and $p^{\bot}$ are
$\mathcal{M}_{\boldsymbol{y}}$-orthogonal for all fixed $\boldsymbol{y}\in U$.
Therefore the highest order derivative
$\partial^{\boldsymbol{\nu}}\lambda_{1}$ is not present in this equation.
Next, we apply the Leibniz general product rule and obtain

(5.11)
$\begin{split}\sum_{\boldsymbol{0}\leq\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\int_{D}\chi_{1}^{-1}\partial^{\boldsymbol{m}}a\,\partial^{\boldsymbol{\nu}-\boldsymbol{m}}(\nabla
u_{1})\cdot\nabla
p^{\bot}+\partial^{\boldsymbol{m}}b\,\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}\,p^{\bot}-\lambda_{1}\partial^{\boldsymbol{m}}c\,\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}\,p^{\bot}&\\\
=\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\partial^{\boldsymbol{\nu}-\boldsymbol{m}}\lambda_{1}\sum_{\boldsymbol{0}\leq\boldsymbol{\ell}\leq\boldsymbol{m}}\left(\boldsymbol{m}\atop\boldsymbol{\ell}\right)\int_{D}\,\partial^{\boldsymbol{\ell}}c\,\partial^{\boldsymbol{m}-\boldsymbol{\ell}}u_{1}\,p^{\bot}&.\end{split}$

Observe that the highest order derivatives $\partial^{\boldsymbol{\nu}}u_{1}$
do only appear in the left-hand side in the terms with
$\boldsymbol{m}=\boldsymbol{0}$. For this terms there holds the lower bound

(5.12)
$\begin{split}\int_{D}\chi_{1}^{-1}a\,\partial^{\boldsymbol{\nu}}(\nabla
u_{1})\cdot\nabla
p^{\bot}+b\,\partial^{\boldsymbol{\nu}}u_{1}\,p^{\bot}-\lambda_{1}c\,\partial^{\boldsymbol{\nu}}u_{1}\,p^{\bot}&={\mathcal{A}}_{\boldsymbol{y}}(\partial^{\boldsymbol{\nu}}u_{1},p^{\bot})-\lambda_{1}\mathcal{M}_{\boldsymbol{y}}(\partial^{\boldsymbol{\nu}}u_{1},p^{\bot})\\\
&\geq C_{\mu}\|p^{\bot}\|_{V}^{2}.\end{split}$

To see this we recall the definition of ${\mathcal{A}}_{\boldsymbol{y}}$ and
$\mathcal{M}_{\boldsymbol{y}}$, the bi-orthogonal decomposition (5.8) and
Lemma 4.2. The statement of the Lemma follows now from (5.11) and (5.12) by
the triangle inequality, the Cauchy-Schwarz inequality and the Poincaré
inequality (4.2). ∎

######  Lemma 5.4.

For sufficiently regular solutions of (4.1) and
$\boldsymbol{\nu}\in\mathcal{F}$ with $|\boldsymbol{\nu}|\geq 1$ there holds

(5.13)
$\begin{split}\left\|{p_{1}}\right\|_{V}&\leq\frac{\|u_{1}\|_{V}^{2}}{2}\sum_{\boldsymbol{0}<\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\left\|{\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}}\right\|_{V}\left\|{\partial^{\boldsymbol{m}}c}\right\|_{\infty}\\\
&\quad+\frac{\|u_{1}\|_{V}}{2}\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\left\|{\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}}\right\|_{V}\sum_{\boldsymbol{0}\leq\boldsymbol{\ell}\leq\boldsymbol{m}}\left(\boldsymbol{m}\atop\boldsymbol{\ell}\right)\left\|{\partial^{\boldsymbol{m}-\boldsymbol{\ell}}u_{1}}\right\|_{V}\left\|{\partial^{\boldsymbol{\ell}}c}\right\|_{\infty}.\end{split}$

###### Proof.

Since $p_{1}$ is the projection of $\partial^{\boldsymbol{\nu}}u_{1}$ onto
$E(\lambda_{1}(\boldsymbol{y}))$, we have
$p_{1}=\mathcal{M}_{\boldsymbol{y}}(\partial^{\boldsymbol{\nu}}u_{1},u_{1})\,u_{1}$
and hence

(5.14)
$\left\|{p_{1}}\right\|_{V}\leq\left|\mathcal{M}_{\boldsymbol{y}}(\partial^{\boldsymbol{\nu}}u_{1},u_{1})\right|\,\|u_{1}\|_{V}.$

It remains to estimate the modulus of the coefficient. For this we recall that
$\mathcal{M}_{\boldsymbol{y}}(u_{1},u_{1})=1$. Thus, taking the
$\boldsymbol{\nu}$-th derivatives of this identity we get by the Leibniz
general product rule

$0=\partial^{\boldsymbol{\nu}}\mathcal{M}_{\boldsymbol{y}}(u_{1},u_{1})=\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\int_{D}\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}\,\partial^{\boldsymbol{m}}(c\,u_{1})+\int_{D}\,(\partial^{\boldsymbol{\nu}}u_{1})\,c\,u_{1}+\int_{D}u_{1}\,\partial^{\boldsymbol{\nu}}(c\,u_{1}).$

The highest order derivative $\partial^{\boldsymbol{\nu}}u_{1}$ does only
appear in the last two integrals. Splitting it from the remaining terms
implies

(5.15)
$2\mathcal{M}_{\boldsymbol{y}}(\partial^{\boldsymbol{\nu}}u_{1},u_{1})=-\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\int_{D}\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}\,\partial^{\boldsymbol{m}}(c\,u_{1})-\sum_{\boldsymbol{0}<\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\int_{D}u_{1}\,\partial^{\boldsymbol{m}}c\,\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}$

since $\mathcal{M}_{\boldsymbol{y}}$ is symmetric. The statement of the Lemma
follows now from (5.14) and (5.15) by the triangle inequality, the Cauchy-
Schwarz inequality and the Poincaré inequality (4.2). ∎

Inequality (5.9), Lemma 5.3 and Lemma 5.4 imply the estimate

$\displaystyle C_{\mu}\|\partial^{\boldsymbol{\nu}}u_{1}\|_{V}$
$\displaystyle\leq\sum_{\boldsymbol{0}<\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\left\|{\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}}\right\|_{V}\left(\left\|{\partial^{\boldsymbol{m}}a}\right\|_{\infty}+\left\|{\partial^{\boldsymbol{m}}b}\right\|_{\infty}+(\lambda_{1}+\tfrac{C_{\mu}}{2}\|u_{1}\|_{V}^{2})\left\|{\partial^{\boldsymbol{m}}c}\right\|_{\infty}\right)$
(5.16)
$\displaystyle+\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)(\left|\partial^{\boldsymbol{\nu}-\boldsymbol{m}}\lambda_{1}\right|+\tfrac{C_{\mu}}{2}\|u_{1}\|_{V}\|\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}\|_{V})\sum_{0\leq\boldsymbol{\ell}\leq\boldsymbol{m}}\left(\boldsymbol{m}\atop\boldsymbol{\ell}\right)\left\|{\partial^{\boldsymbol{m}-\boldsymbol{\ell}}u_{1}}\right\|_{V}\left\|{\partial^{\boldsymbol{\ell}}c}\right\|_{\infty}.$

###### Proof of Theorem 5.1.

We now prove the main theorem by induction. First, we establish the base for
first derivative of the eigenvalue. Letting $\boldsymbol{\nu}$ be the unit
multiindex, that is $\boldsymbol{\nu}=\boldsymbol{e}$ with
$|\boldsymbol{e}|=1$ we utilize (5.6) together with (4.6), (4.7), (3.2) and
the identities (5.2) to obtain

$\displaystyle\left|\partial^{\boldsymbol{e}}\lambda_{1}\right|\leq\overline{u_{1}}^{2}\left(\frac{\overline{a}\left[\tfrac{1}{2}\right]_{1}}{\boldsymbol{R}^{\boldsymbol{e}}}+\frac{\overline{b}\left[\tfrac{1}{2}\right]_{1}}{\boldsymbol{R}^{\boldsymbol{e}}}+\frac{\overline{\lambda_{1}}\,\overline{c}\left[\tfrac{1}{2}\right]_{1}}{\boldsymbol{R}^{\boldsymbol{e}}}\right)=\frac{\overline{\lambda_{1}}}{\underline{a}}\left(\overline{a}+\overline{b}+\overline{\lambda_{1}}\,\overline{c}\right)\frac{\left[\tfrac{1}{2}\right]_{1}}{\boldsymbol{R}^{\boldsymbol{e}}}=\overline{\lambda_{1}}\frac{\sigma_{1}\left[\tfrac{1}{2}\right]_{1}}{\boldsymbol{R}^{\boldsymbol{e}}}$

For the eigenfunction we substitute $\boldsymbol{\nu}=\boldsymbol{e}$ into
(5.16) and use (4.6), (4.7), (3.2) to get

$\displaystyle\left\|{\partial^{\boldsymbol{e}}u_{1}}\right\|_{V}$
$\displaystyle\leq\frac{\overline{u_{1}}}{C_{\mu}}\left(\frac{\overline{a}\left[\tfrac{1}{2}\right]_{1}}{\boldsymbol{R}^{\boldsymbol{e}}}+\frac{\overline{b}\left[\tfrac{1}{2}\right]_{1}}{\boldsymbol{R}^{\boldsymbol{e}}}+\left(\overline{\lambda_{1}}+\frac{C_{\mu}\overline{u_{1}}^{2}}{2}\right)\frac{\overline{c}\left[\tfrac{1}{2}\right]_{1}}{\boldsymbol{R}^{\boldsymbol{e}}}\right)=\overline{u_{1}}\left(\frac{\overline{a}+\overline{b}+\overline{\lambda_{1}}\,\overline{c}}{\underline{a}\,\mu}+\frac{\overline{u_{1}}^{2}\overline{c}}{2}\right)\frac{\left[\tfrac{1}{2}\right]_{1}}{\boldsymbol{R}^{\boldsymbol{e}}}.$

By (5.2) the term in parentheses equals to $\sigma$. This concludes the base
of induction. For the inductive step for the eigenvalue derivative bound,
suppose that $\left|\boldsymbol{\nu}\right|\geq 2$ and the bounds (5.3) and
(5.4) hold for all multi-indicies of order strictly less than
$\boldsymbol{\nu}$. In particular, inductive assumptions imply

(5.17)
$\displaystyle\max\left(\frac{\left|\partial^{\boldsymbol{\nu}-\boldsymbol{m}}\lambda_{1}(\boldsymbol{y})\right|}{\overline{\lambda_{1}}},\frac{\left\|{\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}}\right\|_{V}}{\overline{u_{1}}}\right)\leq\frac{(\sigma\rho)^{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}}{\rho\boldsymbol{R}^{{\boldsymbol{\nu}-\boldsymbol{m}}}(\left|\boldsymbol{\nu}-\boldsymbol{m}\right|!)^{1-\delta}},\qquad\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu},\quad|\boldsymbol{\nu}|\geq
2.$

Since $\sigma,\rho\geq 1$ we have for $\boldsymbol{m}>\boldsymbol{0}$ the
trivial upper bound
$(\sigma\rho)^{|\boldsymbol{\nu}-\boldsymbol{m}|}\leq\rho^{-1}(\sigma\rho)^{|\boldsymbol{\nu}|}$,
and thereby

(5.18)
$\displaystyle\left\|{\partial^{\boldsymbol{\nu}-\boldsymbol{m}}u_{1}}\right\|_{V}\leq\frac{\overline{u_{1}}(\sigma\rho)^{\left|\boldsymbol{\nu}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}}{\rho^{2}\boldsymbol{R}^{{\boldsymbol{\nu}-\boldsymbol{m}}}(\left|\boldsymbol{\nu}-\boldsymbol{m}\right|!)^{1-\delta}},\qquad\boldsymbol{0}<\boldsymbol{m}\leq\boldsymbol{\nu},\quad|\boldsymbol{\nu}|\geq
2.$

Notice that this estimate indeed holds for the extended range including
$\boldsymbol{m}=\boldsymbol{\nu}$, since in this case, according to (4.7), we
have
$\|u_{1}\|_{V}\leq\overline{u_{1}}\leq\overline{u_{1}}\rho^{-2}(\sigma\rho)^{|\boldsymbol{\nu}|}$.
The following estimate also follows with similar arguments

(5.19)
$\displaystyle\left\|{\partial^{\boldsymbol{m}-\boldsymbol{\ell}}u_{1}}\right\|_{V}\leq\frac{\overline{u_{1}}(\sigma\rho)^{\left|\boldsymbol{m}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{m}-\boldsymbol{\ell}\right|}}{\rho\boldsymbol{R}^{{\boldsymbol{m}-\boldsymbol{\ell}}}(\left|\boldsymbol{m}-\boldsymbol{\ell}\right|!)^{1-\delta}},\qquad\boldsymbol{0}\leq\boldsymbol{\ell}\leq\boldsymbol{m},\quad\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu},\quad|\boldsymbol{\nu}|\geq
2.$

We are now ready to proceed with the inductive step. For the derivative of the
eigenvalue recall (5.6), the regularity assumptions (3.2) and the bounds
(5.17) – (5.19) to obtain

$\displaystyle\left|\partial^{\boldsymbol{\nu}}\lambda_{1}\right|$
$\displaystyle\leq\overline{u_{1}}\sum_{\boldsymbol{0}<\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\frac{\overline{u_{1}}(\sigma\rho)^{\left|\boldsymbol{\nu}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}}{\rho^{2}\boldsymbol{R}^{{\boldsymbol{\nu}-\boldsymbol{m}}}(\left|\boldsymbol{\nu}-\boldsymbol{m}\right|!)^{1-\delta}}\left(\frac{(\overline{a}+\overline{b})\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{m}\right|}}{\boldsymbol{R}^{\boldsymbol{m}}(\left|\boldsymbol{m}\right|!)^{1-\delta}}+\overline{\lambda_{1}}\frac{\overline{c}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{m}\right|}}{\boldsymbol{R}^{\boldsymbol{m}}(\left|\boldsymbol{m}\right|!)^{\delta-1}}\right)$
$\displaystyle+\overline{u_{1}}\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\frac{\overline{\lambda_{1}}(\sigma\rho)^{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}}{\rho\boldsymbol{R}^{{\boldsymbol{\nu}-\boldsymbol{m}}}(\left|\boldsymbol{\nu}-\boldsymbol{m}\right|!)^{1-\delta}}\sum_{0\leq\boldsymbol{\ell}\leq\boldsymbol{m}}\left(\boldsymbol{m}\atop\boldsymbol{\ell}\right)\frac{\overline{u_{1}}(\sigma\rho)^{\left|\boldsymbol{m}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{m}-\boldsymbol{\ell}\right|}}{\rho\boldsymbol{R}^{{\boldsymbol{m}-\boldsymbol{\ell}}}(\left|\boldsymbol{m}-\boldsymbol{\ell}\right|!)^{1-\delta}}\frac{\overline{c}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\ell}\right|}}{\boldsymbol{R}^{\boldsymbol{\ell}}(\left|\boldsymbol{\ell}\right|!)^{1-\delta}},$

Since $\delta\geq 1$ and
$|\boldsymbol{\nu}-\boldsymbol{m}|!|\boldsymbol{m}|!\leq|\boldsymbol{\nu}|!$
and likewise
$|\boldsymbol{\nu}-\boldsymbol{m}|!|\boldsymbol{m}-\boldsymbol{\ell}|!|\boldsymbol{\ell}|!\leq|\boldsymbol{\nu}|!$
the term $(|\boldsymbol{\nu}|!)^{\delta-1}$ can be factored out by making the
right-hand side larger. Collecting the remaining terms we obtain

$\displaystyle\left|\partial^{\boldsymbol{\nu}}\lambda_{1}\right|$
$\displaystyle\leq\overline{u_{1}}^{2}(\overline{a}+\overline{b}+\overline{\lambda_{1}}\overline{c})\frac{(\sigma\rho)^{\left|\boldsymbol{\nu}\right|}}{\rho^{2}\boldsymbol{R}^{\boldsymbol{\nu}}}(\left|\boldsymbol{\nu}\right|!)^{\delta-1}\sum_{\boldsymbol{0}<\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{m}\right|}$
$\displaystyle+\overline{\lambda_{1}}\overline{u_{1}}^{2}\,\overline{c}\frac{(\sigma\rho)^{\left|\boldsymbol{\nu}\right|}}{\rho^{2}\boldsymbol{R}^{\boldsymbol{\nu}}}(\left|\boldsymbol{\nu}\right|!)^{\delta-1}\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\sum_{0\leq\boldsymbol{\ell}\leq\boldsymbol{m}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\left(\boldsymbol{m}\atop\boldsymbol{\ell}\right)\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{m}-\boldsymbol{\ell}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\ell}\right|}.$

Thank to the bounds (6.9), (6.10) from the appendix, the first sum equals to
$3[\tfrac{1}{2}]_{|\boldsymbol{\nu}|}$, the second is bounded by
$8[\tfrac{1}{2}]_{|\boldsymbol{\nu}|}$. Notice moreover that
$\overline{u_{1}}^{2}=\underline{a}^{-1}\overline{\lambda_{1}}$ as we arrive
at

(5.20)
$\displaystyle\left|\partial^{\boldsymbol{\nu}}\lambda_{1}\right|\leq\overline{\lambda_{1}}\left(\frac{3(\overline{a}+\overline{b}+\overline{\lambda_{1}}\overline{c})}{\underline{a}}+8\overline{u_{1}}^{2}\overline{c}\right)\frac{(\sigma\rho)^{\left|\boldsymbol{\nu}\right|}}{\rho^{2}\boldsymbol{R}^{\boldsymbol{\nu}}}(\left|\boldsymbol{\nu}\right|!)^{\delta-1}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}\right|}.$

According to the definition (5.1) of the coefficient contrast $K_{a}$ and
$K_{c}$ the term in parentheses equals $\rho_{1}$. Since $\rho_{1}\leq\rho$
this implies the statement (5.3) for the derivatives of the eigenvalue.

The inductive step for the eigenfunction is similar. Here we combine (5.16)
with the regularity assumptions (3.2) and the bounds (5.17) – (5.19) to obtain

$\displaystyle\|\partial^{\boldsymbol{\nu}}u_{1}\|_{V}$ $\displaystyle\leq
C_{\mu}^{-1}\sum_{\boldsymbol{0}<\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\frac{\overline{u_{1}}(\sigma\rho)^{\left|\boldsymbol{\nu}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}}{\rho^{2}\boldsymbol{R}^{{\boldsymbol{\nu}-\boldsymbol{m}}}(\left|\boldsymbol{\nu}-\boldsymbol{m}\right|!)^{1-\delta}}\frac{(\overline{a}+\overline{b}+\left(\overline{\lambda_{1}}+\frac{1}{2}C_{\mu}\overline{u_{1}}^{2}\right)\overline{c})\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{m}\right|}}{\boldsymbol{R}^{\boldsymbol{m}}(\left|\boldsymbol{m}\right|!)^{1-\delta}}$
$\displaystyle+C_{\mu}^{-1}\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\frac{\overline{\lambda_{1}}(\sigma\rho)^{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}}{\rho\boldsymbol{R}^{{\boldsymbol{\nu}-\boldsymbol{m}}}(\left|\boldsymbol{\nu}-\boldsymbol{m}\right|!)^{1-\delta}}\sum_{0\leq\boldsymbol{\ell}\leq\boldsymbol{m}}\left(\boldsymbol{m}\atop\boldsymbol{\ell}\right)\frac{\overline{u_{1}}(\sigma\rho)^{\left|\boldsymbol{m}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{m}-\boldsymbol{\ell}\right|}}{\rho\boldsymbol{R}^{{\boldsymbol{m}-\boldsymbol{\ell}}}(\left|\boldsymbol{m}-\boldsymbol{\ell}\right|!)^{1-\delta}}\frac{\overline{c}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\ell}\right|}}{\boldsymbol{R}^{\boldsymbol{\ell}}(\left|\boldsymbol{\ell}\right|!)^{1-\delta}}$
$\displaystyle+\frac{\overline{u_{1}}}{2}\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\frac{\overline{u_{1}}(\sigma\rho)^{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}}{\rho\boldsymbol{R}^{{\boldsymbol{\nu}-\boldsymbol{m}}}(\left|\boldsymbol{\nu}-\boldsymbol{m}\right|!)^{1-\delta}}\sum_{0\leq\boldsymbol{\ell}\leq\boldsymbol{m}}\left(\boldsymbol{m}\atop\boldsymbol{\ell}\right)\frac{\overline{u_{1}}(\sigma\rho)^{\left|\boldsymbol{m}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{m}-\boldsymbol{\ell}\right|}}{\rho\boldsymbol{R}^{{\boldsymbol{m}-\boldsymbol{\ell}}}(\left|\boldsymbol{m}-\boldsymbol{\ell}\right|!)^{1-\delta}}\frac{\overline{c}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\ell}\right|}}{\boldsymbol{R}^{\boldsymbol{\ell}}(\left|\boldsymbol{\ell}\right|!)^{1-\delta}}.$

Extracting the term with the factorial and recalling that
$C_{\mu}=\underline{a}\,\mu$ we find

$\begin{split}\left\|{\partial^{\boldsymbol{\nu}}u_{1}}\right\|_{V}&\leq\overline{u_{1}}\left(\frac{(\overline{a}+\overline{b}+\overline{\lambda_{1}}\,\overline{c})}{\underline{a}\,\mu}+\frac{\overline{u_{1}}^{2}\overline{c}}{2}\right)\frac{(\sigma\rho)^{\left|\boldsymbol{\nu}\right|}}{\rho^{2}\boldsymbol{R}^{\boldsymbol{\nu}}}(\left|\boldsymbol{\nu}\right|!)^{\delta-1}\sum_{\boldsymbol{0}<\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{m}\right|}\\\
&\quad+\overline{u_{1}}\left(\frac{\overline{\lambda_{1}}\overline{c}}{\underline{a}\,\mu}+\frac{\overline{u_{1}}^{2}\overline{c}}{2}\right)\frac{(\sigma\rho)^{\left|\boldsymbol{\nu}\right|}}{\rho^{2}\boldsymbol{R}^{\boldsymbol{\nu}}}(\left|\boldsymbol{\nu}\right|!)^{\delta-1}\sum_{\boldsymbol{0}<\boldsymbol{m}<\boldsymbol{\nu}}\sum_{0\leq\boldsymbol{\ell}\leq\boldsymbol{m}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\left(\boldsymbol{m}\atop\boldsymbol{\ell}\right)\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}-\boldsymbol{m}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{m}-\boldsymbol{\ell}\right|}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\ell}\right|}\end{split}$

By Lemma 6.2 the sums above can be estimated by
$3[\tfrac{1}{2}]_{|\boldsymbol{\nu}|}$ and
$8[\tfrac{1}{2}]_{|\boldsymbol{\nu}|}$, respectively, and hence

$\begin{split}\left\|{\partial^{\boldsymbol{\nu}}u_{1}}\right\|_{V}\leq\overline{u_{1}}\left(\frac{3(\overline{a}+\overline{b})+11\overline{\lambda_{1}}\,\overline{c}}{\underline{a}\,\mu}+11\frac{\overline{u_{1}}^{2}\,\overline{c}}{2}\right)\frac{(\sigma\rho)^{\left|\boldsymbol{\nu}\right|}}{\rho^{2}\boldsymbol{R}^{\boldsymbol{\nu}}}(\left|\boldsymbol{\nu}\right|!)^{\delta-1}\left[\tfrac{1}{2}\right]_{\left|\boldsymbol{\nu}\right|}\end{split}$

The definition (5.1) of the coefficient contrast $K_{a}$ and $K_{c}$ implies
that the term in the parentheses coincides with $\rho_{2}$. Since
$\rho_{2}\leq\rho$, the statement (5.4) follows. This finishes the proof. ∎

##  6\. Numerical Experiments

Let us mention that several existing publications contain numerical studies on
parametric and stochastic eigenvalue problems, see e.g. [12, 14]. These
numerical experiments involve advanced numerical integration approaches like
the Quasi-Monte Carlo and a multilevel version thereof. The sharper regularity
result from our paper (for the analytic case $\delta=1$) has immediate
applications in the context of [12, 14]. The aim of this section is however
not to reproduce the existing tests using advanced numerical techniques, but
rather give an elementary example that demonstrates the propagation of the
analytic and Gevrey-$\delta$ regularity from the diffusion coefficient to the
smallest eigenvalue that preserves the same value $\delta$.

Let the domain $D=(0,1)^{2}$ and consider the problem (1.1) with $b\equiv 0$
and $c\equiv 1$ and $a$ is one of following functions

(6.1) $a^{(1)}(\boldsymbol{x},y)=2+\sin\left(\pi(x_{1}+x_{2}+y)\right)$

and

(6.2)
$a^{(2)}(\boldsymbol{x},y)=1+\left(x_{1}+x_{2}\right)\exp\left(-\frac{1}{\sqrt{y+1}}\right).$

Here $y$ is a scalar real random variable uniformly distributed in $[-1,1]$.
The aim is to compute the expected value of the smallest eigenvalue

(6.3) $\mathbb{E}[\lambda_{1}]=\frac{1}{2}\int_{-1}^{1}\lambda_{1}(y)\,dy.$

Since $\lambda_{1}$ is not available in closed form the integral (6.3) can be
approximated by numerical quadrature, e.g. the Gauss-Legendre quadrature which
will be efficient for the case of a single real-valued parameter. Let
$\\{\xi_{i},w_{i}\\}_{i=1}^{n}$ be the nodes and weights of the $n$-point
Gauss-Legendre quadrature

(6.4) $Q_{n}[\lambda_{1}]=\frac{1}{2}\sum_{i=1}^{n}w_{i}\lambda_{1}(\xi_{i}).$

We are interested in the behaviour of the quadrature error

(6.5) $\varepsilon_{n}=|\mathbb{E}[\lambda_{1}]-Q_{n}[\lambda_{1}]|$

with increasing $n$. It is known that the convergence of $\varepsilon_{n}$ is
strongly related to the regularity of $\lambda_{1}$ with respect to $y$. In
particular [4, Theorem 5.2] implies that

(6.6) $\varepsilon_{n}\leq C\exp(-rn^{1/\delta}).$

with positive constants $C$ and $r$ independent of $n$.

  * (1)

In the case of the analytic diffusion coefficient $a^{(1)}$ as in (6.1)
Theorem 3.2 implies that $\lambda_{1}$ is analytic in $y$, i.e. $\delta=1$,
and therefore we expect

(6.7) $\varepsilon_{n}^{(1)}\leq C\exp(-rn).$

  * (2)

The diffusion coefficient $a^{(2)}$ is not analytic near $y=-1$, but is
Gevrey-$\delta$ uniformly for all $y\in[-1,1]$ with $\delta\geq 3$, see [3]
and [4, Section 6]. Theorem 3.2 implies that $\lambda_{1}$ is Gevrey-$\delta$
with the same $\delta=3$ and hence we expect

(6.8) $\varepsilon_{n}^{(2)}\leq C\exp(-rn^{1/3}).$

In order to observe the behaviour predicted in (6.7) and (6.8) we solve
deterministic equations (4.1) in every quadrature point $y=\xi_{i}$ on a very
fine finite element grid having $16.129$ degrees of freedom. Since
$\mathbb{E}[\lambda_{1}]$ is not available in closed form, we approximate it
by the very fine Gauss-Legendre quadrature $Q_{n^{*}}[\lambda_{1}]$ with
$n^{*}=40$ quadrature nodes for $a^{(1)}$ and $n^{*}=123$ quadrature nodes for
$a^{(2)}$. As an eigensolver we use the RQ algorithm with an absolute error
tolerance of $10^{-14}$.

In Figure 1, we plot the relative error $\varepsilon_{n}^{(1)}$ against the
number of quadrature points $n$ in the semi-logarithmic scale. The reference
line clearly shows the linear trend of the type $-rn+\log C$ and thereby
confirms (6.7).

In Figure 2, we plot the relative error $\varepsilon_{n}^{(2)}$ with respect
to the third root of the number of quadrature points $m:=n^{1/3}$ in the semi-
logarithmic scale. Here we can also observe the linear trend of the type
$-rm+\log C$. This confirms (6.8) and thereby demonstrates the meaning and
validity of Theorem 3.2.

Figure 1. Quadrature error $\varepsilon^{(1)}_{n}$ with respect to the number
$n$ of quadrature points.

Figure 2. Quadrature error $\varepsilon^{(2)}_{n}$ with respect to
$m=n^{1/3}$, where $n$ is the number of quadrature points.

## References

  * [1] Roman Andreev and Christoph Schwab, _Sparse tensor approximation of parametric eigenvalue problems_ , Lecture Notes in Computational Science and Engineering, Springer Berlin Heidelberg, August 2011, pp. 203–241.
  * [2] Louis Boutet de Monvel and Paul Krée, _Pseudo-differential operators and Gevrey classes_ , Ann. Inst. Fourier (Grenoble) 17 (1967), no. fasc. 1, 295–323. MR 226170
  * [3] Hua Chen and Luigi Rodino, _General theory of PDE and Gevrey classes_ , General theory of partial differential equations and microlocal analysis (Trieste, 1995), Pitman Res. Notes Math. Ser., vol. 349, Longman, Harlow, 1996, pp. 6–81. MR 1429633
  * [4] Alexey Chernov and Christoph Schwab, _Exponential convergence of Gauss-Jacobi quadratures for singular integrals over simplices in arbitrary dimension_ , SIAM J. Numer. Anal. 50 (2012), no. 3, 1433–1455. MR 2970750
  * [5] Albert Cohen, Ronald DeVore, and Christoph Schwab, _Convergence rates of best $N$-term Galerkin approximations for a class of elliptic sPDEs_, Found. Comput. Math. 10 (2010), no. 6, 615–646. MR 2728424
  * [6] Albert Cohen, Ronald Devore, and Christoph Schwab, _Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDE’s_ , Anal. Appl. (Singap.) 9 (2011), no. 1, 11–47. MR 2763359
  * [7] Tiangang Cui, Hans De Sterck, Alexander D. Gilbert, Stanislav Polishchuk, and Robert Scheichl, _Multilevel monte carlo methods for stochastic convection-diffusion eigenvalue problems_ , 2023.
  * [8] E. B. Davies, _Spectral theory and differential operators_ , Cambridge Studies in Advanced Mathematics, vol. 42, Cambridge University Press, Cambridge, 1995. MR 1349825
  * [9] Jürgen Dölz and David Ebert, _On uncertainty quantification of eigenpairs with higher multiplicity_ , 2022.
  * [10] S. Giani and I. G. Graham, _Adaptive finite element methods for computing band gaps in photonic crystals_ , Numerische Mathematik 121 (2011), no. 1, 31–64.
  * [11] David Gilbarg and Neil S. Trudinger, _Elliptic partial differential equations of second order_ , Classics in Mathematics, Springer-Verlag, Berlin, 2001, Reprint of the 1998 edition. MR 1814364
  * [12] A. D. Gilbert, I. G. Graham, F. Y. Kuo, R. Scheichl, and I. H. Sloan, _Analysis of quasi-monte carlo methods for elliptic eigenvalue problems with stochastic coefficients_ , Numerische Mathematik 142 (2019), no. 4, 863–915.
  * [13] Alexander D. Gilbert and Robert Scheichl, _Multilevel quasi-monte carlo for random elliptic eigenvalue problems i: Regularity and error analysis_ , 2020\.
  * [14] Alexander D. Gilbert and Robert Scheichl, _Multilevel quasi-monte carlo for random elliptic eigenvalue problems ii: Efficient algorithms and numerical results_ , 2022.
  * [15] Luka Grubišić, Mikael Saarikangas, and Harri Hakula, _Stochastic collocation method for computing eigenspaces of parameter-dependent operators_ , Numer. Math. 153 (2023), no. 1, 85–110. MR 4530202
  * [16] Antoine Henrot, _Extremum problems for eigenvalues of elliptic operators_ , Frontiers in Mathematics, Birkhäuser Verlag, Basel, 2006. MR 2251558
  * [17] V. Kaarnioja, F. Y. Kuo, and I. H. Sloan, _Uncertainty quantification using periodic random variables_ , SIAM Journal on Numerical Analysis 58 (2020), no. 2, 1068–1091.
  * [18] Tosio Kato, _Perturbation theory for linear operators_ , Springer Berlin Heidelberg, 1966.
  * [19] Steven G. Krantz and Harold R. Parks, _A primer of real analytic functions_ , Birkhäuser Boston, 2002.
  * [20] F. Y. Kuo, Ch. Schwab, and I. H. Sloan, _Correction to “Quasi-Monte Carlo methods for high-dimensional integration: the standard (weighted Hilbert space) setting and beyond” [mr2928989]_ , ANZIAM J. 54 (2013), no. 3, 216–219. MR 3085814
  * [21] Frances Y. Kuo and Dirk Nuyens, _Application of quasi-Monte Carlo methods to elliptic PDEs with random diffusion coefficients: a survey of analysis and implementation_ , Found. Comput. Math. 16 (2016), no. 6, 1631–1696. MR 3579719
  * [22] F. Rellich and J. Berkowitz, _Perturbation theory of eigenvalue problems_ , New York University. Institute of Mathematical Sciences, Gordon and Breach, 1969\.
  * [23] Luigi Rodino, _Linear partial differential operators in Gevrey spaces_ , World Scientific Publishing Co., Inc., River Edge, NJ, 1993. MR 1249275

## APPENDIX

Here we collect some combinatorial results needed in the proof of Theorem 5.1.

######  Lemma 6.1.

For all $\boldsymbol{\nu}\in\mathbb{N}^{d}_{0}$ and $0\leq
r\leq\left|\boldsymbol{\nu}\right|$, we have

$\sum_{\left|\boldsymbol{m}\right|=r\atop\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)=\left(\left|\boldsymbol{\nu}\right|\atop
r\right).$

###### Proof.

For any $z\in\mathbb{R}$, we choose functions
$h(\boldsymbol{x}),g(\boldsymbol{x})$ and $f(\boldsymbol{x})$ as following

$\displaystyle h(\boldsymbol{x})=\exp\left(z\sum_{i=1}^{d}x_{i}\right),\quad
g(\boldsymbol{x})=\exp\left(\sum_{i=1}^{d}x_{i}\right)\text{ and
}f(\boldsymbol{x})=h(\boldsymbol{x})\,g(\boldsymbol{x})=\exp\left((z+1)\sum_{i=1}^{d}x_{i}\right)$

where $\boldsymbol{x}\in\mathbb{R}^{d}$. On the one hand, for every
$\boldsymbol{\nu}\in\mathbb{N}^{d}_{0}$, take $\boldsymbol{\nu}$-derivative
with respect to $\boldsymbol{x}$ and evaluate at
$\boldsymbol{x}=\boldsymbol{0}$, we derive

$\displaystyle\partial^{\boldsymbol{\nu}}f(\boldsymbol{0})=(z+1)^{\left|\boldsymbol{\nu}\right|}f(\boldsymbol{0})=\sum_{r=0}^{\left|\boldsymbol{\nu}\right|}\left(\left|\boldsymbol{\nu}\right|\atop
r\right)z^{r}.$

On the other hand, by Leibniz product rule, we get

$\displaystyle\partial^{\boldsymbol{\nu}}f(\boldsymbol{0})=\sum_{\boldsymbol{0}\leq\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\partial^{\boldsymbol{\nu}-\boldsymbol{m}}g(\boldsymbol{0})\partial^{\boldsymbol{m}}h(\boldsymbol{0})=\sum_{\boldsymbol{0}\leq\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)z^{\left|\boldsymbol{m}\right|}=\sum_{r=0}^{\left|\boldsymbol{\nu}\right|}\sum_{\left|\boldsymbol{m}\right|=r\atop\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)z^{r}.$

We finish the proof by comparing the coefficients in front of $z^{r}$ in the
above equations. ∎

######  Lemma 6.2.

For two multiindices
$\boldsymbol{\nu},\boldsymbol{m},\boldsymbol{\ell}\in\mathcal{F}$ and
$\delta\geq 1$ we have

(6.9)
$\sum_{0<\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)[\tfrac{1}{2}]_{|\boldsymbol{\nu}-\boldsymbol{m}|}[\tfrac{1}{2}]_{|\boldsymbol{m}|}\leq
3[\tfrac{1}{2}]_{|\boldsymbol{\nu}|},$ (6.10)
$\sum_{0<\boldsymbol{m}<\boldsymbol{\nu}}\sum_{0\leq\boldsymbol{\ell}\leq\boldsymbol{m}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)\left(\boldsymbol{m}\atop\boldsymbol{\ell}\right)[\tfrac{1}{2}]_{|\boldsymbol{\nu}-\boldsymbol{m}|}[\tfrac{1}{2}]_{|\boldsymbol{m}-\boldsymbol{\ell}|}[\tfrac{1}{2}]_{|\boldsymbol{m}|}\leq
8[\tfrac{1}{2}]_{|\boldsymbol{\nu}|}.$

###### Proof.

For the first term, appying Lemma 6.1 and Lemma 2.3, we arrive at

$\displaystyle\sum_{0<\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop\boldsymbol{m}\right)[\tfrac{1}{2}]_{|\boldsymbol{\nu}-\boldsymbol{m}|}[\tfrac{1}{2}]_{|\boldsymbol{m}|}=\sum_{r=1}^{\left|\boldsymbol{\nu}\right|}\sum_{\left|\boldsymbol{m}\right|=r\atop\boldsymbol{m}\leq\boldsymbol{\nu}}\left(\boldsymbol{\nu}\atop
r\right)[\tfrac{1}{2}]_{|\boldsymbol{\nu}|-r}[\tfrac{1}{2}]_{r}=\sum_{r=1}^{\left|\boldsymbol{\nu}\right|}\left(|\boldsymbol{\nu}|\atop
r\right)[\tfrac{1}{2}]_{|\boldsymbol{\nu}|-r}[\tfrac{1}{2}]_{r}\leq
3[\tfrac{1}{2}]_{|\boldsymbol{\nu}|}.$

The second inequality is proved similarly. We apply Lemma 2.3 twice to obtain
the final result. ∎