# Improvements to the theoretical estimates of the
Schwarz preconditioner with $\Delta$-GenEO coarse space for the
indefinite Helmholtz problem

Victorita Dolean111CASA, TU Eindhoven<EMAIL_ADDRESS>Mark
Fry222Department of Mathematics and Statistics, University of Strathclyde,
Glasgow, G1 1XH, UK<EMAIL_ADDRESS><EMAIL_ADDRESS>Ivan G.
Graham333Department of Mathematical Sciences, University of Bath, Bath, BA2
7AY, UK<EMAIL_ADDRESS>Matthias Langer†

###### Abstract

The purpose of this work is to improve the estimates for the $\Delta$-GenEO
method from the paper [2] when applied to the indefinite Helmholtz equation.
We derive $k$-dependent estimates of quantities of interest ensuring the
robustness of the method.

MSC Classes: 65N22, 65N55, 65F10
Keywords: elliptic PDEs, finite element methods, domain decomposition,
preconditioning, GMRES, Helmholtz equation, spectral coarse space, robustness

##  1 Introduction

In this work we focus on the indefinite version of Helmholtz boundary value
problems with highly variable coefficients defined on
$\Omega\subset\mathbb{R}^{d}$, given by:

$\displaystyle-\operatorname{div}(A\nabla u)-k^{2}u$ $\displaystyle=f\qquad$
$\displaystyle\text{in }\Omega,$ (1.1) $\displaystyle u$
$\displaystyle=0\qquad$ $\displaystyle\text{on }\partial\Omega.$

We use the same setting as in [2], namely we work on a computational domain
$\Omega\subset\mathbb{R}^{d}$ ($d=2,3$) which is supposed to be polygonal or
Lipschitz polyhedral. The matrix $A$ is positive definite, and we assume that
(1.1) has a unique weak solution $u\in H^{1}(\Omega)$ for all $f\in
L^{2}(\Omega)$.

The global domain $\Omega$ is covered by a set of overlapping subdomains
$\Omega_{i}$, $i=1,\ldots,N$, and the classical one-level additive Schwarz
preconditioner is built from partial solutions on each subdomain. Since this
preconditioner is, in general, not scalable as the number of subdomains grows,
an additional global coarse solve is usually added to enhance scalability, as
well as robustness with respect to coefficient heterogeneity or the increase
in the wavenumber. As in [2] we use local generalised eigenvalue problems that
are obtained by a shift of the left-hand side of (1.1). After discretisation
by finite elements, the linear systems arising from (1.1) are symmetric but
indefinite. For this reason, we use GMRES as the iterative solver, and the
convergence analysis relies on the ‘Elman theory’ [3], which requires an upper
bound for the norm of the preconditioned matrix and a lower bound on the
distance of its field of values from the origin. As a result, the number of
GMRES iterations to achieve a given error tolerance is a function of these two
bounds.

The main results of this preprint. Our main theoretical result provides
rigorous and $k$-explicit upper bounds on the coarse mesh diameter $H$ and on
the ‘eigenvalue tolerance’ $\tau$ (for the local generalised eigenvalue
problems) which ensure that GMRES enjoys robust and mesh-independent
convergence when applied to the preconditioned problem. This rigorous explicit
bound for (1.1) follows the methodology from [2] and provides an improvement
of the $k$-dependent estimates. As a reminder, GenEO coarse spaces are usually
based on the $m_{i}$ (dominant) eigenfunctions corresponding to the smallest
eigenvalues
$\lambda^{i}_{1}\leq\lambda^{i}_{2}\leq\cdots\leq\lambda_{m_{i}}^{i}$ of the
generalised eigenvalue problem on the subdomain $\Omega_{i}$. To obtain a
robust rate of convergence for GMRES that depends only on $\Lambda$ (the
maximum number of times any point is overlapped by the subdomains
$\Omega_{i}$), we need some conditions on $H$ and the eigenvalue tolerance
$\tau\mathrel{\mathop{:}}=\min_{i=1}^{N}\lambda^{i}_{m_{i}+1}$, where $N$ is
the number of subdomains. In the estimates $C_{\textup{{stab}}}>0$ the
stability constant for the problem (1.1) also appears. The hidden constants
usually depend only on $\Lambda$.

More explicitly, in [2, Theorem 4.1] it is proved that a robust rate is
achieved if

$H\lesssim
k^{-2}\qquad\text{and}\qquad(1+C_{\textup{{stab}}})^{2}\,k^{8}\lesssim\tau.$
(1.2)

In the current paper, we show that the bounds (1.2) can be improved:

$H\lesssim
k^{-1}\qquad\text{and}\qquad(1+C_{\textup{{stab}}})^{2}k^{4}\lesssim\tau.$
(1.3)

Although, this is a first improvement, these bounds remain quite pessimistic,
showing a strong dependence on the wavenumber $k$.

##  2 Useful Background

###  2.1 Problem formulation and discretisation

The weak formulation of (1.1) is to find $u\in H^{1}_{0}(\Omega)$ such that

$b(u,v)=(f,v)\qquad\text{for all}\ v\in H^{1}_{0}(\Omega),$ (2.1)

where $f\in L^{2}(\Omega)$ and $b(\cdot,\cdot):H^{1}_{0}(\Omega)\times
H^{1}_{0}(\Omega)\rightarrow\mathbb{R}$ is defined as

$b(u,v)=\int_{\Omega}(A\nabla u\cdot\nabla v-k^{2}uv)\,\mathrm{d}x.$

We shall be making use of the positive definite bilinear form
$a(\cdot,\cdot):H^{1}_{0}(\Omega)\times
H^{1}_{0}(\Omega)\rightarrow\mathbb{R}$.

$a(u,v)=\int_{\Omega}A\nabla u\cdot\nabla vdx,$

If $a$ and $b$ are defined on a subdomain $\Omega^{\prime}$ of $\Omega$ we
employ the notations $a_{\Omega^{\prime}}$ and $b_{\Omega^{\prime}}$. The
following weak regularity assumptions are going to be used throughout this
work.

######  Assumption 2.1.

The coefficient $A$ in problem (1.1) satisfies the following conditions.

  1. (i)

$A:\Omega\rightarrow\mathbb{R}^{d\times d}$ is symmetric with
$0<a_{\textup{{min}}}\leq a_{\textup{{max}}}$ s.t.

$a_{\textup{{min}}}|\bm{\xi}|^{2}\leq A(x)\bm{\xi}\cdot{\bm{\xi}}\leq
a_{\textup{{max}}}|\bm{\xi}|^{2}\qquad\text{for all}\
x\in\Omega,\bm{\xi}\in\mathbb{R}^{d}.$ (2.2)

  2. (ii)

Without loss of generality, $a_{\textup{{min}}}=1$ and the diameter,
$D_{\Omega}$, of the domain $\Omega$ is such that $D_{\Omega}\leq 1$. If this
is not the case, then the problem can be scaled accordingly.

######  Notation 2.2.

For any subdomain $\Omega^{\prime}\subset\Omega$, we use
$(\cdot,\cdot)_{\Omega^{\prime}}$ to denote the $L^{2}(\Omega^{\prime})$ inner
product, with the norm denoted by $\|\cdot\|_{\Omega^{\prime}}$. When
$\Omega^{\prime}=\Omega$, the inner product is written as $(\cdot,\cdot)$ with
norm $\|\cdot\|$. The norm induced by the positive bilinear form $a$ is
denoted by $\|u\|_{a,\Omega^{\prime}}=\sqrt{a_{\Omega^{\prime}}(u,u)}$. When
$\Omega^{\prime}=\Omega$ we abandon the subscript $\Omega$. We see that
$b(u,v)=a(u,v)-k^{2}(u,v)$.

The bilinear form $a(\cdot,\cdot)$ is positive definite (SPD), whereas
$b(\cdot,\cdot)$ is symmetric, but in general indefinite. In all what follows,
solvability of (2.1) is assumed.

######  Assumption 2.3.

For any $f\in L^{2}(\Omega)$, it is assumed that the problem 2.1 has a unique
solution $u\in H_{0}^{1}(\Omega)$ and there exists a constant
$C_{\textup{{stab}}}>0$ such that

$\|u\|_{a}\leq C_{\textup{{stab}}}\|f\|\qquad\text{for all}\ f\in
L^{2}(\Omega).$ (2.3)

Let $\mathcal{T}_{h}$ be any shape regular triangular mesh over the domain,
$\Omega$. For the purpose of this work, 2 or 3-dimensional simplices are being
considered, but this could easily be applied to $d$-dimensional simplices,
where the maximum diameter is $h$. Let $V^{h}\subset H^{1}_{0}(\Omega)$ be any
conforming finite element space. The Galerkin approximation of (2.1) is to
find $u_{h}\in V^{h}$ such that

$b(u_{h},v)=(f,v)\qquad\text{for all}\ v\in V^{h}.$ (2.4)

If $n$ denotes the dimension of $V^{h}$, with a basis given by
$\\{\phi_{i}\\}_{i=1}^{n}$, then (2.4) can be represented by the linear system

$\mathbf{B}\mathbf{u}=\mathbf{f}.$ (2.5)

The matrix $\mathbf{B}$ is defined in terms of the basis functions as
$(\mathbf{B})_{ij}=b(\phi_{j},\phi_{i})$ and $(\mathbf{f})_{i}=(f,\phi_{i})$.
It is also possible to define the matrix $\mathbf{A}$ using the same basis
functions as $(\mathbf{A})_{ij}=a(\phi_{j},\phi_{i})$.

The solvability of (2.4) is assured by the following Lemma from [5, Theorem
2]. This is required due to the indefiniteness of (2.4).

######  Lemma 2.4 (Schatz and Wang, 1996).

Let Assumptions 2.1 and 2.3 hold. Then there exists an $h_{0}>0$ such that,
for each $h$ with $0<h<h_{0}$, the problem (2.4) has a unique solution
$u_{h}\in V_{H}$. Moreover, let $u$ be the unique solution of (2.1). Then, for
every $\varepsilon>0$ there exists $h_{1}=h_{1}(\varepsilon)>0$ such that, for
every $h\in(0,h_{1})$,

$\|u-u_{h}\|\leq\varepsilon\|u-u_{h}\|_{{H^{1}}(\Omega)}$ (2.6)

and

$\|u-u_{h}\|_{H^{1}(\Omega)}\leq\varepsilon\|f\|.$ (2.7)

######  Remark 2.5.

The introduction of (2.6) differs from [2], allowing for improvements to be
made in later Lemmas.

By combining the Friedrichs inequality [4, Theorem 13.19] with Assumption 2.1,
it possible to state that, for any subdomain $\Omega^{\prime}\subset\Omega$
with diameter $H$, we obtain the estimate

$\|u\|_{\Omega^{\prime}}\leq\frac{H}{\sqrt{2}\,}\|\nabla
u\|_{\Omega^{\prime}}\leq\frac{H}{\sqrt{2}\,}\|u\|_{a_{\Omega^{\prime}}}\qquad\text{for
all}\ u\in H^{1}_{0}(\Omega^{\prime}).$ (2.8)

This inequality will be used extensively in what follows.

###  2.2 Domain decomposition

In order to construct the two-level Schwarz preconditioner, the first-level
preconditioner needs to be formulated. This is achieved by first partitioning
the global domain, $\Omega$, into a set of $N$ subdomains,
$\\{\Omega_{i}^{\prime}\\}_{i=1}^{N}$. It is being assumed that the global
mesh, $\mathcal{T}_{h}$, is sufficiently fine to resolve the subdomains. Then
each subdomain, $\Omega^{\prime}$, is expanded by one or more layers of mesh
elements, in the sense of Definition 2.6, creating the desired overlapping set
of $N$ subdomains, $\\{\Omega_{i}\\}_{i=1}^{N}$.

######  Definition 2.6.

Given a subdomain $\Omega^{\prime}\subset\Omega$, which is resolved by the
chosen mesh, the extension of $\Omega^{\prime}$ by a layer of elements is

$\Omega_{\textsf{e}}=\operatorname{Int}\Biggl{(}\bigcup_{\\{\ell\mid\operatorname{supp}(\phi_{\ell})\cap\Omega^{\prime}\not=\emptyset\\}}\operatorname{supp}(\phi_{\ell})\Biggr{)}$

where $\operatorname{Int}(\cdot)$ is the interior of domain. Extension by
multiple layers can then be achieved by applying this recursively.

For the domains $i=1,\ldots,N$, it is possible to define the spaces,

$\widetilde{V}_{i}=\\{v|_{\Omega_{i}}:v\in V^{h}\\}\subset
H^{1}(\Omega_{i})\quad\text{and}\quad
V_{i}=\\{v\in\widetilde{V}_{i}:v|_{\partial\Omega_{i}}=0\\}\subset
H^{1}_{0}(\Omega_{i}).$ (2.9)

For each subdomain, $\Omega_{i}$, we denote its diameter by $H_{i}$, and we
set $H\mathrel{\mathop{:}}=\max\\{H_{i}\\}$. For any
$u,v\in\widetilde{V}_{i}$, the local bilinear forms can be defined

$a_{\Omega_{i}}(u,v)=\int_{\Omega_{i}}A\nabla u\cdot\nabla
v\,\mathrm{d}x,\qquad b_{\Omega_{i}}(u,v)=\int_{\Omega_{i}}(A\nabla
u\cdot\nabla v-k^{2}uv)\,\mathrm{d}x.$

######  Remark 2.7.

Although $b_{\Omega_{i}}(\cdot,\cdot)$ is generally indefinite, for small
enough $H$ the form restricted to $V_{i}$ is positive definite. This is
discussed in Lemma 3.8.

For any $v_{i}\in V_{i}$, let $E_{i}v_{i}$ denote its zero extension to the
whole of the domain $\Omega$. Then,

$E_{i}:V_{i}\rightarrow V^{h},\quad i=1,\ldots,N.$ (2.10)

The $L^{2}(\Omega)$ adjoint of the extension operator is called restriction
operator,

$R_{i}:V^{h}\rightarrow V_{i}.$

By use of the extension operator, the restriction of the bilinear forms to
$V_{i}$ can be given as

$a_{\Omega_{i}}(u,v)=a(E_{i}u,E_{i}v),\qquad
b_{\Omega_{i}}(u,v)=b(E_{i}u,E_{i}v),\qquad(u,v)_{\Omega_{i}}=(E_{i}u,E_{i}v)$

for all $u,v\in V_{i}$. The one-level additive Schwarz preconditioner can now
be given in matrix form as

$\mathbf{M}^{-1}_{AS,1}=\sum^{N}_{i=1}\mathbf{E}_{i}\mathbf{B}^{-1}_{i}\mathbf{R}_{i},\qquad\text{where
}\mathbf{B}_{i}=\mathbf{R}_{i}\mathbf{B}\mathbf{E}_{i}.$ (2.11)

Here, $\mathbf{E}_{i}$ and $\mathbf{R}_{i}$ denote the matrix representations
of $E_{i}$ and $R_{i}$ with respect to the basis functions
$\\{\phi_{i}\\}_{i=1}^{n}$ and some basis in $V_{i}$.

In order to improve the effectiveness of the preconditioner, a coarse space is
added. This improves the global exchange of information between the
subdomains. Let $V_{0}\subset V^{h}$ be such a coarse space. Let
$E_{0}:V_{0}\rightarrow V^{h}$ be the natural embedding, and let $R_{0}$ be
the $L^{2}$ adjoint of $E_{0}$,

$(v_{0},R_{0}w)=(E_{0}v_{0},w)\qquad\text{for all}\ w\in V^{h},v_{0}\in
V_{0}.$

The two-level additive Schwarz preconditioner can now be given in matrix form
as,

$\mathbf{M}^{-1}_{AS,2}=\sum^{N}_{i=0}\mathbf{E}_{i}\mathbf{B}^{-1}_{i}\mathbf{R}_{i},\qquad\text{where
}\mathbf{B}_{i}=\mathbf{R}_{i}\mathbf{B}\mathbf{E}_{i}.$ (2.12)

The preconditioned linear system from (2.5) reads as,

$\mathbf{M}^{-1}_{AS,2}\mathbf{B}\mathbf{u}=\mathbf{M}_{AS,2}^{-1}\mathbf{f}$
(2.13)

It is now possible to define the projectors used in the analysis. For each
$i=0,\ldots,N$, the projectors, $T_{i}:V^{h}\rightarrow V_{i}$, are defined
by,

$b_{\Omega_{i}}(T_{i}u,v)=b(u,E_{i}v)\qquad\text{for all}\ v\in V_{i},$ (2.14)

where $\Omega_{0}=\Omega$. The existence of the $T_{i}$ operators are
guaranteed by Lemma 3.8. Given the operators $T_{i}$, the operator
$T:V^{h}\rightarrow V^{h}$ is defined as

$T=\sum_{i=0}^{N}E_{i}T_{i}.$ (2.15)

This allows for the two-level additive Schwarz preconditioner to be
represented in terms of the projector operator, $T$, as follows. We recall
here the result from [2].

######  Proposition 2.8.

For any $u,v\in V^{h}$, with corresponding nodal vectors
$\mathbf{u},\mathbf{v}\in\mathbb{R}^{d}$,

$\langle\mathbf{M}_{AS,2}^{-1}\mathbf{B}\mathbf{u},\mathbf{v}\rangle_{\mathbf{A}}=a(Tu,v),$
(2.16)

where $\langle\cdot,\cdot\rangle_{\mathbf{A}}$ is the inner product on
$\mathbb{R}^{d}$, using the matrix $\mathbf{A}$.

###  2.3 The $\Delta$-GenEO coarse space

In order to use the $\Delta$-GenEO coarse space, we need to recall some
definitions from [6].

######  Definition 2.9 ([6] Definition 3.2).

Given a subdomain $\Omega_{i}$, which is formed from a union of elements from
the whole domain, let

$\overline{\operatorname{dof}}(\Omega_{i})\mathrel{\mathop{:}}=\bigl{\\{}\ell\mid
1\leq\ell\leq n\text{ and
}\operatorname{supp}(\phi_{\ell})\cap\Omega_{i}\neq\emptyset\bigr{\\}}$

denote the set of degrees of freedom that are active in the subdomain
$\Omega_{i}$, including ones on the boundary. In a similar manner, let

$\operatorname{dof}(\Omega_{i})\mathrel{\mathop{:}}=\bigl{\\{}\ell\mid
1\leq\ell\leq n\text{ and
}\operatorname{supp}(\phi_{\ell})\subset\overline{\Omega_{i}}\bigr{\\}}$

denote the internal degrees of freedom.

######  Definition 2.10 (Partition of unity).

Let $\operatorname{dof}(\Omega_{i})$ be as in Definition 2.9. For any degree
of freedom, $j\in\\{1,\ldots,n\\}$, let $\mu_{j}$ denote the number of
subdomains for which $j$ is an internal degree of freedom, i.e.

$\mu_{j}\mathrel{\mathop{:}}=\\#\bigl{\\{}i\mid 1\leq i\leq
N,\,j\in\operatorname{dof}(\Omega_{i})\bigr{\\}}.$

Then, for $i\in\\{1,\ldots,N\\}$, the local partition of unity operator,
$\Xi_{i}:\widetilde{V}_{i}\rightarrow V_{i}$ is defined by

$\Xi_{i}(v)\mathrel{\mathop{:}}=\sum_{j\in\operatorname{dof}(\Omega_{i})}\frac{1}{\mu_{j}}v_{j}\phi_{j}^{i}\qquad\text{for}\
v=\sum_{j\in\overline{\operatorname{dof}}(\Omega_{i})}v_{j}\phi_{j}^{i}\in\widetilde{V}_{i}.$

The local generalised eigenvalue problem that is going to form the basis of
the coarse space can now be introduced.

######  Definition 2.11 ([2], Definition 2.8).

For each $j=1,\ldots,N$, we define the following generalised eigenvalue
problem. Find
$(p_{\ell}^{j},\lambda_{j})\in\widetilde{V}_{j}\setminus\\{0\\}\times\mathbb{R}$
such that

$a_{\Omega_{j}}(p_{\ell}^{j},v)=\lambda_{j}a_{\Omega_{j}}\bigl{(}\Xi_{j}(p_{\ell}^{j}),\Xi_{j}(v)\bigr{)}\qquad\text{for
all}\ v\in\widetilde{V}_{j}.$

######  Remark 2.12.

As $a_{\Omega_{j}}(\Xi_{j}(\cdot),\Xi_{j}(\cdot))$ is symmetric positive
definite on $V_{j}$, the finite eigenvectors, $p_{\ell}^{j}$, can be
normalised with respect to $a_{\Omega_{j}}(\Xi_{j}(\cdot),\Xi_{j}(\cdot))$ to
satisfy the orthogonality conditions

$a_{\Omega_{j}}\bigl{(}\Xi_{j}(p_{m}^{j}),\Xi_{j}(p_{\ell}^{j})\bigr{)}=\delta_{m}^{\ell}\quad\text{
and }\quad
a_{\Omega_{j}}\bigl{(}p_{m}^{j},p_{\ell}^{j}\bigr{)}=\lambda_{j}\delta_{m}^{\ell}.$
(2.17)

For the treatment of the eigenvalue infinity see the discussion in [6].

######  Definition 2.13 ($\Delta$-GenEO Coarse space).

For each $j\in 1,\ldots,N$, let $\\{p_{\ell}^{j}\\}^{m_{j}}_{\ell=1}$ be the
eigenfunctions corresponding to the $m_{j}$ smallest eigenvalues from
Definition 2.11. The value of $m_{j}$ is to be chosen, with more details on
this later. The coarse space, $V_{0}$, is given by

$V_{0}\mathrel{\mathop{:}}=\operatorname{span}\\{E_{j}\Xi_{j}(p_{\ell}^{j})\mid
l=1,\ldots,m_{j}\text{ and }j=1,\ldots,N\\}.$ (2.18)

We will also need the following notations

$\Lambda\mathrel{\mathop{:}}=\max_{T\in\mathcal{T}_{h}}\\!\left(\\#\\{\Omega_{j}\mid
1\leq j\leq
N,T\subset\Omega_{j}\\}\right),\qquad\tau\mathrel{\mathop{:}}=\min_{1\leq
j\leq N}\lambda^{j}_{m_{j}+1}.$

##  3 Statement of the main result and theoretical tools

We start this section with the statement of the main result and continue by
providing a few technical lemmas needed in the proof of this result, which
will be given later on. For convenience we set
$\Theta\mathrel{\mathop{:}}=\frac{1}{\tau}$.

######  Theorem 3.1 (GMRES convergence of the two-level preconditioned
system).

Assume that $h\in(0,h_{1})$, where $h_{1}$ is as in Lemma 3.11, then select
$H$ and $\tau$ such that

$\displaystyle s$
$\displaystyle\mathrel{\mathop{:}}=(1+\Lambda^{2}\Theta)16\sqrt{2}k^{2}\Lambda^{\frac{3}{2}}\Theta^{\frac{1}{2}}(1+C_{\textup{{stab}}})<1,$
(3.1) $\displaystyle t$
$\displaystyle\mathrel{\mathop{:}}=32H^{2}k^{2}(1+\Lambda^{2}\Theta)\Lambda<1.$

If GMRES is applied in the $\langle\cdot,\cdot\rangle_{\mathbf{A}}$-inner
product to solve the preconditioned system given by (2.13), then after $m$
iterations, the norm of the residual, ${\bf r}^{(m)}$, is bounded as

$\|{\bf
r}^{(m)}\|_{\mathbf{A}}^{2}\leq\biggl{(}1-\frac{c_{1}^{2}}{c_{2}^{2}}\biggr{)}^{m}\|{\bf
r}^{(0)}\|_{\mathbf{A}}^{2},$ (3.2)

where $c_{1}$ and $c_{2}$ are given by

$c_{1}\mathrel{\mathop{:}}=\\\
bigl(4(1+\Lambda^{2}\Theta)\bigr{)}^{-1}\bigl{(}1-\max\\{t,s\\}\bigr{)},\qquad
c_{2}\mathrel{\mathop{:}}=18+18\Lambda^{2}.$ (3.3)

######  Corollary 3.2.

If (3.1) in Theorem 3.1 are satisfied, then this leads to the conditions

$H\lesssim
k^{-1}\quad\text{and}\quad(1+C_{\textup{{stab}}})^{2}k^{4}\lesssim\tau.$ (3.4)

If these conditions are satisfied, then $c_{1}$ and $c_{2}$ are both
independent on problem parameters including the heterogeneity and the
wavenumber $k$, leading to a robust GMRES convergence.

###### Proof.

It is known that $\Lambda\geq 1$; so we have $1\leq\Lambda\leq\Lambda^{3/2}$
and

$16\sqrt{2}k^{2}\Lambda\Theta^{\frac{1}{2}}(1+C_{\textup{{stab}}})\leq(1+\Lambda^{2}\Theta)16\sqrt{2}k^{2}\Lambda^{\frac{3}{2}}\Theta^{\frac{1}{2}}(1+C_{\textup{{stab}}})<1,$

thus leading to the inequality

$2k^{2}\Lambda\Theta^{\frac{1}{2}}(1+C_{\textup{{stab}}})<\frac{1}{8\sqrt{2}}.$
(3.5)

If we use the definition of $\tau$, then (3.5) translates into

$(1+C_{\textup{{stab}}})^{2}k^{4}\lesssim\tau.$

Now, using the second inequality from (3.1), we obtain

$64H^{2}k^{2}\leq 64H^{2}k^{2}(1+\Lambda^{2}\Theta)\Lambda<1,$

which leads to the condition $H\lesssim k^{-1}$. ∎

In practice, conditions (3.4) will introduce a constraint in the size of the
subdomains, depending on $k$ and on the number of modes to be added in the
coarse space.

######  Remark 3.3.

Note that the condition (3.1) is different from the one obtained in [2], and
allows for the improved conditions on $H$ and $\tau$. These results follow
from the modifications made to the later lemmas.

###  3.1 Properties of the $\Delta$-GenEO coarse space

The following three lemmas and Proposition 3.7 are also given in [2] and
adapted from [6]. The first gives an error estimate for the local projection
operator, which is used to approximate a function $v\in\widetilde{V}_{i}$ in
the space being spanned by the eigenfunctions from Definition 2.11.

######  Lemma 3.4 (Projection operator onto the coarse space).

Let $i\in\\{1,\ldots,N\\}$ and $\\{(p_{\ell}^{i},\lambda_{\ell}^{i})\\}$ be
the eigenpairs of the generalised eigenproblem, given in Definition 2.11.
Suppose that $m_{i}\in\\{1,\ldots,\dim(V_{i})-1\\}$ is such that
$0<\lambda_{m_{i}+1}^{i}<\infty$, and the eigenvalues can be ordered such that
$\lambda_{1}^{i}\leq\lambda_{2}^{i}\leq\ldots\leq\lambda_{\dim(\widetilde{V}_{i})}^{i}$.
Then the local projector, $\Pi_{m_{i}}^{i}$, defined by

$\Pi_{m_{i}}^{i}v\mathrel{\mathop{:}}=\sum_{l=1}^{m_{i}}a_{\Omega_{i}}\bigl{(}\Xi_{i}(v),\Xi_{i}(p_{\ell}^{i})\bigr{)}p_{\ell}^{i},$
(3.6)

satisfies

$\|w\|_{a_{\Omega_{i}}}^{2}\leq\|v\|^{2}_{a_{\Omega_{i}}}\quad\text{ and
}\quad\|\Xi(w)\|^{2}_{a_{\Omega_{i}}}\leq\frac{1}{\lambda_{m_{i}+1}^{i}}\|w\|_{a_{\Omega_{i}}}^{2}\qquad\text{for
all}\ v\in\widetilde{V}_{i},$ (3.7)

where $w=v-\Pi_{m_{i}}^{i}v$.

From these local error estimates it is possible to build a global
approximation property.

######  Lemma 3.5 (Global approximation property).

Assume that the conditions in Lemma 3.4 are satisfied. Let $v\in V^{h}$; then

$\inf_{z\in
V_{0}}\left\|v-z\right\|_{a}^{2}\leq\left\|v-z_{0}\right\|_{a}^{2}\leq\Lambda^{2}\Theta\|v\|^{2}_{a},$
(3.8)

where

$z_{0}=\sum_{i=1}^{N}E_{i}\Xi_{i}\\!\left(\Pi_{m_{i}}^{i}v|_{\Omega_{i}}\right)\qquad\text{and}\qquad\Theta=\max_{1\leq
i\leq N}\frac{1}{\lambda_{m_{i}+1}^{i}}.$

The following lemma shows that the GenEO coarse space, in combination with the
local finite element space, allows for a stable decomposition. This property
is a key component for bounding the condition number in the two-level Schwarz
preconditioner for positive definite cases, as used in [6].

######  Lemma 3.6 (Stable decomposition).

Let $v\in V^{h}$. Then the decomposition

$z_{0}=\sum_{i=1}^{N}E_{i}\Xi_{i}\bigl{(}\Pi_{m_{i}}^{i}{v|_{\Omega_{i}}}\bigr{)}\quad\text{
and }\quad
z_{i}=\Xi_{i}\bigl{(}v|_{\Omega_{i}}-\Pi_{m_{i}}^{i}v|_{\Omega_{i}}\bigr{)}\qquad\text{for}\
i=1,\ldots,N$

satisfies $v=\sum_{i=0}^{N}E_{i}z_{i}$ and

$\|z_{0}\|_{a}^{2}+\sum_{i=1}^{N}\|z_{i}\|_{a_{\Omega_{i}}}^{2}\leq
4(1+\Lambda^{2}\Theta)\|v\|^{2}_{a}.$

It is convenient at this time to introduce the projection operators
$P_{i}:V^{h}\rightarrow V_{i}$ such that, for each $i=0,\ldots,N$,

$a_{\Omega_{i}}(P_{i}u,v)=a(u,E_{i}v)\text{for all}\ v\in V_{i}.$ (3.9)

In [7, Section 2.2] it is proved that these operators are well defined. Using
these operator $P_{i}$ we can defines $P:V^{h}\rightarrow V^{h}$ as

$P=\sum_{i=0}^{N}E_{i}P_{i}.$ (3.10)

######  Proposition 3.7.

Under the same assumptions as in Lemma 3.6, any $u\in V^{h}$ satisfies

$\|u\|^{2}_{a}\leq 4(1+\Lambda^{2}\Theta)(Pu,u)_{a}$ (3.11)

and

$\sum^{N}_{i=0}\|P_{i}u\|^{2}_{a_{\Omega_{i}}}\leq(\Lambda+1)\|u\|^{2}_{a}.$
(3.12)

###  3.2 Solvability and stability of $T_{i},i=1,\ldots,N$

The next two lemmas are reformulations of those used in [2] when applied to
(1.1), which is a particular case of the problem from [2]. These lemmas are
necessary even if there are no improvements of the estimates at this stage.

######  Lemma 3.8 ($T_{i}$ is well defined for $i=1,\ldots,N$).

If $Hk<\sqrt{2}$, then $b_{\Omega_{i}}(\cdot,\cdot)$ is positive definite on
$H_{0}^{1}(\Omega_{i})$ and the operators $T_{i}$, $i=1,\ldots,N$ are well
defined.

###### Proof.

Using the definition of $b_{\Omega_{i}}(\cdot,\cdot)$ and the Friedrichs
inequality we obtain

$b_{\Omega_{i}}(u,u)=a_{\Omega_{i}}(u,u)-k^{2}(u,u)_{\Omega_{i}}\geq\frac{2}{H^{2}}\|u\|_{\Omega_{i}}^{2}-k^{2}\|u\|_{\Omega_{i}}^{2}=\frac{2-H^{2}k^{2}}{H^{2}}\|u\|_{\Omega_{i}}^{2}.$

If the condition $Hk<\sqrt{2}$ is satisfied, then
$b_{\Omega_{i}}(\cdot,\cdot)$ is positive definite. By the Lax–Milgram lemma,
this ensures that $T_{i}$ is well defined. ∎

######  Remark 3.9.

Whilst this a sufficient condition for $T_{i}$ to be well defined, it is not a
necessary condition. In general, if $k^{2}$ is not an eigenvalue of the matrix
corresponding to $a_{\Omega_{i}}$, then $T_{i}$ is well defined.

In order to show the robustness of the two-level method additive Schwarz
method, it is necessary to prove the stability estimates for the $T_{i}$
operators.

######  Lemma 3.10 (Stability of $T_{i}$, $i=1,\ldots,N$).

Suppose that $Hk\leq 1$. Then, for all $u\in V^{h}$,

$\|T_{i}u\|_{a_{\Omega_{i}}}\leq
2\|u\|_{a_{\Omega_{i}}}+k\sqrt{2}\|u\|_{\Omega_{i}}.$ (3.13)

###### Proof.

Using the defining relation
$b_{\Omega_{i}}(T_{i}u,T_{i}u)=b_{\Omega_{i}}(u,T_{i}u)$ and the
Cauchy–Schwarz inequality we obtain

$\displaystyle\|T_{i}u\|^{2}_{a_{\Omega_{i}}}$
$\displaystyle=a_{\Omega_{i}}(T_{i}u,T_{i}u)=b_{\Omega_{i}}(T_{i}u,T_{i}u)+k^{2}(T_{i}u,T_{i}u)_{\Omega_{i}}$
$\displaystyle=b_{\Omega_{i}}(u,T_{i}u)+k^{2}(T_{i}u,T_{i}u)_{\Omega_{i}}$
$\displaystyle=a_{\Omega_{i}}(u,T_{i}u)-k^{2}(u,T_{i}u)_{\Omega_{i}}+k^{2}(T_{i}u,T_{i}u)_{\Omega_{i}}$
$\displaystyle\leq\|u\|_{a_{\Omega_{i}}}\|T_{i}u\|_{a_{\Omega_{i}}}+k^{2}\|u\|_{\Omega_{i}}\|T_{i}u\|_{\Omega_{i}}+k^{2}\|T_{i}u\|^{2}_{\Omega_{i}}$
$\displaystyle\leq\|u\|_{a_{\Omega_{i}}}\|T_{i}u\|_{a_{\Omega_{i}}}+k^{2}\frac{H}{\sqrt{2}\,}\|u\|_{\Omega_{i}}\|T_{i}u\|_{a_{\Omega_{i}}}+k^{2}\frac{H^{2}}{2}\|T_{i}u\|^{2}_{a_{\Omega_{i}}}.$

After simplification, this can be rewritten as

$\biggl{(}1-k^{2}\frac{H^{2}}{2}\biggr{)}\|T_{i}u\|_{a_{\Omega_{i}}}\leq\|u\|_{a_{\Omega_{i}}}+k^{2}\frac{H}{\sqrt{2}\,}\|u\|_{\Omega_{i}}.$

Together with the assumption $Hk\leq 1$, this implies the claimed inequality
(3.13). ∎

###  3.3 Solvability and stability of $T_{0}$

In this section we derive sufficient conditions for solvability and stability
of the coarse space operator $T_{0}$. The conditions obtained from the
subsequent Lemmas 3.11 and 3.12 are an improvement on those found in their
counterparts from [2].

First, to ensure that $T_{0}$ is well defined, a condition on $\Theta$ is
required.

######  Lemma 3.11 ($T_{0}$ is well defined).

Suppose that

$2k\Lambda\Theta^{\frac{1}{2}}(1+C_{\textup{{stab}}})<1.$ (3.14)

Then there exists $h_{1}>0$ such that, for all $h\in(0,h_{1})$, the operator
$T_{0}$ is well defined.

###### Proof.

Assume that there exists a $w_{0}\in V_{0}\backslash\\{0\\}$ such that

$b(w_{0},z)=0\qquad\text{for all}\ z\in V_{0}.$ (3.15)

Let $w\in H^{1}_{0}(\Omega)$ be the solution of

$b(w,v)=(w_{0},v)\qquad\text{for all}\ v\in H^{1}_{0}(\Omega).$

Let $\varepsilon>0$ such that $k^{2}\varepsilon^{2}\leq\frac{1}{2}$. We shall
impose an additional constraint on $\varepsilon$ later. Then Lemma 2.4 implies
that there exists $h_{1}>0$ such that, for all $h\in(0,h_{1})$, there exists a
solution $w_{h}\in V^{h}$ of

$b(w_{h},v)=(w_{0},v)\qquad\text{for all}\ v\in V^{h},$

and (2.6) and (2.7) hold. As $w_{0}\in V_{0}\backslash\\{0\\}$, we can choose
$v=w_{0}$. Then, for all $z\in V_{0}$,

$\displaystyle\|w_{0}\|^{2}$
$\displaystyle=b(w_{h},w_{0})=b(w_{h},w_{0})-b(z,w_{0})=b(w_{h}-z,w_{0})$
$\displaystyle=a(w_{h}-z,w_{0})-k^{2}(w_{h}-z,w_{0})\leq|a(w_{h}-z,w_{0})|+k^{2}|(w_{h}-z,w_{0})|$
$\displaystyle\leq\|w_{h}-z\|_{a}\|w_{0}\|_{a}+k^{2}\|w_{h}-z\|\|w_{0}\|\leq\|w_{h}-z\|_{a}\|w_{0}\|_{a}+k^{2}\varepsilon\|w_{h}-z\|_{H^{1}(\Omega)}\|w_{0}\|$
$\displaystyle\leq\|w_{h}-z\|_{a}\|w_{0}\|_{a}+k^{2}\varepsilon^{2}\|w_{0}\|^{2}$
$\displaystyle\leq\|w_{h}-z\|_{a}\|w_{0}\|_{a}+\frac{1}{2}\|w_{0}\|^{2}$

and hence

$\|w_{0}\|^{2}\leq 2\|w_{h}-z\|_{a}\|w_{0}\|_{a}.$

As this is true for all $z\in V_{0}$, it follows that

$\|w_{0}\|^{2}\leq 2\|w_{0}\|_{a}\inf_{z\in V_{0}}\|w_{h}-z\|_{a}.$

Combined with Lemma 3.5 this yields

$\|w_{0}\|^{2}\leq 2\Lambda\Theta^{\frac{1}{2}}\|w_{0}\|_{a}\|w_{h}\|_{a}.$
(3.16)

Choosing $z=w_{0}$ in (3.15) we have

$0=b(w_{0},w_{0})=\|w_{0}\|^{2}_{a}-k^{2}\|w_{0}\|^{2}$ (3.17)

and hence $\|w_{0}\|_{a}=k\|w_{0}\|$. Together with (3.16) this gives

$\displaystyle\|w_{0}\|^{2}\leq
2k\Lambda\Theta^{\frac{1}{2}}\|w_{0}\|\|w_{h}\|_{a}.$ (3.18)

We can combine (2.3) and (2.7) to obtain

$\displaystyle\|w_{h}\|_{a}$
$\displaystyle\leq\|w\|_{a}+\|w-w_{h}\|_{a}\leq\|w\|_{a}+a_{\textup{{max}}}\|w-w_{h}\|_{H^{1}(\Omega)}$
$\displaystyle\leq
C_{\textup{{stab}}}\|w_{0}\|+a_{\textup{{max}}}\varepsilon\|w_{0}\|\leq(C_{\textup{{stab}}}+1)\|w_{0}\|$
(3.19)

if $\varepsilon$ is chosen small enough so that
$a_{\textup{{max}}}\varepsilon\leq 1$. We now obtain from (3.19) with (3.18)
that

$\|w_{0}\|^{2}\leq
2k\Lambda\Theta^{\frac{1}{2}}(1+C_{\textup{{stab}}})\|w_{0}\|^{2}.$

With the assumption (3.14) this leads to a contradiction. ∎

Another requirement is to find the stability conditions for $T_{0}$.

######  Lemma 3.12 (Stability of $T_{0}$).

Suppose that (3.14) is satisfied. Then there exists $h_{1}>0$ such that, for
$h\in(0,h_{1})$,

$\|T_{0}u-u\|\leq
2\Lambda\Theta^{\frac{1}{2}}(1+C_{\textup{{stab}}})\|T_{0}u-u\|_{a}\qquad\text{for
all}\ u\in V_{h}.$ (3.20)

Suppose, in addition, that

$2k^{2}\Lambda\Theta^{\frac{1}{2}}(1+C_{\textup{{stab}}})\leq\frac{1}{2}\,.$
(3.21)

Then

$\|u-T_{0}u\|\leq\sqrt{2}\|u\|_{a}\qquad\text{for all}\ u\in V^{h}.$ (3.22)

###### Proof.

Under condition (3.14), Lemma 3.11 ensures the existence of $h_{1}>0$ such
that, for $h\in(0,h_{1})$, the operator $T_{0}:V^{h}\rightarrow V_{0}$ is well
defined. Consider the auxiliary problem

$\text{find}\;\;w_{h}\in V^{h}\;\;\text{such that}\qquad
b(w_{h},v)=(T_{0}u-u,v)\quad\text{for all}\ v\in V^{h},$ (3.23)

which has a unique solution by Lemma 2.4. The definition of $T_{0}$ implies
that $b(T_{0}u-u,z)=0$ for all $z\in V_{0}$. We can choose $v=T_{0}u-u$ in
(3.23), which yields

$\displaystyle\|T_{0}u-u\|^{2}$
$\displaystyle=b(w_{h},T_{0}u-u)=b(w_{h},T_{0}u-u)-b(z,T_{0}u-u)=b(w_{h}-z,T_{0}u-u)$
$\displaystyle=a(w_{h}-z,T_{0}u-u)-k^{2}(w_{h}-z,T_{0}u-u)$
$\displaystyle\leq|a(w_{h}-z,T_{0}u-u)|+k^{2}|(w_{h}-z,T_{0}u-u)|$
$\displaystyle\leq\|w_{h}-z\|_{a}\|T_{0}u-u\|_{a}+k^{2}\|w_{h}-z\|\,\|T_{0}u-u\|$
$\displaystyle\leq\|w_{h}-z\|_{a}\|T_{0}u-u\|_{a}+k^{2}\varepsilon\|w_{h}-z\|_{H^{1}(\Omega)}\|T_{0}u-u\|$
$\displaystyle\leq\|w_{h}-z\|_{a}\|T_{0}u-u\|_{a}+k^{2}\varepsilon^{2}\|T_{0}u-u\|^{2},$

where we used Lemma 2.4 for the last two inequalities. Choosing
$\varepsilon>0$ such that $k^{2}\varepsilon^{2}\leq\frac{1}{2}$ we obtain

$\|T_{0}u-u\|\leq 2\|w_{h}-z\|_{a}\|T_{0}u-u\|_{a}.$ (3.24)

As this is true for all $z\in V_{0}$, it follows Lemma 3.5 that

$\|T_{0}u-u\|^{2}\leq 2\|T_{0}u-u\|_{a}\inf_{z\in V_{0}}\|w_{h}-z\|_{a}\leq
2\Lambda\Theta^{\frac{1}{2}}\|T_{0}u-u\|_{a}\|w_{h}\|_{a}.$ (3.25)

If we choose $\varepsilon$ small enough so that
$a_{\textup{{max}}}\varepsilon\leq 1$, we can use a similar calculation as in
(3.19) to get $\|w_{h}\|_{a}\leq(1+C_{\textup{{stab}}})\|T_{0}u-u\|$. Combined
with (3.25) this proves (3.20).

We now come to the proof of (3.22). From the definition of $P_{0}$ in (3.9) we
have

$a(P_{0}u-u,v)=0\qquad\text{for all}\ u\in V^{h},v\in V_{0}.$

Taking $v=T_{0}u\in V_{0}$, we get $a(T_{0}u,P_{0}u-u)=0$. Moreover, since
$P_{0}u-T_{0}u\in V_{0}$, the definition of $T_{0}$ yields

$b(u-T_{0}u,P_{0}u-T_{0}u)=0.$ (3.26)

Using (3.26), the link between the bilinear forms $a$ and $b$, and the
Cauchy–Schwarz inequality, we obtain

$\displaystyle\|u-T_{0}u\|_{a}^{2}$
$\displaystyle=a(u-T_{0}u,u-T_{0}u)=b(u-T_{0}u,u-T_{0}u)-k^{2}(u-T_{0}u,u-T_{0}u)$
$\displaystyle=b(u-T_{0}u,u-T_{0}u)-b(u-T_{0}u,P_{0}u-T_{0}u)-k^{2}(u-T_{0}u,u-T_{0}u)$
$\displaystyle=b(u-T_{0}u,u-P_{0}u)-k^{2}(u-T_{0}u,u-T_{0}u)$
$\displaystyle=a(u-T_{0}u,u-P_{0}u)+k^{2}(u-T_{0}u,u-P_{0}u)-k^{2}(u-T_{0}u,u-T_{0}u)$
$\displaystyle=a(u-T_{0}u,u-P_{0}u)+k^{2}(u-T_{0}u,T_{0}u-P_{0}u)$
$\displaystyle=a(u-T_{0}u,u-P_{0}u)+a(T_{0}u,u-P_{0}u)+k^{2}(u-T_{0}u,T_{0}u-P_{0}u)$
$\displaystyle=a(u,u-P_{0}u)+k^{2}(u-T_{0}u,T_{0}u-P_{0}u)$
$\displaystyle\leq\|u\|_{a}\|u-P_{0}u\|_{a}+k^{2}\|u-T_{0}u\|\,\|P_{0}u-T_{0}u\|.$
(3.27)

Since our domain satisfies $D_{\Omega}\leq 1$, we can use the definition of
$P_{0}$ and the Friedrichs inequality to get

$\displaystyle\|u-P_{0}u\|_{a}$ $\displaystyle\leq\|u\|_{a},$
$\displaystyle\|P_{0}u-T_{0}u\|$
$\displaystyle\leq\|P_{0}u-T_{0}u\|_{a}\leq\|P_{0}(u-T_{0}u)\|_{a}\leq\|u-T_{0}u\|_{a}.$

These relations, together with (3.27) and (3.20) lead to

$\displaystyle\|u-T_{0}u\|_{a}^{2}$
$\displaystyle\leq\|u\|_{a}^{2}+k^{2}\|u-T_{0}u\|\,\|u-T_{0}u\|_{a}$
$\displaystyle\leq\left\|u\right\|_{a}^{2}+k^{2}2\Lambda\Theta^{\frac{1}{2}}(1+C_{\textup{{stab}}})\left\|T_{0}u-u\right\|_{a}^{2}.$
(3.28)

Rearranging (3.28) we arrive at

$\Bigl{(}1-2k^{2}\Lambda\Theta^{\frac{1}{2}}(1+C_{\textup{{stab}}})\Bigr{)}\|T_{0}u-u\|_{a}^{2}\leq\|u\|_{a}^{2}.$

Using assumption (3.21) we can finish the proof of (3.22). ∎

##  4 Proof of the main result

Before we can prove Theorem 3.1, we need to state and prove the following
lemma.

######  Lemma 4.1.

Suppose that the assumptions in Theorem 3.1 are satisfied. Then, for all $u\in
V_{h}$,

$c_{1}\|u\|_{a}^{2}\leq(Tu,u)_{a},$ (4.1)

and

$\|Tu\|_{a}\leq c_{2}\|u\|_{a}^{2},$ (4.2)

where $c_{1}$ and $c_{2}$ are as in Theorem 3.1.

###### Proof.

To prove (4.1), let $u\in V^{h}$. We proceed in several steps.

Step 1. This is the same procedure as in [2], but adapted for application to
(1.1) directly. We start with (3.11) and then use (3.10), (3.11), and the
definition of the $T_{i}$ projection operators to get a preliminary estimate:

$\displaystyle\bigl{(}4(1+\Lambda^{2}\Theta)\bigr{)}^{-1}\|u\|_{a}^{2}$
$\displaystyle\leq
a(Pu,u)=a\Biggl{(}\sum_{i=0}^{N}E_{i}P_{i}u,u\Biggr{)}=\sum_{i=0}^{N}a(E_{i}P_{i}u,u)$
$\displaystyle=\sum_{i=0}^{N}\Bigl{(}b(u,E_{i}P_{i}u)+k^{2}(u,E_{i}P_{i}u)\Bigr{)}$
$\displaystyle=\sum_{i=0}^{N}\Bigl{(}b(E_{i}T_{i}u,E_{i}P_{i}u)+k^{2}(u,E_{i}P_{i}u)\Bigr{)}$
$\displaystyle=\sum_{i=0}^{N}\Bigl{(}a(E_{i}T_{i}u,E_{i}P_{i}u)-k^{2}(E_{i}T_{i}u,E_{i}P_{i}u)+k^{2}(u,E_{i}P_{i}u)\Bigr{)}$
$\displaystyle=\sum_{i=0}^{N}a(E_{i}T_{i}u,E_{i}P_{i}u)-k^{2}\sum_{i=0}^{N}\bigl{(}E_{i}T_{i}u-u,E_{i}P_{i}u\bigr{)}.$

Since, by definition
$a(E_{i}T_{i}u,E_{i}P_{i}u)=a_{\Omega_{i}}(T_{i}u,P_{i}u)$ and
$a_{\Omega_{i}}(P_{i}u,v)=a(u,E_{i}v)$ for all $v\in V_{i}$, we have
$a(E_{i}T_{i}u,E_{i}P_{i}u)=a_{\Omega_{i}}(T_{i}u,P_{i}u)=a(E_{i}T_{i}u,u)$,
which, in turn, implies

$\displaystyle\bigl{(}4(1+\Lambda^{2}\Theta)\bigr{)}^{-1}\|u\|_{a}^{2}$
$\displaystyle\leq\sum_{i=0}^{N}a(E_{i}T_{i}u,u)-k^{2}\sum_{i=0}^{N}(E_{i}T_{i}u-u,E_{i}P_{i}u)$
$\displaystyle=a(Tu,u)-k^{2}\sum_{i=0}^{N}(E_{i}T_{i}u-u,E_{i}P_{i}u).$

This can be rewritten as (the sum can be split into terms related to $i=0$ and
$i\geq 1$)

$\|u\|_{a}^{2}\leq
4(1+\Lambda^{2}\Theta)\Biggl{(}a(Tu,u)-k^{2}(E_{0}T_{0}u-u,E_{0}P_{0}u)-k^{2}\sum_{i=1}^{N}(E_{i}T_{i}u-u,E_{i}P_{i}u)\Biggr{)}.$
(4.3)

Step 2. The next stage is to bound the second and third terms inside the
bracket of the preliminary estimate (4.3). It is at this stage that we take
advantage of the improved estimates obtained in Lemmas 3.11 and 3.12. We start
with the second term in (4.3), where $E_{0}$ should be considered as a natural
embedding like the identity operator ($E_{i}T_{0}u=T_{0}u$ and
$E_{i}P_{0}u=P_{0}u$). We also use (3.20) and (3.22) to get

$\displaystyle-k^{2}(E_{0}T_{0}u-u,E_{0}P_{0}u)$ $\displaystyle\leq
k^{2}\|E_{0}T_{0}u-u\|\,\|E_{0}P_{0}u\|=k^{2}\|T_{0}u-u\|\|P_{0}u\|$
$\displaystyle\leq k^{2}\|T_{0}u-u\|\,\|P_{0}u\|_{a}\leq
k^{2}2\Lambda\Theta^{\frac{1}{2}}(1+C_{\textup{{stab}}})\|T_{0}u-u\|_{a}\|P_{0}u\|_{a}$
$\displaystyle\leq
k^{2}2\sqrt{2}\Lambda\Theta^{\frac{1}{2}}(1+C_{\textup{{stab}}})\|u\|_{a}\|P_{0}u\|_{a}.$

With this the second term in (4.3) can be estimated as follows, where $s$ is
as in Theorem 3.1,

$-4(1+\Lambda^{2}\Theta)k^{2}\bigl{(}E_{0}T_{0}u-u,E_{0}P_{0}u\bigr{)}\leq\frac{s}{2\Lambda^{\frac{1}{2}}}\|u\|_{a}\|P_{0}u\|_{a}\leq\frac{s}{\sqrt{2}(\Lambda+1)^{\frac{1}{2}}}\|u\|_{a}\|P_{0}u\|_{a}.$
(4.4)

Let us now consider the sum in (4.3). Note that, for each
$i\in\\{1,\ldots,N\\}$, the operator $E_{i}$ is just an extension by $0$
outside the domain ($a(u,E_{i}v)=a_{\Omega_{i}}(u,v)$ for $u\in V^{h},v\in
V_{i}$). Using Lemma 3.10 we obtain

$\displaystyle-k^{2}(E_{i}T_{i}u-u,E_{i}P_{i}u)\leq\big{|}k^{2}(E_{i}T_{i}u-u,P_{i}u)_{\Omega_{i}}\big{|}\leq
k^{2}\|E_{i}T_{i}u-u\|_{\Omega_{i}}\|P_{i}u\|_{\Omega_{i}}$ $\displaystyle\leq
k^{2}\|E_{i}T_{i}u-u\|_{\Omega_{i}}\|P_{i}u\|_{\Omega_{i}}\leq
k^{2}\bigl{(}\|T_{i}u\|_{\Omega_{i}}+\|u\|_{\Omega_{i}}\bigr{)}\|P_{i}u\|_{\Omega_{i}}$
$\displaystyle\leq
k^{2}\biggl{(}\frac{H}{\sqrt{2}\,}\|T_{i}u\|_{a_{\Omega_{i}}}+\|u\|_{\Omega_{i}}\biggr{)}\|P_{i}u\|_{\Omega_{i}}\leq
k^{2}\Bigl{(}H\sqrt{2}\|u\|_{a_{\Omega_{i}}}+kH\|u\|+\|u\|_{\Omega_{i}}\Bigr{)}\|P_{i}u\|_{\Omega_{i}}$
$\displaystyle\leq
Hk^{2}\Bigl{(}H\sqrt{2}\|u\|_{a_{\Omega_{i}}}+3\|u\|_{\Omega_{i}}\Bigr{)}\|P_{i}u\|_{a_{\Omega_{i}}}\leq
Hk^{2}\biggl{(}H\sqrt{2}\|u\|_{a_{\Omega_{i}}}+\frac{3H}{\sqrt{2}\,}\|u\|_{a_{\Omega_{i}}}\biggr{)}\|P_{i}u\|_{a_{\Omega_{i}}}$
$\displaystyle<4H^{2}k^{2}\|u\|_{a_{\Omega_{i}}}\|P_{i}u\|_{a_{\Omega_{i}}}.$

Taking the sum over $i$ and applying the Cauchy–Schwarz inequality and the
overlap property,
$\sum_{i=1}^{N}\|u\|_{a_{\Omega_{i}}}^{2}\leq\Lambda\|u\|_{a}^{2}$, we obtain

$-k^{2}\sum_{i=1}^{N}(E_{i}T_{i}u-u,E_{i}P_{i}u)\leq\sum_{i=1}^{N}4H^{2}k^{2}\|u\|_{a_{\Omega_{i}}}\|P_{i}u\|_{a_{\Omega_{i}}}\\\
\leq
4H^{2}k^{2}\Biggl{(}\sum_{i=1}^{N}\|u\|_{a_{\Omega_{i}}}^{2}\Biggr{)}^{\frac{1}{2}}\Biggl{(}\sum_{i=1}^{N}\|P_{i}u\|_{a_{\Omega_{i}}}^{2}\Biggr{)}^{\frac{1}{2}}\leq
4H^{2}k^{2}\Lambda^{\frac{1}{2}}\|u\|_{a}\Biggl{(}\sum_{i=1}^{N}\|P_{i}u\|_{a_{\Omega_{i}}}^{2}\Biggr{)}^{\frac{1}{2}}.$

Using assumption (3.1) and $t$ as in Theorem 3.1 we arrive at

$\displaystyle-4(1+\Lambda^{2}\Theta)k^{2}\sum_{i=1}^{N}\bigl{(}E_{i}T_{i}u-u,E_{i}P_{i}u\bigr{)}\leq
16(1+\Lambda^{2}\Theta)H^{2}k^{2}\Lambda^{\frac{1}{2}}\|u\|_{a}\Biggl{(}\sum_{i=1}^{N}\|P_{i}u\|_{a_{\Omega_{i}}}^{2}\Biggr{)}^{\frac{1}{2}}\quad$
$\displaystyle\leq\frac{t}{2\Lambda^{\frac{1}{2}}}\|u\|_{a}\Biggl{(}\sum_{i=1}^{N}\|P_{i}u\|_{a_{\Omega_{i}}}^{2}\Biggr{)}^{\frac{1}{2}}\leq\frac{t}{\sqrt{2}(\Lambda+1)^{\frac{1}{2}}}\|u\|_{a}\Biggl{(}\sum_{i=1}^{N}\|P_{i}u\|_{a_{\Omega_{i}}}^{2}\Biggr{)}^{\frac{1}{2}}.$
(4.5)

Step 3. We are now ready to combine (4.4) and (4.5) in order to get an
estimate for the second and third terms in (4.3) in terms of $\|u\|^{2}_{a}$.
This again follows the same procedure as used in [2] with no further
refinements (note that we also use the well-known inequality $(a+b)^{2}\leq
2a^{2}+2b^{2}\;\Rightarrow\;a+b\leq\sqrt{2a^{2}+2b^{2}}$),

$\displaystyle-4(1+\Lambda^{2}\Theta)\Biggl{(}k^{2}(E_{0}T_{0}u-u,E_{0}P_{0}u)+k^{2}\sum_{i=1}^{N}(E_{i}T_{i}u-u,E_{i}P_{i}u)\Biggr{)}$
$\displaystyle\leq\frac{s}{\sqrt{2}(\Lambda+1)^{\frac{1}{2}}}\|u\|_{a}\|P_{0}u\|_{a}+\frac{t}{\sqrt{2}(\Lambda+1)^{\frac{1}{2}}}\|u\|_{a}\Biggl{(}\sum_{i=1}^{N}\|P_{i}u\|_{a_{\Omega_{i}}}^{2}\Biggr{)}^{\frac{1}{2}}$
$\displaystyle\leq\frac{\max\\{s,t\\}}{\sqrt{2}(\Lambda+1)^{1/2}}\|u\|_{a}\Biggl{[}\|P_{0}u\|_{a}+\Biggl{(}\sum_{i=1}^{N}\|P_{i}u\|_{a_{\Omega_{i}}}^{2}\Biggr{)}^{\frac{1}{2}}\Biggr{]}$
$\displaystyle\leq\frac{\max\\{s,t\\}}{(\Lambda+1)^{1/2}}\|u\|_{a}\Biggl{(}\sum_{i=0}^{N}\|P_{i}u\|_{a_{\Omega_{i}}}^{2}\Biggr{)}^{\frac{1}{2}}$
$\displaystyle\leq\frac{\max\\{s,t\\}}{(\Lambda+1)^{1/2}}\|u\|_{a}\bigl{(}\|u\|_{a}^{2}(\Lambda+1)\bigr{)}^{\frac{1}{2}}=\max\\{s,t\\}\|u\|^{2}_{a}.$
(4.6)

In the last inequality we used (3.12).

Step 4. Although no further improvements to the results are obtained at this
stage, (3.22) is a direct result of earlier enhancements to the estimations
being found. Following the same procedure as used in [2], we exploit the
estimate (4.6) and proceed to the final result. In a first instance we use
(4.6) together with (4.3) to get

$\|u\|^{2}_{a}\leq 4(1+\Lambda^{2}\Theta)a(Tu,u)+\max\\{s,t\\}\|u\|^{2}_{a},$

which leads to the desired estimate

$c_{1}\|u\|^{2}_{a}=\bigl{(}4(1+\Lambda^{2}\Theta)\bigr{)}^{-1}\bigl{(}1-\max\\{s,t\\}\bigr{)}\|u\|^{2}_{a}\leq
a(Tu,u).$

We now come to the proof of (4.2). We start with the following relation, where
we use $E_{0}T_{0}u=T_{0}u$,

$\|Tu\|_{a}^{2}=\Bigg{\|}T_{0}u+\sum_{i=1}^{N}E_{i}T_{i}u\Bigg{\|}_{a}^{2}\leq
2\|T_{0}u\|_{a}^{2}+2\Bigg{\|}\sum_{i=1}^{N}E_{i}T_{i}u\Bigg{\|}_{a}^{2}.$
(4.7)

Using the Cauchy–Schwarz inequality and (3.22) we find an upper bound for the
first term on the right-hand side of (4.7):

$\displaystyle\|T_{0}u\|_{a}^{2}$
$\displaystyle=a(T_{0}u,T_{0}u)=a(T_{0}u,T_{0}u)-a(u,T_{0}u)+a(u,T_{0}u)$
$\displaystyle\leq
a(T_{0}u-u,T_{0}u)+a(u,T_{0}u)\leq\|T_{0}u-u\|_{a}\|T_{0}u\|_{a}+\|u\|_{a}\|T_{0}u\|_{a}$
$\displaystyle\leq 3\|u\|_{a}\|T_{0}u\|_{a},$

which yields

$\|T_{0}u\|_{a}^{2}\leq 9\|u\|_{a}^{2}.$ (4.8)

For the second term on the right-hand side of (4.7) we use Lemma 3.10 to get

$\displaystyle\Bigg{\|}\sum_{i=1}^{N}E_{i}T_{i}u\Bigg{\|}_{a}^{2}$
$\displaystyle\leq\Lambda\sum_{i=1}^{N}\|T_{i}u\|_{a_{\Omega_{i}}}^{2}\leq\Lambda\sum_{i=1}^{N}\Bigl{(}2\|u\|_{a_{\Omega_{j}}}+k\sqrt{2}\|u\|_{\Omega_{j}}\Bigr{)}^{2}$
$\displaystyle\leq\Lambda\sum_{i=1}^{N}\Bigl{(}2\|u\|_{a_{\Omega_{j}}}+kH\|u\|_{a_{\Omega_{j}}}\Bigr{)}^{2}$
$\displaystyle\leq 9\Lambda\sum_{i=1}^{N}\|u\|_{a_{\Omega_{j}}}^{2}\leq
9\Lambda^{2}\|u\|_{a}^{2}.$ (4.9)

Combining (4.8) and (4.9) with (4.7) we arrive at

$\|Tu\|_{a}\leq
2(9\|u\|_{a}^{2})+2(9\Lambda^{2}\|u\|_{a}^{2})=(18+18\Lambda^{2})\|u\|_{a}^{2},$
(4.10)

which finishes the proof of (4.2). ∎

We are now ready to prove Theorem 3.1.

###### Proof of Theorem 3.1.

Using (4.1) with (2.16) we obtain

$c_{1}\|u\|_{a}^{2}\leq
a(Tu,u)=\langle\mathbf{M}_{AS,2}^{-1}\mathbf{B}\mathbf{u},\mathbf{u}\rangle_{\mathbf{A}},$

which can be written as

$c_{1}\leq\frac{\langle\mathbf{M}_{AS,2}^{-1}\mathbf{B}\mathbf{u},\mathbf{u}\rangle_{\mathbf{A}}}{\|\mathbf{u}\|_{\mathbf{A}}^{2}}.$

In a similar way we use (4.2), again with (2.16), to get

$\|\mathbf{M}_{AS,2}^{-1}\mathbf{B}\mathbf{u}\|_{\mathbf{A}}^{2}=\|Tu\|_{a}^{2}\leq
c_{2}\|u\|_{a}^{2}=c_{2}\|\mathbf{u}\|^{2}_{\mathbf{A}},$

which implies

$\|\mathbf{M}_{AS,2}^{-1}\mathbf{B}\|_{\mathbf{A}}^{2}\leq c_{2}.$

Now the result follows directly from Elman theory [3]. ∎

## References

  * [1] N. Bootland, V. Dolean, I. G. Graham, C. Ma, R. Scheichl, GenEO coarse spaces for heterogeneous indefinite elliptic problems. In Domain Decomposition Methods in Science and Engineering XXVI, Lect. Notes Comput. Sci. Eng., vol. 145, Springer, Cham, Switzerland, 2022, pp. 117–125.
  * [2] N. Bootland, V. Dolean, I. G. Graham, C. Ma, R. Scheichl, Overlapping Schwarz methods with GenEO coarse spaces for indefinite and non-self-adjoint problems. IMA J. Numer. Anal. 43, 1899–1936, 2023.
  * [3] S. C. Eisenstat, H. C. Elman and M. H. Schultz, Variational iterative methods for nonsymmetric systems of linear equations, SIAM J. Numer. Anal. 20, 345–357, 1983.
  * [4] G. Leoni, A First Course in Sobolev Spaces, Grad. Stud. Math. 181, Amer. Math. Soc., Providence, RI, 2017.
  * [5] A. H. Schatz and J. P. Wang, Some new error estimates for Ritz–Galerkin methods with minimal regularity assumptions, Math. Comput., 213(65), 19–27, 1996.
  * [6] N. Spillane, F. Nataf, V. Dolean, P. Hauret, C. Pechstein and R. Scheichl, Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps, Numer. Math. 4(126), 741–740, 2014.
  * [7] A. Toselli and O. Widlund, Domain Decomposition Methods — Algorithms and Theory, Springer Ser. Comput. Math. 34, Springer-Verlag, Berlin, 2005.


# Anomalous Gilbert Damping and Duffing Features of the SFS $\varphi_{0}$
Josephson Junction

Yu. M. Shukrinov1,2 I. R. Rahmonov1,3 A. Janalizadeh4 M. R. Kolahchi4 1 BLTP,
JINR, Dubna, Moscow Region, 141980, Russia
2 Dubna State University, Dubna, 141980, Russia
3 Umarov Physical Technical Institute, TAS, Dushanbe 734063, Tajikistan
4 Department of Physics, Institute for Advanced Studies in Basic Sciences,
P.O. Box 45137-66731, Zanjan, Iran

(August 27, 2024)

###### Abstract

We demonstrate unusual features of phase dynamics, IV-characteristics and
magnetization dynamics of the $\varphi_{0}$ Josephson junction at small values
of spin-orbit interaction, ratio of Josephson to magnetic energy and Gilbert
damping. In particular, an anomalous shift of the ferromagnetic resonance
frequency with an increase of Gilbert damping is found. The ferromagnetic
resonance curves show the Duffing oscillator behaviour, reflecting the
nonlinear nature of Landau-Lifshitz-Gilbert (LLG) equation. Based on the
numerical analysis of each term in LLG equation we obtained an approximated
equation demonstrated both damping effect and Duffing oscillator features. The
resulting Duffing equation incorporates the Gilbert damping in a special way
across the dissipative term and the restoring force. A resonance method for
the determination of spin-orbit interaction in noncentrosymmetric materials
which play the role of barrier in $\varphi_{0}$ junctions is proposed.

## Introduction.

The Josephson junctions (JJ) with the current-phase relation
$I=I_{c}\sin(\varphi-\varphi_{0})$, where the phase shift $\varphi_{0}$ is
proportional to the magnetic moment of ferromagnetic layer determined by the
parameter of spin-orbit interaction, demonstrate a number of unique features
important for superconducting spintronics, and modern information technology
linder15 ; shukrinov-ufn21 ; bobkova20 ; mazanik-pra20 ; shukrinov-apl17 ;
nashaat-pepan20 . The phase shift allows one to manipulate the internal
magnetic moment using the Josephson current, and the reverse phenomenon which
leads to the appearance of the DC component in the superconducting current
buzdin-prl08 ; konschelle-prl09 ; shukrinov-prb19 .

Interactive fields can bring nonlinear phenomena of both classical, and
quantum nature. A basic example is the magnons strongly interacting with
microwave photons zha14 . As a result we could name Bose-Einstein condensation
of such quasiparticles, i.e. magnons ser14 ; dem08 , and synchronization of
spin torque nano-oscillators as they coherently emit microwave signals in
response to d.c. current sheh05 . It is interesting that (semi)classical
anharmonic effects in the magnetodynamics described by the Landau-Lifshitz-
Gilbert (LLG) model in thin films or heterostructures shen20 ; nik20 , and the
quantum anharmonicity in the cavity mangnonics elyasi20  can well be modeled
by so simple a nonlinear oscillator as Duffing. The corresponding Duffing
equation contains a cubic term and describes the oscillations of the various
nonlinear systems moon15 .

Despite the fact that nonlinear features of LLG are studied often during a
long time and in different systems, manifestation of the Duffing oscillator
behavior in the framework of this equation is still not completely studied.
Closer to our present investigation, in the study of the dynamics of
antiferromagnetic bimeron under an alternating current, Duffing equation forms
a good model, and this has applications in weak signal detection shen20 ;
wan99 ; alm07 . As another application with Duffing oscillator at work, we can
mention the ultra thin Co20Fe60B20 layer, and its large angle magnetization
precession under microwave voltage. There are also ‘foldover’ features,
characteristic of the Duffing spring, in the magnetization dynamics of the
Co/Ni multilayer excited by a microwave current nik20 ; nay79 ; chen09 . But
nonlinear features of $\varphi_{0}$ Josephson junctions have not been
carefully studied yet. In this Letter, we show that the Duffing oscillator
helps in the understanding of the nonlinear features of $\varphi_{0}$
Josephson junctions at small values of system parameters.

Coupling of superconducting current and magnetization and its manifestation in
the IV-characteristics and magnetization dynamics opens the door for the
resonance method determination of spin-orbit intensity in noncentrosymmetric
materials playing the role of barrier in $\varphi_{0}$ junctions. As it is
well known, the spin-orbit interaction plays an important role in modern
physics, so any novel method for its determination in real materials would be
very important. There are a series of recent experiments demonstrating the
modification of Gilbert damping by the superconducting correlations (see
Ref.silaev2020  and citations therein). In particular, the pronounced peaks in
the temperature dependence of Gilbert damping have been observed for the
ferromagnetic insulator/superconductor multilayers yao2018  which might be
explained by the presence of spin relaxation mechanisms like the spin-orbit
scattering silaev2020 . Here, we use the noncentrosymmetric ferromagnetic
material as a weak link in $\varphi_{0}$ junctions. The suitable candidates
may be MnSi or FeGe, where the lack of inversion center comes from the
crystalline structure konschelle-prl09 .

The Gilbert damping determines the magnetization dynamics in ferromagnetic
materials but its origin is not well understood yet. Effect of nonlinearity on
damping in the system is very important for application of these materials in
fast switching spintronics devices. Our study clarifies such effects. In
Ref.zhao16  the authors discuss the experimental study of temperature-
dependent Gilbert damping in permalloy (Py) thin films of varying thicknesses
by ferromagnetic resonance, and provide an important insight into the physical
origin of the Gilbert damping in ultrathin magnetic films.

In this Letter we demonstrate an anomalous dependence of the ferromagnetic
resonance frequency with an increase of the Gilbert damping. We find that the
resonance curves demonstrate features of Duffing oscillator, reflecting the
nonlinear nature of LLG equation. The damped precession of the magnetic moment
is dynamically driven by the Josephson supercurrent, and the resonance
behavior is given by the dynamics of the Duffing spring. The resonance methods
for the determination of spin-orbit interaction in the $\varphi_{0}$ junction
are proposed.

## Model and Methods.

In the considered SFS $\varphi_{0}$ junction (see Fig.1) the superconducting
phase difference $\varphi$ and magnetization $\mathbf{M}$ of the F layer are
two coupled dynamical variables.

Figure 1: Schematic view of SFS $\varphi_{0}$ Josephson junction. The external
current applied along $x$ direction, ferromagnetic easy axis is along $z$
direction.

Based on the LLG equation for the magnetic moment $\mathbf{M}$ with effective
magnetic field ${\bf H_{eff}}$, resistively capacitively shunted junction
(RCSJ) model, and Josephson relation for the phase difference $\varphi$, we
describe dynamics of the SFS $\varphi_{0}$ junction by the system of equations
in normalized variables

$\displaystyle\displaystyle\frac{d{\bf m}}{dt}$ $\displaystyle=$
$\displaystyle\omega_{F}{\bf h_{eff}}\times{\bf m}+\alpha\bigg{(}{\bf
m}\times\frac{d{\bf m}}{dt}\bigg{)},\vspace{0.1 cm}$
$\displaystyle\displaystyle{\bf h_{eff}}$ $\displaystyle=$ $\displaystyle
Gr\sin(\varphi-rm_{y}){\bf\widehat{y}}+m_{z}{\bf\widehat{z}},\vspace{0.1 cm}$
(1) $\displaystyle\displaystyle\frac{dV}{dt}$ $\displaystyle=$
$\displaystyle\frac{1}{\beta_{c}}[I-V+r\frac{dm_{y}}{dt}-\sin(\varphi-
rm_{y})],\vspace{0.1 cm}$ $\displaystyle\displaystyle\frac{d\varphi}{dt}$
$\displaystyle=$ $\displaystyle V,\vspace{0.1 cm}$

where ${\bf m}$ is vector of magnetization with components $m_{x,y,z}$,
normalized to the $M_{0}=\|{\bf M}\|$ and and satisfying the constraint
$\sum_{i=x,y,z}m_{i}^{2}(t)=1$, $\omega_{F}=\Omega_{F}/\omega_{c}$,
$\Omega_{F}=\gamma K/\nu$ is ferromagnetic resonance frequency, $\gamma$ is
the gyromagnetic ratio, $K$ is an anisotropic constant, $\nu$ is the volume of
the ferromagnetic F layer, $\alpha$ is the phenomenological damping constant
(Gilbert damping), ${\bf h_{eff}}$ is the vector of effective magnetic field,
normalized to the $K/M_{0}$ (${\bf h_{eff}}={\bf H_{eff}}M_{0}/K$),
$G=E_{J}/(K\nu)$ relation of Josephson energy to magnetic one, $r$ is a
parameter of spin-orbit coupling, $\varphi$ is phase difference of JJ, $V$ is
voltage normalized to the $V_{c}=I_{c}R$, $I_{c}$ critical current of JJ, $R$
resistance of JJ, $\beta_{c}=2eI_{c}CR^{2}/\hbar$ is McCumber parameter, $C$
is capacitance of JJ, $I$ is bias current normalized to the $I_{c}$. In this
system of equation time $t$ is normalized to the $\omega_{c}^{-1}$, where
$\omega_{c}=2eI_{c}R/\hbar$ is characteristic frequency. In the chosen
normalization, the average voltage corresponds to the Josephson frequency
$\omega_{J}$.

## Ferromagnetic resonance in $\varphi_{0}$ junction.

The ferromagnetic resonance features are demonstrated by average voltage
dependence of the maximal amplitude of the $m_{y}$ component ($m_{y}^{max}$),
taken at each value of bias current. To stress novelty and importance of our
finding, we first present the analytical results for average voltage
dependence of $m_{y}^{max}$ along IV-characteristics in the ferromagnetic
resonance region. As it was discussed in Refs.konschelle-prl09 ; shukrinov-
pepan20 ; shukrinov-ltp20 , in case $Gr\ll 1$, $m_{z}\approx 1$, and
neglecting quadratic terms $m_{x}$ and $m_{y}$, we get

$\left\\{\begin{array}[]{ll}\displaystyle\dot{m}_{x}=\xi[-m_{y}+Gr\sin\omega_{J}t-\alpha
m_{x}]\vspace{0.2 cm}\\\ \displaystyle\dot{m}_{y}=\xi[m_{x}-\alpha
m_{y}],\end{array}\right.$ (2)

where $\xi=\omega_{F}/(1+\alpha^{2})$. This system of equations can be written
as the second order differential equation with respect to the $m_{y}$

$\displaystyle\ddot{m}_{y}=-2\alpha\xi\dot{m}_{y}-\xi^{2}(1+\alpha^{2})m_{y}+\xi^{2}Gr\sin\omega_{J}t.$
(3)

Corresponding solution for $m_{y}$ has the form

$m_{y}(t)=\frac{\omega_{+}-\omega_{-}}{r}\sin\omega_{J}t-\frac{\alpha_{+}+\alpha_{-}}{r}\cos\omega_{J}t,$
(4)

where

$\omega_{\pm}=\frac{Gr^{2}\omega_{F}}{2}\frac{\omega_{J}\pm\omega_{F}}{((\omega_{J}\pm\omega_{F})^{2}+(\alpha\omega_{J})^{2})},$
(5)

and

$\alpha_{\pm}=\frac{Gr^{2}\omega_{F}}{2}\frac{\alpha\omega_{J}}{((\omega_{J}\pm\omega_{F})^{2}+(\alpha\omega_{J})^{2})}.$
(6)

So, $m_{y}$ demonstrates resonance with dissipation when Josephson frequency
is approaching the ferromagnetic one ($\omega_{J}\rightarrow\omega_{F}$). The
maximal amplitude $m_{y}^{max}$ as a function of voltage (i.e., Josephson
frequency $\omega_{J}$) at different $\alpha$, calculated using (4), is
presented in Fig.2 (a). We see the usual characteristic variation of the
resonance curve with an increase in dissipation parameter when the maximal
amplitude and position of resonance pick corresponds to the damped resonance.
We note that the analytical result (4) were obtained in the case $Gr\ll 1$.

Figure 2: (a) Analytical results for maximal amplitude $m_{y}^{max}$ in the
ferromagnetic resonance region for different $\alpha$; (b) Numerical results
for maximal amplitude of magnetization $m_{y}-$component at each values of
bias current and voltage along IV-characteristics of the $\varphi_{0}$
junction in the ferromagnetic resonance region for various $\alpha$. Inset
shows the manifestation of the resonance subharmonics. Parameters are:
$\beta_{c}=25$, G=0.05, r=0.05,$\omega_{F}=0.5$.

Presented in Fig.2(b) results of numerical simulations $m_{y}^{max}(V)$
dependence at different values of dissipation parameter $\alpha$ demonstrate
the essential differences with the results followed from the analytical
consideration (4). We note also that the strong coupling of the
superconducting phase difference $\varphi$ and magnetization $\mathbf{M}$ of
the F layer manifests itself by appearance of subharmonics of the resonance at
$\omega=1/2,1/3,1/4$ demonstrated in the inset to Fig.2(b).

We stress two important features followed from the presented results. First,
the ferromagnetic resonance curves show the foldover effect, i.e., the
features of Duffing oscillator. Different from a linear oscillator, the
nonlinear Duffing demonstrates a bistability under external periodic force
kovacic11 . Second, the ferromagnetic resonance curves demonstrate an unusual
dependence of the resonance frequency as a function of Gilbert damping
$\alpha$. As shown in Fig. 3(a), an increase in damping leads to a nonuniform
change in the resonant frequency, i.e., with an increase in damping the
resonance maximum shifts to $\omega_{F}$ at small $\alpha$, but then moves to
the opposite side, demonstrating the usual damped resonance. So, with an
increase in $\alpha$, unusual dependence of the resonance voltage transforms
to the usual one. For the parameters chosen, the critical value of this
transformation is around $\alpha=0.02-0.03$. We call this unusual behaviour of
the resonance maximum of $m_{y}^{max}$ as an “$\alpha$-effect”.

Figure 3: (a) $\alpha$-dependence of the resonance curve $m_{y}^{max}(V)$ peak
presented in Fig.2 in the damping parameter interval [0.006 – 0.2]. Dashed
line indicates ferromagnetic resonance position; (b) Comparison of the
resonance curves $m_{y}^{max}(V)$ calculated by full LLG equation (Model and
Methods.) and the approximate equation (The generalized Duffing equation for
$\varphi_{0}$ junction.).

Both the $\alpha-$effect and Duffing features in our system appear due to the
nonlinear features of the system dynamics at small $G,r,\alpha\ll 1$. To prove
it, we have carried out the numerical analysis of each term of LLG full
equation (first two equations in (Model and Methods.)) for the set of model
parameters $G=0.05$, $r=0.05$ $\alpha=0.005$. After neglecting the terms of
order $10^{-6}$, we have

$\displaystyle\frac{\dot{m}_{x}}{\xi}$ $\displaystyle=$ $\displaystyle-
m_{y}m_{z}+Grm_{z}\sin(\varphi-rm_{y})-\alpha m_{x}m_{z}^{2},$
$\displaystyle\frac{\dot{m}_{y}}{\xi}$ $\displaystyle=$ $\displaystyle
m_{x}m_{z}-\alpha m_{y}m_{z}^{2},$ (7) $\displaystyle\frac{\dot{m}_{z}}{\xi}$
$\displaystyle=$ $\displaystyle-Grm_{x}\sin(\varphi-rm_{y})+\alpha
m_{z}(m_{x}^{2}+m_{y}^{2}),$

In this approximation we observe both the “$\alpha$–effect” and Duffing
oscillator features. Neglecting here the last term $\alpha
m_{z}(m_{x}^{2}+m_{y}^{2})$ in third equation for $\dot{m}_{z}$, which is
order of $10^{-4}$, leads to the losing of the Duffing oscillator features,
but still keeps alpha-effect. We note that equation (Ferromagnetic resonance
in $\varphi_{0}$ junction.) keeps the time invariance of the magnetic moment,
so that term plays an important role for manifestation of Duffing oscillator
features by LLG equation.

## The generalized Duffing equation for $\varphi_{0}$ junction.

The LLG is a nonlinear equation and in case of simple effective field it can
be transformed to the Duffing equation moon15 ; shen20 . Such transformation
was used in Ref.moon15  to demonstrate the nonlinear dynamics of the magnetic
vortex state in a circular nanodisk under a perpendicular alternating magnetic
field that excites the radial modes of the magnetic resonance. They showed
Duffing-type nonlinear resonance and built a theoretical model corresponding
to the Duffing oscillator from the LLG equation to explore the physics of the
magnetic vortex core polarity switching for magnetic storage devices.

The approximated LLG system of equations (Ferromagnetic resonance in
$\varphi_{0}$ junction.) demonstrates both $\alpha$-effect and features of
Duffing oscillator. As demonstrated in the Supplemental Materials supplement ,
the generalized Duffing equation for the $\varphi_{0}$ junction,

$\displaystyle\ddot{m}_{y}+2\xi\alpha\dot{m}_{y}+\xi^{2}(1+\alpha^{2})m_{y}$
$\displaystyle-\xi^{2}(1+\alpha^{2})m_{y}^{3}$ $\displaystyle=$
$\displaystyle\xi^{2}Gr\sin\omega_{J}t.$ (8)

can be obtained directly from the LLG system of equations.

As we see, for small enough $G$ and $r$, it is only the dimensionless damping
parameter $\alpha$ in LLG that plays a role in the dynamics of the system. We
can think of a harmonic spring with a constant that is hardened or softened by
the nonlinear term. For a usual Duffing spring, with independent coefficients
of the various terms, the resonance peak relative to the harmonic (linear)
resonant frequency folds over to the smaller (softening) or larger (hardening)
frequencies. In the frequency response, the interplay of the specific
dependence of each coefficient on $\alpha$ plays an important role and as
Fig.3(a) shows, there is a particular $\alpha$ that brings the resonant
frequency closest to ferromagnetic resonance.

Simulations of the $m_{y}$ dynamics in the framework of Duffing equation can
explain observed foldover effect in the frequency dependence of $m_{y}^{max}$.
Comparison the results followed from analytical approximate equation (The
generalized Duffing equation for $\varphi_{0}$ junction.) and results of full
equation (Model and Methods.) for maximal amplitude of $m_{y}^{max}$ in the
ferromagnetic resonance region is presented in Fig.3(b). So, the magnetization
dynamics in the SFS $\varphi_{0}$-junction due to the voltage oscillations can
effectively be described by a scalar Duffing oscillator, synchronizing the
precession of the magnetic moment with the Josephson oscillations.

## Effect of spin-orbit interactions.

As we mentioned above, the spin-orbit interaction plays an important role in
different fields of modern physics. Here we have suggested a novel method for
its determination in real noncentrosymmetric ferromagnetic materials like MnSi
or FeGe, where the lack of inversion center comes from the crystalline
structure Ref.konschelle-prl09  and which play role a weak link in
$\varphi_{0}$ junctions. Based on the obtained results, presented in Fig.4, we
propose different versions of the resonance method for the determination of
spin-orbit interaction in these materials.

Figure 4: (a) Voltage dependence of $m_{y}^{max}$ in the ferromagnetic
resonance region at different values of spin-orbit interaction based on (Model
and Methods.) at $G=0.05$, $\alpha=0.01$. Inset enlarges the main harmonic;
(b) The same as in (a) for $\alpha=0.1$; (c) Shift of $m_{y}^{max}$ peak as a
function of spin-orbit interaction at two values of Gilbert damping; (d)
r-dependence of the main harmonic and subharmonics peaks in case (a); (e) The
same as in (d) for the case (b).

Particularly, in Fig.4(a) we present the simulation results of maximal
amplitude $m_{y}^{max}$ based on (Model and Methods.) at $G=0.05$,
$\alpha=0.01$ at different values of spin-orbit parameter $r$ in the
ferromagnetic resonance region. This case corresponds to the nonlinear
approximation leading to the Duffing equation (The generalized Duffing
equation for $\varphi_{0}$ junction.). The same characteristics calculated by
equation (Model and Methods.) for larger value $\alpha=0.1$, i.e.
corresponding to the linear approximation (3) are presented in Fig.4(b). As it
was expected, in case $\alpha=0.01$ the foldover effect is more distinct.

In Fig.4(c) the $r$-dependence of the resonance peak position, obtained from
the simulation results of full equation at $\alpha=0.01$ and $\alpha=0.1$ for
the same set of model and simulation parameters is demonstrated. We stress
here that nonlinear features of LLG equation leading to the Duffing’s shift of
the $m_{y}^{max}$ peak of main harmonic with r presented in Fig.4(c) show the
manifestation of nonlinearity.

Despite the noted differences between results for $\alpha=0.01$ and
$\alpha=0.1$ , we see in both cases a monotonic linear increase of
$m_{y}^{max}$ peak of main harmonic and subharmonics with $r$ demonstrated in
Fig.4(d) and Fig.4(e). Such linear dependence can be noted from Eq. (6) of
Ref. shen20 , but the authors did not discuss it. This dependence might serve
as a calibrated curve for spin-orbit interaction intensity, thus creating the
resonance methods for $r$ determination.

## Conclusions.

Based on the reported features of the $\varphi_{0}$ Josephson junction at
small values of spin-orbit interaction, ratio of Josephson to magnetic energy
and Gilbert damping, we have demonstrated that the coupled superconducting
current and the magnetic moments in the $\varphi_{0}$-junction result in the
current phase relation intertwining with the ferromagnetic LLG dynamics. The
ferromagnetic resonance clearly shows this interplay. In particular, an
anomalous shift of the ferromagnetic resonance frequency with an increase of
Gilbert damping is found. The ferromagnetic resonance curves demonstrate
features of Duffing oscillator, reflecting the nonlinear nature of LLG
equation. The obtained approximated equation demonstrates both damping effect
and Duffing oscillator features. We have shown that due to the nonlinearity,
as modeled by the generalized Duffing equation, the parameters of the system
can compensate each other resulting in unusual response. The position of the
maximum can shift towards and then away from the expected resonant frequency,
as the damping is decreased. There are also foldover effects that was
explained by the proposed model. A resonance method for the determination of
spin-orbit interaction in noncentrosymmetric materials which play the role of
barrier in $\varphi_{0}$ junctions was proposed.

The experimental testing of our results would involve SFS structures with
ferromagnetic material having enough small value of Gilbert damping. Potential
candidate for experimental realization could be ferromagnetic metals or
insulators which have small values of damping parameter ($\alpha\sim
10^{-3}-10^{-4}$). In Ref.schoen16  the authors report on a binary alloy of
cobalt and iron that exhibits a damping parameter approaching $10^{-4}$, which
is comparable to values reported only for ferrimagnetic insulators onbasli14 ;
kelly13 . Using superconductor-ferromagnetic insulator-superconductor on a 3D
topological insulator might be a way to have strong spin-orbit coupling needed
for $\varphi_{0}$ JJ and small Gilbert dissipation for $\alpha$-effect
bobkova20 . We note in this connection that the yttrium iron garnet YIG is
especially interesting because of its small Gilbert damping ($\alpha\sim
10^{-5}$). The interaction between the Josephson current and magnetization is
determined by the ratio of the Josephson to the magnetic anisotropy energy
$G=E_{J}/(K\nu)$ and spin-orbit interaction $r$. The value of the Rashba-type
parameter $r$ in a permalloy doped with $Pt$ hrabec16  and in the ferromagnets
without inversion symmetry, like MnSi or FeGe, is usually estimated to be in
the range $0.1-1$. The value of the product $Gr$ in the material with weak
magnetic anisotropy $K\sim 4\times 10^{-5}KA^{-3}$ rusanov04 , and a junction
with a relatively high critical current density of $(3\times 10^{5}-5\times
10^{6})A/cm^{2}$ robinson12  is in the range $1-100$. It gives the set of
ferromagnetic layer parameters and junction geometry that make it possible to
reach the values used in our numerical calculations for the possible
experimental observation of the predicted effect.

Numerical simulations were funded by the project 18-71-10095 of the Russian
Scientific Fund. A.J. and M.R.K. are grateful to IASBS for financial support.

## References

  * (1) Jacob Linder and W. A. Jason Robinson, Nature Physics 11, 307 (2015).
  * (2) Yu. M. Shukrinov, Accepted for UFN. DOI:https://doi.org/10.3367/UFNe.2020.11.038894
  * (3) A.A. Mazanik, I.R. Rahmonov, A.E. Botha, and Yu.M. Shukrinov, Phys. Rev. Applied 14, 014003 (2020).
  * (4) M. Nashaat and Yu. M. Shukrinov, Physics of Particles and Nuclei Letters, 17, 79. (2020).
  * (5) I. V. Bobkova , A. M. Bobkov, I. R. Rahmonov, A. A. Mazanik , K. Sengupta, and Yu. M. Shukrinov, Phys. Rev. B 102, 134505 (2020).
  * (6) Yu. M. Shukrinov, I. R. Rahmonov, K. Sengupta and A. Buzdin, Applied Physics Letters, 110, 182407, (2017).
  * (7) A. Buzdin, Phys. Rev. Lett. 101, 107005 (2008).
  * (8) F. Konschelle, A. Buzdin, Phys. Rev. Lett. 102, 017001 (2009).
  * (9) Yu. M. Shukrinov, I. R. Rahmonov, and K. Sengupta, Phys. Rev. B 99, 224513 (2019).
  * (10) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Phys. Rev. Lett. 113, 156401 (2014).
  * (11) Alexander A. Serga, Vasil S. Tiberkevich, Christian W. Sandweg, Vitaliy I. Vasyuchka, Dmytro A. Bozhko, Andril V. Chumak, Timo Neumann, Bjorn Obry, Gennadii A. Melkov, Andrei N. Slavin, and Burkard Hillebrands, Nat. Commun. 5, 3452 (2014).
  * (12) V. E. Demidov, O. Dzyapko, S. O. Demokritov, G. A. Melkov, and A. N. Slavin, Phys. Rev. Lett. 100, 047205 (2008).
  * (13) Shehzaad Kaka, Matthew R. Pufall, William H. Rippard, Thomas J. Silva, Stephen E. Russek, and Jordan A. Katine, Nature 437, 389 (2005).
  * (14) Laichuan Shen, Jing Xia, Xichao Zhang, Motohiko Ezawa, Oleg A. Tretiakov, Xiaoxi Liu, Guoping Zhao, and Yan Zhou, Phys. Rev. Lett. 124, 037202 (2020).
  * (15) A. I. Nikitchenko and N. A. Pertsev, Phys. Rev. Appl. 14, 034022 (2020).
  * (16) Mehrdad Elyasi, Yaroslav M. Blanter, and Gerrit E. W. Bauer, Phys. Rev. B 101, 054402 (2020).
  * (17) K. W. Moon, B. S. Chun,W. Kim, Z. Q. Qiu, and C. Hwang, Sci. Rep. 4, 6170 (2015).
  * (18) G. Wang, D. Chen, J. Lin, and X. Chen, IEEE Transactions on Industrial Electronics 46, 440 (1999).
  * (19) R. Almog, S. Zaitsev, O. Shtempluck, E. Buks, Phys. Rev. Lett. 98, 078103 (2007).
  * (20) A. H. Nayfeh and D. T. Mook, Nonlinear Oscillations (Wiley, New York, 1979).
  * (21) W. Chen, G. de Loubens, J.-M. L. Beaujour, J. Z. Sun, and A. D. Kent, Appl. Phys. Lett. 95, 172513 (2009).
  * (22) M. A. Silaev, Phys. Rev. B 102, 144521 (2020).
  * (23) Y. Yao, Q. Song, Y. Takamura, J. P. Cascales, W. Yuan, Y. Ma, Y. Yun, X. C. Xie, J. S. Moodera, and W. Han, Phys. Rev. B 97, 224414 (2018).
  * (24) Yuelei Zhao, Qi Song, See-Hun Yang, Tang Su1, Wei Yuan, Stuart S. P. Parkin, Jing Shi and Wei Han, Sci. Rep. 6, 22890 (2016).
  * (25) Shukrinov, Y.M., Rahmonov, I.R., Phys. Part. Nuclei 51, 816 (2020).
  * (26) Yu. M. Shukrinov, I. R. Rahmonov, and A. E. Botha., Low Temp. Phys. 46, 932 (2020).
  * (27) Ivana Kovacic, Michael J Brennan. The Duffing Equation : Nonlinear Oscillators and their Behaviour. — John Wiley and Sons, 2011.
  * (28) See Supplemental Material for details of our procedure to get the generalized Duffing equation from Landau-Lifshitz-Gilbert system of equations.
  * (29) M. A. W. Schoen, D. Thonig, M. L. Schneider, T. J. Silva, H. T. Nembach, O. Eriksson, O. Karis and J. M. Shaw, Nature Physics 12, 842 (2016).
  * (30) O. A. Kelly, A. Anane, R. Bernard, J. B. Youssef, C. Hahn, A. H. Molpeceres, C. Carrétéro, E. Jacquet, C. Deranlot, P. Bortolotti, R. Lebourgeois, J.-C. Mage, G. de Loubens, O. Klein, V. Cros, and A. Fert, Appl. Phys. Lett. 103, 082408 (2013).
  * (31) M. C. Onbasli, A. Kehlberger, D. H. Kim, G. Jakob, M. Kläui, A. V. Chumak, B. Hillebrands, and C. A. Ross, APL Mater. 2, 106102 (2014).
  * (32) A. Hrabec, F. J. T. Goncalves, C. S. Spencer, E. Arenholz, A. T. N’Diaye, R. L. Stamps, and C. H. Marrows, Phys. Rev. B 93, 014432 (2016).
  * (33) A. Yu. Rusanov, M. Hesselberth, J. Aarts, and A. I. Buzdin, _Phys. Rev. Lett._ 93, 057002 (2004).
  * (34) J. W. A. Robinson, F. Chiodi, M. Egilmez, G. B. Halász and M. G. Blamire, _Scientific Report_ 2, 699 (2012).