# Constant Curvature 3-branes in 5-D f(R) Bulk

Shafaq Gulzar Elahi<EMAIL_ADDRESS>Northwestern
University, Evanston, Illinois-60201, USA  Soumya Samrat Mandal
<EMAIL_ADDRESS>Purdue University, West Lafayette, Indiana-47907, USA
Soumitra SenGupta<EMAIL_ADDRESS>Indian Association for the Cultivation
of Science (IACS), Kolkata-700032, India

## Abstract

Braneworld models remain the most promising candidates to address several
important questions in low-energy particle phenomenology and cosmology. The
role of the moduli field(s) and its stabilization is an integral part of this
question. In this work, we show that a 5-dimensional warped braneworld model
with higher curvature gravity in bulk admits de-Sitter and anti de-Sitter
solutions on the branes. The remarkable feature of having a positive vacuum
energy on the visible brane is the presence of a metastable minimum and a
global minimum for the modulus potential. While the metastable minimum leads
to a consistent cosmological model of a bouncing universe, the concomitant
existence of the global minimum provides a vacuum for the modulus to roll down
to stability. Further, this model is shown to be consistent with the swampland
conjecture to qualify as a viable candidate in the low energy description of
string landscape

##  I Introduction

Ever since the proposal of Kaluza and Klein, there have been extensive studies
to explore the possible roles of the extra spatial dimension in our
3-dimensional universe. Various implications of the higher dimensional
geometry on lower dimensional physics range from gauge hierarchy problem [1],
the origin of dark matter/ dark energy [2]-[6], neutrino mass
hierarchy[7],[8], gauge/gravity duality[9], big bang or bouncing
cosmology[10]-[17], etc. One of the generic features of all kinds of extra-
dimensional models is the presence of one or more moduli fields associated
with the geometry of the internal manifold. While extracting an acceptable
theory for our 3+ 1 dimensional universe continues to be a big challenge, it
is also very hard to propose a mechanism that points towards a unique solution
to the geometry of the extra dimensions. While string theory helped to screen
out to some extent the nature of the internal manifold using the swampland
conjecture[18]-[21], however, the observational signatures of the dynamics of
the moduli field in cosmological/astrophysical scenarios are expected to play
a crucial role in exploring the possible signatures of extra dimensions. In
this context, the braneworld models are the most important candidates to bear
signatures of the extra dimensions and have been explored extensively. The
simple models for braneworld physics assume the existence of one extra
dimension and look for its signature on the lower dimensional brane. Although
a proper resolution of the origin of hierarchical scales in our brane namely
the electroweak and Planck scale continued to elude us for a long time, the
presence of extra dimension is found to offer a possible resolution to this
problem. However, it was realized that such models often bring in new
hierarchies in disguise through the scale associated with internal geometry. A
typical example of this was encountered in large extra-dimensional
models[11],[22]-[24]. As an alternative, a warped geometry model such as the
Randall-Sundrum model (RS) [25] could address different issues without
introducing any hierarchical scale in the model. The model is based on a
space-time which is fundamentally 5-dimensional with the background being
$AdS_{5}$ and the extra spatial dimension $y=\phi r_{c}$ (where $r_{c}$ is the
compactification radius) being angular and subject to $S_{1}/Z_{2}$
compactification. Consequently, there are two 3-branes at the orbifold fixed
points, the TeV/Visible and the Planck/Hidden Brane. The metric under
consideration is a non-factorizable one i.e. the 4-D part of the metric is
multiplied by an exponential function of the extra-dimension.The five-
dimensional action is:

$\mathcal{S}=\int{d^{5}x\,\sqrt{-g}(\text{M}^{3}\text{R}-\Lambda_{5})}-\int{d^{4}x\,\sqrt{-g_{i}}\mathcal{V}_{i}}$
(1)

where M is the fundamental 5-D Planck scale (we are working in natural units),
R is the 5-D Ricci Scalar,$\Lambda_{5}$ is the bulk cosmological constant,
$\cal{V}$i is the tension of the ith brane(i=hid(vis)) and $\eta_{\mu\nu}$is
the 4D metric.
The RS solution implies a negative bulk cosmological constant
$\Lambda_{5}=-12\text{M}^{3}k^{2}$ where $k=\sqrt{-\Lambda_{5}/12M^{3}}$ and
brane tensions
$\cal{V}$${}_{\text{hid}}=-$$\cal{V}$${}_{\text{vis}}=12\text{M}^{3}k^{2}$.
The Solution is:

$ds^{2}=e^{-2kr_{c}|\phi|}\eta_{\mu\nu}dx^{\mu}dx^{\nu}+r_{c}^{2}d\phi^{2}$
(2)

A natural extension to the RS model is to include the higher curvature terms
in the action which is motivated by the fact that the 5-D is endowed with a
large cosmological constant $\mathcal{O}$$\sim\rm M$ which in turn ensures
that no intermediate scale unlike large extra-dimensional models develop in
the theory. However, the gravity action with such a large background curvature
most naturally should be described higher curvature correction to Einstein’s
action. In [27], we considered the simplest extension to the Einstein-Hilbert
Action in the bulk in the form of $\rm f(R)=R+\alpha R^{2}-\Lambda_{5}$ since
the bulk has the scalar curvature $\rm R\sim\mathcal{O}(M^{2})$ and found the
corresponding modification to the warp factor in a perturbative manner. We
also showed that the presence of additional gravitational degrees of freedom
in the bulk due to f(R) gravity could automatically generate a stabilizing
potential for the modulus/radion field. Meanwhile, in an earlier work [28],
the RS model was generalized to include non-flat maximally symmetric 3-branes.
This work produced a wide class of warped solutions where the RS model was
just a special case when the brane cosmological constant was zero. The de-
Sitter generalization was particularly relevant to address various
cosmological scenarios with the accelerated expansion of the universe. One
interesting outcome of considering a non-flat brane was subsequently shown
[29] in the context of modulus stabilization. The outcome indicated the
generation of a nonzero modulus potential where the ground state of the
modulus indicated a metastable vacuum. The potential and the vacuum value of
the modulus were proportional to the brane cosmological constant which
disappears if the brane cosmological constant is set to zero as in the
original RS case. It was further shown in [17] that the existence of such a
modulus potential with a metastable ground state can lead to a transient
phantom phase resulting in a bouncing cosmological model with a low scalar
tensor ratio consistent with Planck data. However, the possibility of having
an appropriate global minimum for the radion field could not be resolved in
this model. Moreover, despite several interesting outcomes, this non-flat
generalization was done with Einstein action in the bulk though the bulk
continued to have a vacuum curvature of the order of the Planck scale. The
presence of higher curvature terms in the action should be a more appropriate
description in such a scenario. Motivated by our earlier work [27]
generalizing the warped geometry model with higher curvature gravity with flat
3-branes, in this work, we want to analyze the effect of considering a non-
zero 4-D cosmological constant in our universe i.e. the 4-D metric is no
longer $\eta_{\mu\nu}$, instead, it has a constant curvature with a positive
or negative brane cosmological constant leading to a de-Sitter or anti de-
Sitter 3-brane. Our goal is to explore the nature of the modulus potential and
the consequent modulus stabilization in such a scenario.

##  II The Model

Throughout the paper, we will adapt the mostly positive convention for the
metric i.e. $\eta_{\mu\nu}=(-,+,+,+)$, the Greek indices $\alpha,\beta$..
denote 4-D indices whereas the Latin indices a,b, A, B… denote the 5-D
indices.
Consider the metric ansatz:

$ds^{2}=e^{-2\text{A(y)}}\overline{g}_{\mu\nu}(x)dx^{\mu}dx^{\nu}+\text{B(y)}^{2}dy^{2}$
(3)

Here $\overline{g}_{\mu\nu}$ is the 4-D metric.The non-zero Christoffel
symbols for this metric are :
$\Gamma^{5}_{55}=\text{B(y)}^{\prime}/\text{B(y)}$,
$\Gamma^{5}_{\mu\nu}=\overline{g}_{\mu\nu}e^{-2\text{A(y)}}\text{A(y)}^{\prime}/\text{B(y)}^{2}$
, $\Gamma^{\mu}_{5\nu}=-\delta^{\mu}_{\nu}\text{A(y)}^{\prime}$and
${}^{4}\Gamma^{\alpha}_{\mu\nu}$. Here, prime denotes the derivative w.r.t.
the extra-dimensional coordinate y. Also,

$\displaystyle\text{R}={}^{4}\mathcal{R}e^{2\text{A(y)}}+\frac{-4}{\text{B}(y)^{3}}\left[2\text{A(y)}^{\prime}\text{B}^{\prime}(y)+\text{B}(y)\left(5\text{A(y)}^{\prime
2}-2\text{A(y)}^{\prime\prime}\right)\right]$ (4)
$\displaystyle\text{R}_{55}=-4\left[\frac{\text{A}^{\prime}(y)\text{B}^{\prime}(y)}{\text{B}(y)}-\text{A}^{\prime\prime}(y)+\text{A}^{\prime}(y)^{2}\right]$
(5)
$\displaystyle\text{R}_{\mu\nu}={}^{4}\rm\mathcal{R}_{\mu\nu}+\frac{1}{\text{B}(y)^{3}}\overline{g}_{\mu\nu}e^{-2\text{A(y)}}\left[\text{B}(y)\left(\text{A}^{\prime\prime}(y)-4\text{A}^{\prime}(y)^{2}\right)-\text{A}^{\prime}(y)\text{B}^{\prime}(y)\right]$
(6)

where ${}^{4}\mathcal{R}$ is the 4-D scalar curvature. We will consider that
the 5-D scalar curvature is only a function of the extra-dimension i.e., $\rm
R=R(y)$. This is because we chose ${}^{4}\mathcal{R}=4\Omega$ and
${}^{4}\mathcal{R}_{\mu\nu}=\Omega\overline{g}_{\mu\nu}$, where $\Omega$ is a
constant and can be both positive or negative. We will set $\rm M=1$ for
simplification. We want to explore constant curvature 3-brane solutions of RS-
like space-time when the bulk action is given by:

$\mathcal{S}=\int{d^{5}x\,\sqrt{-{g}}\leavevmode\nobreak\
\text{f(R)}}-\int{d^{4}x\,\sqrt{-g_{i}}\mathcal{V}_{i}}$ (7)

where $\rm f(R)=R+\alpha R^{2}-\Lambda_{5}$. The field equation for a general
f(R) action is:

$\text{R}_{\text{MN}}\text{f}_{\text{R}}(\text{R})-\frac{1}{2}g_{\text{MN}}\text{f(R)}+g_{\text{MN}}\square\text{f}_{\text{R}}(\text{R})-\nabla_{\text{M}}\nabla_{\text{N}}\text{f}_{\text{R}}(\text{R})=\frac{1}{2}\textbf{T}_{\text{MN}}$
(8)

where $\textbf{T}_{\text{MN}}$ is the energy-momentum tensor corresponding to
matter of some kind. Substituting 4, 5 and 6 in 8, we obtain the following
equations (for the bulk):

$\displaystyle\rm-\frac{1}{B(y)^{4}}\leavevmode\nobreak\
2\left(-3B(y)^{4}A^{\prime}(y)^{2}+\Omega
e^{2A(y)}B(y)^{6}+3k^{2}B(y)^{6}\right)$
$\displaystyle\rm-\frac{2\alpha}{B(y)^{4}}\leavevmode\nobreak\
\Biggl{\\{}80A^{\prime}(y)^{2}B^{\prime}(y)^{2}+8\Omega
e^{2A(y)}B(y)^{4}A^{\prime}(y)^{2}+32B(y)A^{\prime}(y)\left(2B^{\prime}(y)\left(2A^{\prime}(y)^{2}\right.\right.$
$\displaystyle\rm\left.\left.-A^{\prime\prime}(y)\right)-A^{\prime}(y)B^{\prime\prime}(y)\right)+4B(y)^{2}\left(-4A^{\prime\prime}(y)^{2}+5A^{\prime}(y)^{4}+8A^{\prime\prime\prime}(y)A^{\prime}(y)\right.$
$\displaystyle\rm\left.-32A^{\prime}(y)^{2}A^{\prime\prime}(y)\right)+4\Omega^{2}e^{4A(y)}B(y)^{6}\Biggr{\\}}=0$
(9)
$\displaystyle\rm-\frac{1}{B(y)^{7}}\left(3B(y)^{5}A^{\prime\prime}(y)-3B(y)^{4}A^{\prime}(y)B^{\prime}(y)-6B(y)^{5}A^{\prime}(y)^{2}+B(y)^{7}\left(\Omega
e^{2A(y)}+6k^{2}\right)\right)$
$\displaystyle\rm-\frac{1}{B(y)^{7}}8\alpha\Biggl{\\{}\Omega
e^{2A(y)}B(y)^{4}A^{\prime}(y)B^{\prime}(y)+30A^{\prime}(y)B^{\prime}(y)^{3}+10B(y)B^{\prime}(y)\left(3B^{\prime}(y)\left(2A^{\prime}(y)^{2}\right.\right.$
$\displaystyle\rm\left.\left.-A^{\prime\prime}(y)\right)-2A^{\prime}(y)B^{\prime\prime}(y)\right)+\Omega
e^{2A(y)}B(y)^{5}\left(A^{\prime}(y)^{2}-A^{\prime\prime}(y)\right)+B(y)^{2}\left(12A^{\prime\prime\prime}(y)B^{\prime}(y)+8A^{\prime\prime}(y)B^{\prime\prime}(y)\right.$
$\displaystyle\rm\left.-16A^{\prime}(y)^{2}B^{\prime\prime}(y)+37A^{\prime}(y)^{3}B^{\prime}(y)+2A^{\prime}(y)\left(B^{\prime\prime\prime}(y)-36A^{\prime\prime}(y)B^{\prime}(y)\right)\right)+B(y)^{3}\left(-2A^{\prime\prime\prime\prime}(y)+12A^{\prime\prime}(y)^{2}\right.$
$\displaystyle\rm\left.\left.+5A^{\prime}(y)^{4}+16A^{\prime\prime\prime}(y)A^{\prime}(y)-37A^{\prime}(y)^{2}A^{\prime\prime}(y)\right)\right)\Biggr{\\}}=0$
(10)

###  II.1 Perturbative Solution

Consider $\alpha$ to be small such that $\text{R}^{2}$ in f(R) is a small
correction over R solution, In this backdrop, consider the following ansatz
for A(y) and B(y):

$\displaystyle\rm e^{-A(y)}=\omega\leavevmode\nobreak\
sinh\left[ln\left(\frac{c_{2}}{\omega}\right)-k\leavevmode\nobreak\
y\right]+\alpha A_{1}(y)$ (11) $\displaystyle\rm B(y)=1+\alpha
b_{0}\leavevmode\nobreak\ y$ (12)

where $\omega^{2}=\Omega/3k^{2}$ for $\Omega>0$ and
$\rm\omega\leavevmode\nobreak\ \rm
sinh\left[ln\left(\frac{c_{2}}{\omega}\right)-k\leavevmode\nobreak\ y\right]$
represents the unperturbed solution as per [28] i.e. by setting $\alpha\to 0$
in IIand II. The constant $c_{2}=\rm\sqrt{1+\sqrt{1+\omega^{2}}}$ is obtained
by setting the boundary condition at the hidden brane $\rm e^{-A(y=0)}=1$.We
have chosen to fix $\rm B(y)$ as per [27] and find a solution for $\rm A(y)$.
$\cal{O}$$(\alpha^{0})$ is trivially satisfied and at $\cal{O}$$(\alpha)$, II
and II reduce to:

$\displaystyle\rm-\frac{1}{c_{2}\leavevmode\nobreak\ k}\leavevmode\nobreak\
3\Omega
e^{ky}\left(kA_{1}(y)-A_{1}(y)^{\prime}\right)+9c_{2}\leavevmode\nobreak\
ke^{-ky}\left(A_{1}(y)^{\prime}+k\leavevmode\nobreak\ A_{1}(y)\right)$
$\displaystyle\rm+\frac{1}{6c_{2}^{2}\leavevmode\nobreak\
k^{2}}\Omega^{2}e^{2ky}\left(3b_{0}\leavevmode\nobreak\
y+10k^{2}\right)+\frac{3}{2}c_{2}^{2}\leavevmode\nobreak\
k^{2}\leavevmode\nobreak\ e^{-2ky}\left(3b_{0}\leavevmode\nobreak\
y+10k^{2}\right)+3b_{0}\leavevmode\nobreak\ y\Omega-10k^{2}\Omega=0$ (13)
$\displaystyle\rm 54c_{2}^{3}\leavevmode\nobreak\ k^{3}\leavevmode\nobreak\
e^{ky}\left(-A_{1}(y)^{\prime\prime}+2k\leavevmode\nobreak\
A_{1}(y)^{\prime}+3k^{2}A_{1}(y)\right)-18c_{2}\leavevmode\nobreak\ k\Omega
e^{3ky}\left(-A_{1}(y)^{\prime\prime}-2k\leavevmode\nobreak\
A_{1}(y)^{\prime}\right.$ $\displaystyle\rm\left.\rm+3\leavevmode\nobreak\
k^{2}A_{1}(y)\right)+3b_{0}\left(9c_{2}^{4}\leavevmode\nobreak\
k^{4}\leavevmode\nobreak\ (4ky-1)+\Omega^{2}e^{4ky}\leavevmode\nobreak\
(4ky+1)\right)+40k^{3}\left(\Omega e^{2ky}-3c_{2}^{2}\leavevmode\nobreak\
k^{2}\right)^{2}=0$ (14)

The solution is:

$\rm A_{1}(y)=\frac{1}{12}c_{2}\leavevmode\nobreak\
ke^{-ky}\left(-3b_{0}\leavevmode\nobreak\
y^{2}+20(c_{2}-2)k-20k^{2}\leavevmode\nobreak\ y\right)+\frac{\Omega
e^{ky}\left(20c_{2}\leavevmode\nobreak\ k-y\left(3b_{0}\leavevmode\nobreak\
y+20k^{2}\right)\right)}{36c_{2}\leavevmode\nobreak\ k}$ (15)

As $\Omega\to 0$,we recover the solution obtained in [27] and as both
$\alpha\leavevmode\nobreak\ \rm and\leavevmode\nobreak\ \Omega\to 0$, 15
reduces to the RS solution [25]. We can perform a similar analysis for
$\Omega<0$ considering the hyperbolic cosine function instead in the ansatz
solution.

###  II.2 Brane Tensions

The field equations II and II are supplemented by the boundary conditions (up
to leading order in $\alpha$)

$\left[A^{\prime}(y)\right]_{i}=\frac{\epsilon_{i}}{12}\mathcal{V}_{i}$ (16)

where $\epsilon_{\mathrm{hid}}=-\epsilon_{\mathrm{vis}}=1$. We want to
calculate the leading order corrections to the brane tensions in this setup.
Using 15 in 11, we have the leading order correction to the warp factor.
Substitute this expression in 16, we obtain the following expressions for
brane-tensions:

$\displaystyle\rm\mathcal{V}_{\text{hid}}=$
$\displaystyle\rm\frac{12k\left[1+\frac{\omega^{2}}{\text{c2}^{2}}+\frac{10}{3}\alpha
k^{2}(c_{2}-1)\left(1-\frac{\omega^{2}}{c_{2}^{2}}\right)\right]}{\left(1-\frac{\omega^{2}}{c_{2}^{2}}\right)+\frac{10}{3}\alpha
k^{2}\left[\frac{\omega^{2}}{c_{2}}+(c_{2}-2)\right]}$ (17)
$\displaystyle\rm\mathcal{V}_{vis}=$
$\displaystyle\rm-\frac{12k\Bigl{\\{}1+\frac{\omega^{2}}{c_{2}^{2}}e^{2\pi
kr_{c}}+\alpha\left[\mathcal{B}(r_{c})+\frac{\mathbf{G}(r_{c})\omega^{2}}{3c_{2}^{2}}e^{2\pi
kr_{c}}\right]\Bigr{\\}}}{1-\frac{\omega^{2}}{c_{2}^{2}}e^{2\pi
kr_{c}}+\alpha\left[\frac{\mathcal{H}(r_{c})\Omega}{18c_{2}^{2}k}e^{2\pi
kr_{c}}+\mathcal{J}(r_{c})\right]}$ (18)

where

$\displaystyle\mathcal{B}(r_{c})=\pi r_{c}\left(-\frac{1}{2}\pi
b_{0}kr_{c}+b_{0}-\frac{10k^{3}}{3}\right)+\frac{10}{3}(c_{2}-1)k^{2}$ (19)
$\displaystyle\mathbf{G}(r_{c})=\pi r_{c}\left(\frac{3}{2}\pi
b_{0}kr_{c}+3b_{0}+10\text{c2}k^{3}\right)-10(c_{2}-1)k^{2}$ (20)
$\displaystyle\mathcal{H}(r_{c})=20c_{2}k-\pi r_{c}\left(3\pi
b_{0}r_{c}+20k^{2}\right)$ (21)
$\displaystyle\mathcal{J}(r_{c})=\frac{1}{6}k\left(-3\pi^{2}b_{0}r_{c}^{2}+20(c_{2}-2)k-20\pi
k^{2}r_{c}\right)$ (22)

In the limit $\alpha\to 0$,

$\displaystyle\rm\mathcal{V}_{\text{hid}}=12k\left[\frac{1+\frac{\omega^{2}}{c_{2}^{2}}}{1-\frac{\omega^{2}}{c_{2}^{2}}}\right];\hskip
14.22636pt\mathcal{V}_{\text{vis}}=-12k\left[\frac{1+\frac{\omega^{2}}{c_{2}^{2}}e^{2k\pi
r_{c}}}{1-\frac{\omega^{2}}{c_{2}^{2}}e^{2k\pi r_{c}}}\right]$ (23)

consistent with the results of [28]. In this case, $\cal{V}$hid is always
positive and $\cal{V}$vis is negative for the entire range of $\Omega>0$ for
$\alpha$ and $b_{0}$ positive.

##  III Modulus Stabilization

As has been discussed earlier modulus stabilization is an essential criterion
for warped braneworld scenarios. Without a stable value for the moduli, no
low-energy phenomenology can be extracted from the higher dimensional models.
As a result, the study of the stabilization mechanism of higher dimensional
models is of utmost importance to make the model viable and relevant. The
stabilized modulus field has wide phenomenological implications in the domain
of cosmology and collider physics. Goldberger and Wise [30]-[31] came up with
a stabilization mechanism for the original Randall Sundrum model, which
resolved the gauge hierarchy problem. In their work, a bulk scalar field was
introduced which in turn generated a potential to stabilize the modulus field
to the desired vev. However, they ignored the back reaction of the bulk scalar
field onto the background geometry, which was later taken into account in
[32]. Further, keeping in mind the scalar-tensor equivalence in $f(R)$ gravity
theories [33]-[34], efforts were made to achieve modulus stabilization with
the conformally transformed metric[35]-[38]. In our previous work, we could
achieve modulus stabilization in an original RS-like scenario in a $f(R)$
gravity background, which included higher curvature correction to the Einstein
Hilbert action. On the other hand, it was shown that upon generalizing the
original RS model to a non-flat braneworld scenario, with Einstein Hilbert
action in the bulk, the modulus potential had a metastable minimum [29]. In
our present work, we are interested in examining the nature of the modulus
potential in a generalized non-flat braneworld scenario in $f(R)$ gravity
background. The vacuum expectation value of the modulus field, $T$, will yield
the compactification radius, $r_{c}$, in our model.

Considering only the bulk action from 1, after integrating out the extra
dimension to obtain an effective 4-D theory, we will get a kinetic term for
the modulus field, it’s non-minimal coupling to the scalar curvature and its
potential. The analytic expression for V(T) in terms of basic functions is not
possible. Therefore, we have obtained a numerical solution for the potential
for different values of the model parameters depicted by the plots below.

$\rm V(T)=\rm 2\int_{0}^{\pi}{d\phi\,T\,e^{-4\text{A({\hbox{\phi\,
,T}})}}\,\text{B({\hbox{\phi\,
,T}})}\,(\text{R}+\alpha\text{R}^{2}-\Lambda_{5})}\Bigg{|}_{\rm
potential\leavevmode\nobreak\ part}$ (24)

(a) V(T) vs T for $\alpha=0.2$ and $b_{0}=10$.

(b) V(T) vs T for $\alpha=0.2$ and $b_{0}=5$

(c) V(T) vs T for $\alpha=0.16$ and $b_{0}=10$

(d) V(T) vs T for $\alpha=0.22$ and $b_{0}=10$

Figure 1: Plots of modulus potential for different values of $\alpha$ and
$b_{0}$.

Stabilization can be achieved for $0<\alpha<1$, suggesting that our
perturbative approach works. This stabilization criterion is the same as the
one that we had obtained for the flat brane scenario [27]. As
$\alpha\xrightarrow{}0$, only a metastable minimum occurs. This in turn
further reinstates the fact that stable minima mainly arose due to the
$\alpha$ correction to the warp factor, and as this correction term tends
towards zero, the result is the same as the result shown in [29].
To be consistent with $f^{\prime}(R)>0$, we obtain the following bound on the
model parameters : $\alpha>0$ and $b_{0}>0$. Hence, our model can be applied
to the cosmology sector as $\alpha>0$. Therefore, we have $0<\alpha<1$ and
$b_{0}>0$ to obtain modulus stabilisation in ghost-free quadratic f(R)
gravity.
We can further see that for a fixed value of $\alpha$, as the value of $b_{0}$
increases, the minima of the potential shifts and now the stabilization occurs
for a smaller value of $T(x)$. In 1(a), stabilization occurs at $T(x)\approx
0.4$( for $b_{0}=10$ ), whereas in 1(b), stabilization occurs at $T(x)\approx
0.55$( for $b_{0}=5$ ), for the same value of $\alpha$ ( $\alpha$ = 0.2 ).For
a fixed value of $b_{0}$, as the value of $\alpha$ increases, stabilization
occurs for a smaller value of $T(x)$. For instance, in 1(c), the stabilization
occurs at $T(x)\approx 0.5$ ($\alpha=0.16$), while in 1(d), the stabilization
occurs at $T(x)\approx 0.35$ ($\alpha=0.22$), for the same value of $b_{0}$ (
$b_{0}=10$ ).

###  III.1 Swampland Conjecture

As string theory predicts the existence of extra dimensions for its
consistency[20], the connection between any effective higher dimensional model
on low energy brane and a class of string vacuum is always worthwhile to
explore[18, 19]. It is important to explore whether the proposed higher
dimensional model satisfies the swampland conjecture and thereby can provide a
consistent low-energy description of the vast string landscape. In this
context, we focus our attention on the restrictions that the modulus potential
needs to satisfy to obey the swampland conjecture. In our model,
$\rm|\partial_{T}V(T)/V(T)|\sim$$\cal{O}$$(1)$ near the vacuum expectation
value of the modulus field $T$. For $\alpha=0.2$ and $b_{0}=10$,
$\rm|\partial_{T}V(T)/V(T)|\sim$$\cal{O}$$(1)$ for $0.35<T<0.45$. Whereas, for
$\alpha=0.2$ and $b_{0}=5$, $\rm|\partial_{T}V(T)/V(T)|\sim$$\cal{O}$$(1)$ for
$0.4<T<0.7$. It can be shown that for a certain choice of model parameters,
$V(T)>0$. Thus the conditions imply that the model qualifies to be a
consistent low-energy description of string landscape.

##  IV Conclusion

Branes with a non-zero cosmological constant are natural generalizations of
the warped geometry model to explore the cosmological implications of
braneworld models. Such a scenario was earlier shown to originate when bulk
curvature contribution fails to cancel exactly the brane tension due to the
presence of additional matter. The scenario with a de-Sitter brane was also
shown to exhibit the presence of a metastable minimum for the modulus field
which was absent for a pure RS-like scenario with flat branes. It was also
exhibited that such a metastable minimum plays a crucial role in triggering a
bouncing cosmological scenario leading to possible answers to avoid
cosmological singularity. The question however remained what happens to the
moduli field once it crosses the metastable minimum and rolls down? In this
work, we have demonstrated that the presence of a higher curvature term in the
bulk leads to a stable minimum alongside a metastable minimum. This eventually
leads to an auto-modulus stabilization by the geometry itself and brings in
interesting phenomenological consequences of a radion field with non-zero vev.
The evolution of the scale factor in the presence of such a modulus potential
with a non-minimal kinetic term for the modulus field may now lead to a
bouncing scenario for an appropriate choice of the parameters. We have also
verified that the swampland criteria are satisfied with the proper choice of
model parameters and therefore the model proposes a consistent effective
description of string landscape with viable ultraviolet description in the
quantum regime.

## Acknowledgement

SGE and SSM would like to extend their gratitude to IACS where they carried
out this work.

## References

  * [1] H. Hatanaka, T. Inami and C. S. Lim,“The Gauge hierarchy problem and higher dimensional gauge theories,” Mod. Phys. Lett. A 13, 2601-2612 (1998) [arXiv:hep-th/9805067 [hep-th]].
  * [2] J. R. Primack,“Dark matter and structure formation,”[arXiv:astro-ph/9707285 [astro-ph]].
  * [3] Frenk, Carlos S., and Simon DM White. "Dark matter and cosmic structure." Annalen der Physik 524.9-10 (2012): 507-534.
  * [4] Joyce, Austin, et al. "Beyond the cosmological standard model." Physics Reports 568 (2015): 1-98.
  * [5] Nojiri, Shin’ichi, and Sergei D. Odintsov. "Dark energy, inflation and dark matter from modified F (R) gravity." arXiv preprint arXiv:0807.0685 (2008).
  * [6] Sudhakantha Girmohanta and Robert Shrock."Extra-dimensional model of dark matter". Phys. Rev. D 104, 115021.
  * [7] Angel, Paul W., Nicholas L. Rodd, and Raymond R. Volkas. "Origin of neutrino masses at the LHC:$\Delta$ L= 2 effective operators and their ultraviolet completions." Physical Review D 87.7 (2013): 073007.
  * [8] Huang, Wei-Chih, and Frank F. Deppisch. "Dark matter origins of neutrino masses." Physical Review D 91.9 (2015): 093011.
  * [9] Klebanov, Igor R., and Matthew J. Strassler. "Supergravity and a confining gauge theory: Duality cascades and $\chi$SB-resolution of naked singularities." Journal of High Energy Physics 2000.08 (2000): 052.
  * [10] Chialva, Diego, Gary Shiu, and Bret Underwood. "Warped reheating in multi-throat brane inflation." Journal of High Energy Physics 2006.01 (2006): 014.
  * [11] Starkman, Glenn D., Dejan Stojkovic, and Mark Trodden. "Large extra dimensions and cosmological problems." Physical Review D 63.10 (2001): 103511.
  * [12] Mazumdar, Anupam. "Extra dimensions and inflation." Physics Letters B 469.1-4 (1999): 55-60.
  * [13] Senatore, Leonardo. "Lectures on inflation." Theoretical Advanced Study Institute in Elementary Particle Physics: new frontiers in fields and strings (2016): 447-543.
  * [14] Bamba, Kazuharu, et al. "Bounce cosmology from F (R) gravity and F (R) bigravity." Journal of Cosmology and Astroparticle Physics 2014.01 (2014): 008.
  * [15] Das, Ashmita, et al. "Bouncing cosmology from warped extra-dimensional scenario." The European Physical Journal C 77.12 (2017): 1-9.
  * [16] Ijjas, Anna, and Paul J. Steinhardt. "Bouncing Cosmology made simple." Classical and Quantum Gravity 35.13 (2018): 135004.
  * [17] I. Banerjee, T. Paul and S. SenGupta,“Bouncing cosmology in a curved braneworld,”JCAP 02, 041 (2021)
  * [18] H. Ooguri and C. Vafa,“On the Geometry of the String Landscape and the Swampland,” Nucl. Phys. B 766, 21-33 (2007)
  * [19] T. D. Brennan, F. Carta and C. Vafa,“The String Landscape, the Swampland, and the Missing Corner,”PoS TASI2017, 015 (2017) [arXiv:1711.00864 [hep-th]].
  * [20] E. Palti,“The Swampland: Introduction and Review,”Fortsch. Phys. 67, no.6, 1900037 (2019) [arXiv:1903.06239 [hep-th]].
  * [21] G. Montefalcone, P. J. Steinhardt and D. H. Wesley,“Dark energy, extra dimensions, and the Swampland,”JHEP 06, 091 (2020)
  * [22] V. A. Rubakov,“Large and infinite extra dimensions: An Introduction,”Phys. Usp. 44, 871-893 (2001) [arXiv:hep-ph/0104152 [hep-ph]].
  * [23] Dvali, Gia, and Alexei Yu Smirnov. "Probing large extra dimensions with neutrinos." Nuclear Physics B 563.1-2 (1999): 63-81.
  * [24] N. Arkani-Hamed, S. Dimopoulos and G. R. Dvali,“The Hierarchy problem and new dimensions at a millimeter,” Phys. Lett. B 429 (1998), 263-272
  * [25] L. Randall and R. Sundrum, “A Large mass hierarchy from a small extra dimension,” Phys. Rev. Lett. 83 (1999), 3370-3373
  * [26] Riley, Bernard. "Standard Model scales from warped extra dimensions." arXiv:0809.0111 (2008).
  * [27] S. G. Elahi, S. S. Mandal and S. SenGupta,“A novel modulus stabilization mechanism in higher dimensional f(R) Gravity,” [arXiv:2212.03276 [gr-qc]].
  * [28] S. Das, D. Maity and S. SenGupta, “Cosmological constant, brane tension and large hierarchy in a generalized Randall-Sundrum braneworld scenario,” JHEP 05, 042 (2008)
  * [29] I. Banerjee and S. SenGupta,“Modulus stabilization in a non-flat warped braneworld scenario,” Eur. Phys. J. C 77, no.5, 277 (2017)
  * [30] W. D. Goldberger and M. B. Wise,“Modulus stabilization with bulk fields,” Phys. Rev. Lett. 83, 4922-4925 (1999)
  * [31] W. D. Goldberger and M. B. Wise,“Phenomenology of a stabilized modulus,”Phys. Lett. B 475 (2000), 275-279
  * [32] Csaki, Csaba, Michael L. Graesser, and Graham D. Kribs. "Radion dynamics and electroweak physics." Physical Review D 63.6 (2001): 065002.
  * [33] De Felice, A., Tsujikawa, S. f(R) Theories. Living Rev. Relativ. 13, 3 (2010). https://doi.org/10.12942/lrr-2010-3
  * [34] S. Nojiri, S. D. Odintsov and V. K. Oikonomou,“Modified Gravity Theories on a Nutshell: Inflation, Bounce and Late-time Evolution,”Phys. Rept. 692 (2017), 1-104 [arXiv:1705.11098 [gr-qc]].
  * [35] Banerjee, Indrani, Tanmoy Paul, and Soumitra SenGupta. "Critical analysis of modulus stabilization in a higher dimensional F (R) gravity." Physical Review D 104.10 (2021): 104018.
  * [36] Paul, Tanmoy, and Soumitra SenGupta. "Radion tunneling in modified theories of gravity." The European Physical Journal C 78.4 (2018): 1-9.
  * [37] Das, Ashmita, et al. "Radion stabilization in higher curvature warped spacetime." The European Physical Journal C 78.2 (2018): 1-11.
  * [38] Chakraborty, Sumanta and Soumitra Sengupta. “Solving higher curvature gravity theories.” The European Physical Journal C 76 (2016): 1-13.


# Addition to the Gnevyshev-Ohl rule and prediction of solar cycle 25

Yu. A. Nagovitsyn1,2 1Central Astronomical Observatory at Pulkovo, Saint-
Petersburg, Russia
2State University of Aerospace Instrumentation, St. Petersburg, Russia V. G.
Ivanov1 1Central Astronomical Observatory at Pulkovo, Saint-Petersburg, Russia
2State University of Aerospace Instrumentation, St. Petersburg, Russia

###### Abstract

In addition to the Gnevyshev-Ohl rule (GOR), the relation of the odd cycle
with the subsequent even one in the 22-year Hale solar cycle was found. It is
shown that 3 years before the 11-year minimum $m$, the value of the relative
sunspot number $S\\!N$ in an odd cycle is closely related to the value of the
maximum in the next even cycle (correlation coefficient $\rho=0.94$), and the
same relation of an odd cycle with the previous even one is weaker. Like GOR,
cycles are linked in pairs, but opposite to the Rule.

Based on this result, we propose to use $S\\!N_{m-3}$ on the descending phase
of the previous odd cycle as a precursor of the subsequent EVEN cycle (Figure
3a) — a precursor called MI3E. For the prediction of an odd cycle or a
prediction without consideration of parity (as in the article by brajsa22),
this method gives less reliable results.

To predict the amplitude of an ODD cycle, we propose to use the precursor of
the seventh year to its maximum $M$ MA7O — $S\\!N_{M-7}$ on the descending
phase of the previous even cycle (Figure 3b). It turned out that in this case,
we can also predict the years near the maximum with a high correlation
coefficient ($\rho=0.90{-}0.94$).

Thus, the proposed approaches allow us to predict cycles of different parity.
According to our prediction, the current solar Cycle 25 in 2023 will reach a
maximum of 154 units with a prediction error of $\pm 25$ (68% confidence) and
$\pm 53$ (95% confidence). In 2024, $S\\!N$ will be almost as high as in 2023
— 147 units, so with smaller time averaging scales, the maximum will fall at
the end of 2023.

###### keywords:

Solar activity; Sunspots; Solar cycle prediction; Gnevyshev-Ohl rule

##  1 Introduction

As noted in (nag09), “The empirical Gnevyshev-Ohl rule (1948), below referred
to as GOR or the Rule, is one of the most puzzling properties of the solar
cyclicity.” gnevohl48, based on the consideration of a relative sunspot
numbers $S\\!N$ from 1700 to 1944, found that the area under the curve (total
power) of the 11-year cycle $\sum S\\!N$ has a high correlation for the pair
of even – subsequent odd cycles and weakly correlates for the pair of odd –
subsequent even cycles. Thus, they concluded that the 22-year Hale cycle
begins with an even cycle. This is a strange circumstance since the resulting
rule divides cycles into pairs, and neighboring pairs are weakly related. One
of the questions we raise in this article is: in addition to the parameter
$\sum S\\!N$, are there cycle parameters that connect consecutive cycles in
different ways, including in the odd – subsequent even cycle pair?

In addition, we touch on a purely practical question of solar cycle
prediction: how will look the recently started Cycle 25, and what will be its
amplitude when the moment of maximum comes. As it turned out, in the context
under consideration, the prediction depends, and quite critically, on the
parity of cycles, i.e., the global organization of the Sun’s magnetic field in
the light of the Hale cycle.

##  2 Gnevyshev-Ohl rule

Almost three-quarters of a century has passed since the publishing of the
article by Gnevyshev and Ohl. Since that time, new 11-year cycles appeared,
adding new statistical data, and the version of $S\\!N$ used by these authors
was replaced by a version called $S\\!N$ 2.0 (clette14, 2016). Therefore,
using SILSO data https://www.sidc.be/silso/ let’s plot the dependencies $\sum
S\\!N_{E}=f_{1}(\sum S\\!N_{O})$, $\sum S\\!N_{O}=f_{2}(\sum S\\!N_{E})$, and
calculate the correlation coefficients $\rho$ for them according to modern
material, starting from the year 1700. The “E” and “O” subscripts hereinafter
indicate even and odd cycles, respectively. As in the original work of
gnevohl48, we exclude from consideration a pair of Cycles 4 and 5. Figure 1
illustrates the results. The confidence interval $\pm 1\sigma$ is selected
(the probability of falling within is 68%).

Figure 1:  Empirical relations (squares) and regressions (solid lines) for two
consecutive cycles: (a) $\sum S\\!N_{O}=(270\pm 99)+(0.84\pm 0.12)\cdot\sum
S\\!N_{E}$, $\rho=0.906$; (b) $\sum S\\!N_{E}=(480\pm 220)+(0.38\pm
0.24)\cdot\sum S\\!N_{O}$, $\rho=0.438$, where $\rho$ is the Pearson
correlation coefficient. Standard deviations are dashed lines, and the
prediction intervals (with 68% confidence) are dotted lines. The circled
point, which was not used in the regressions’ calculation, corresponds to the
prediction of the Cycle 25 (see Figure 6b).

There is an evident confirmation of the GOR. Thus, this rule captures the
causation of the odd cycle with the previous even one for the modern version
of $S\\!N$ 2.0 as well (Figure 1a), while the even cycle within the Rule is
not related to the previous odd one (Figure 1b).

A summary of possible methods for predicting the amplitudes of 11-year cycles
can be found in the review by petrovay20. In particular, it contains the
information about precursor methods.

The precursor methods are based on the idea that the cycle begins before the
sunspot solar cycle minimum between the old and the new cycle. This idea does
not contradict the general approach of the solar dynamo and can be used to
predict the amplitude of the future maximum in the epoch near the previous
minimum. A.I. Ohl (1966) was one of the first to propose such a method: as a
precursor, he used geomagnetic activity. The geomagnetic activity has often
been used for solar activity cycle prediction (for example, pesnell14). sval05
used a polar field. wilson98 proposed a comprehensive prediction for a number
of cycle parameters. “Older” and less-known prediction methods, including
those with the approach of precursors, are contained in Vitinsky’s monograph
(1973).

According to macint20, the international group of experts Solar Cycle 25
Prediction Panel, (SC25PP) concluded that the sunspot solar Cycle 25 will be
similar in amplitude to Cycle 24; the maximum will occur no earlier than 2023
and no later than 2026 with the sunspot number from 95 to 130. However, at the
beginning of this year, Leamon and McIntosh announced that Cycle 24 had
finally ended (Termination Event occurred), and this allows us to look at the
prediction of macint20 differently:
https://spaceweatherarchive.com/2022/02/25/the-termination-event-has-arrived.
According to the new prediction, the maximum of the Cycle 25 will be about 190
(140–240), i.e., close in amplitude to the rather large Cycle 23.

Most recently, an article by brajsa22 was published, in which the maximum of
the started Cycle 25 is predicted as $121\pm 33$ — i.e., low.

Thus, at the moment there is uncertainty about the future amplitude of the
Cycle 25, and the scatter of predictions is large.

##  3 Cycle precursor associated with the minimum phase

In the article by brajsa22, it is proposed to choose values of $S\\!N$ three
years before the minimum of the cycle as a precursor of the maximum amplitude
of the cycle. The maximum smoothed monthly averages were used. We will use
annual averages in our work — this is not a fundamental difference for the
prediction. The interval for the study is the same as that of brajsa22 — from
the year 1749 to the present.

Let’s make some remarks. According to SILSO data, Cycle 20 has a fairly flat
maximum with maximum amplitude in the year 1968. At the same time, $S\\!N$
show a different maximum according to Kislovodsk data
http://www.solarstation.ru/archiv, the sunspot group number $G\\!N$ according
to svalshat16, usos16 data, and the sunspot area according to Greenwich data
https://solarscience.msfc.nasa.gov/greenwch.shtml: in 1970 — see Table 1, the
values in the maxima years are underlined. Based on this, in our work, we will
assume that the maximum of the 20th cycle corresponds to the middle of the
year 1970–1970.5.

Table 1:  Average annual values of various indices of solar activity in Cycle 20.  Year |  $S\\!N$, |  $S\\!N$, |  $G\\!N$, |  $G\\!N$, | Area,
---|---|---|---|---|---
| SILSO | Kislovodsk | Svalgaard | Usoskin |  $\mu$sh,
|  |  | & Schatten, 2016 | el al., 2016 | Greenwich
1964.5 | 15.0 | 12.9 | 0.88 | 0.845 | 54
1965.5 | 22.0 | 20.4 | 1.21 | 1.203 | 113
1966.5 | 66.8 | 61.2 | 3.32 | 3.635 | 593
1967.5 | 132.9 | 133.3 | 7.07 | 7.934 | 1519
1968.5 | 150.0 | 142.2 | 7.14 | 8.129 | 1570
1969.5 | 149.4 | 139.3 | 7.27 | 7.947 | 1450
1970.5 | 148.0 | 148.3 | 8.04 | 8.968 | 1601
1971.5 | 94.4 | 109.6 | 5.49 | 6.081 | 990
1972.5 | 97.6 | 108.8 | 5.33 | 5.955 | 917
1973.5 | 54.1 | 54.0 | 2.95 | 3.235 | 458
1974.5 | 49.2 | 46.7 | 2.69 | 2.824 | 399
1975.5 | 22.5 | 19.9 | 1.19 | 1.252 | 166
1976.5 | 18.4 | 16.7 | 1.08 | 1.122 | 170

Let’s take the values of $S\\!N$ three years before the minimum and calculate
their correlations with the values of different years that are separated by
$\Delta$ years from the next maximum. Negative values of $\Delta$ will
correspond to the years before the maximum, positive — after. Zero is the
maximum itself. The minimum will be denoted by $m$, the maximum by $M$. The
parity of the cycle will be determined by the predicted maximum. We will build
linear relations of the following form using the least squares method (LSM):

$S\\!N_{M+\Delta}=a+b\,S\\!N_{m-3}\,,$ (1)

and weakly non-linear ones:

$S\\!N_{M+\Delta}=a+b\,S\\!N_{m-3}+c\,S\\!N_{m-3}^{2}\,,$ (2)

and also calculate the coefficients of determination for them
$DC\equiv\rho^{2}$. This coefficient is preferable because it has a clear
meaning: it is the proportion of dispersion that can be described using the
appropriate regression. The results are shown in Figure 2 separately for even,
odd, and all cycles together. The dashed line shows the upper limit of the
“acceptable” level of correlation for the prediction from our point of view:
$DC=0.5$ ($\rho=0.7$).

Figure 2:  Coefficients of the determination for relations $S\\!N_{m-3}$ with
$S\\!N_{M+\Delta}$ (a) for the linear form (1), (b) for the nonlinear form
(2). The parity of the cycles and their mixture corresponds to the legend.

The first thing that can be seen in Figure 2: taking into account the
nonlinearity does not lead to a noticeable improvement in correlations. But
that’s not the main thing. It turns out that the tightness of the relation
$S\\!N_{m-3}$ with $S\\!N$ in the maximum and the years closest to it strongly
depends on the parity of the predicted cycle. We can predict an even cycle
with a coefficient of determination D = 0.882 (respectively, with
$\rho=0.94$), while an odd cycle — only with $DC=0.524$ ($\rho=0.72$).

According to the prediction in the style of the article by brajsa22, i.e., if
we do not take into account the parity of cycles, we get for a maximum of
Cycle 25: $S\\!N_{M}=129\pm 34$, $\rho=0.808$, which is quite consistent with
the values obtained in their article. According to the nonlinear prediction
(2), our same values were obtained — taking into account the possible
nonlinearity does not give anything additional: $S\\!N_{M}=(71\pm 17)+(1.45\pm
0.23)S\\!N_{m-3}$. It is also close to the received by brajsa22.

Now let’s see what the regressions and the prediction will look like if we
take into account the parity of the cycles. For even cycles regression

$S\\!N_{M}=(50\pm 14)+(1.60\pm 0.18)S\\!N_{m-3}\,,\quad\rho=0.939\,,$ (3)

If our predicted cycle were even, then its maximum would be $S\\!N_{M}=113\pm
21$. Confidence intervals due to the high correlation coefficient are very
small. However, Cycle 25 is odd. We calculate the regression for odd cycles:

$S\\!N_{M}=(91\pm 29)+(1.32\pm 0.40)S\\!N_{m-3}\,,\quad\rho=0.724$ (4)

and the prediction itself:

$S\\!N_{M}=144\pm 44\,.$ (5)

This implies that brajsa22, having considered even and odd cycles together,
underestimated the predicted value of Cycle 25. The difference in regressions
for even and odd cycles is shown in Figure 3a.

Figure 3:  (a) Dependence of $S\\!N_{M}$ on $S\\!N_{m-3}$ for even (solid
line, gray squares) and odd (dashed line, black circles) cycles. (b)
Dependence of $S\\!N_{M-7}$ on $S\\!N$ for even and odd cycles.

Summarizing, according to the effect proposed by brajsa22, we can successfully
predict the maximum of the next EVEN cycle three years before the minimum,
thereby supplementing the Gnevyshev-Ohl rule with the relation even –
subsequent odd cycles. The success of the prediction of an even cycle is
guaranteed by a high correlation coefficient of the equation (3).

The precursor of the prediction discussed in this section will be called MI3E
for even cycles and MI3O for odd ones (bearing in mind that only MI3E
corresponds to a reliable prediction).

##  4 Cycle precursor associated with the maximum phase

In the paper by slonim84, the idea was expressed to search for cycle
precursors a certain number of years before the maximum. Thus, the possible
predictors turn out to be associated with the maximum phase, and not with the
minimum phase, as in the previous section. Here we will also consider
regressions separately for even, odd cycles and for all cycles together,
choosing for analysis the same time interval as before, since 1749. Figure 4a
illustrates the changes in the coefficients of determination for linear
relationships of the type

$S\\!N_{M}=a+b\cdot S\\!N_{M+\Delta}\,.$ (6)

Here we bear in mind that $\Delta=0$ corresponds to the maximum of the cycle.
Precursors $S\\!N_{M-\Delta}$ correspond to negative $\Delta$ — i.e., years
before maximum, positive ones indicate relations with post-maximum years,
which are also interesting for the prediction. In addition, we take into
account the weak nonlinearity of the relations. For negative $\Delta$ in the
form

$S\\!N_{M}=a+b\cdot S\\!N_{M+\Delta}+c\cdot S\\!N_{M+\Delta}^{2}\,,$ (7)

and for positive ones —

$S\\!N_{M+\Delta}=a+b\cdot S\\!N_{M}+c\cdot S\\!N_{M}^{2}\,.$ (8)

We do this specifically for the possible prediction of not only the maximum
values but also the index values in subsequent years, which, as it turns out,
are quite closely related to the maximum. The determination coefficients for
relations of type (7)–(8) are shown in Figure 4b.

Figure 4:  Coefficients of the determination for relations $S\\!N_{M}$ with
$S\\!N_{M+\Delta}$ for (a) linear form (6), (b) nonlinear forms (7)–(8).

Let’s note two circumstances. The first — in comparison with $S\\!N_{m-3}$ a
more “successful” precursor for the prediction of the maximum is
$S\\!N_{M-7}$: in the case of a linear form (6) for odd cycles $DC=0.849$
($\rho=0.921$), for even $DC=0.863$ ($\rho=0.929$); in the case of a nonlinear
form (7) for odd $DC=0.859$ ($\rho=0.927$), for even $DC=0.903$
($\rho=0.950$). Secondly, when considering all cycles together, the
correlation of $S\\!N_{M-7}$ with $S\\!N_{M}$ is less than for each of the
parity separately, and this means that their relations are different. This is
clearly seen in Figure 3b.

When we predict the amplitude of the future cycle, we are, of course,
primarily interested in the amplitude of the maximum, but we are also
interested in neighboring years. We found that the predictor of the maximum is
$S\\!N_{M-7}$. And how much can this value predict the years $M-1$, $M+1$,
$M+2$, etc.? Figure 5 answers this question: the seventh year before the
maximum is a predictor not only of the maximum but also of the index in the
years near the maximum. This applies to both even and odd cycles, although the
relations between $S\\!N_{M-7}$ and $S\\!N_{M+\Delta}$ are different.

Figure 5:  Determination coefficients for quadratic relations of $S\\!N_{M-7}$
with $S\\!N_{M+\Delta}$.

The prediction precursor discussed in this section we will call MA7E for even
cycles and MA7O for odd ones.

Thus, we have created a “constructor” for the prediction of the solar cycle,
in this case — the odd solar Cycle 25.

##  5 Prediction of the solar activity cycle 25

Let’s use the results of the previous sections to predict Cycle 25. We use the
MA7O predictor. Assuming different years of maximum onset, we get the
$S\\!N_{M}$ values given in the Table 2.

Table 2:  Predicted $S\\!N_{M}$ values depending on the estimated maximum year.  Year | $S\\!N_{M}$
---|---
2022.5 | $199\pm 25$
2023.5 | $159\pm 25$
2024.5 | $133\pm 25$

Then, using regressions, the coefficients of determination of which are shown
in Figure 5, we calculate the values of $S\\!N$ in the years of Cycle 25,
which have prognostic significance, for different maximum years from Table 2.

The results are shown in Figure 6a: light icons — rhombuses, squares, and
circles — connected by a dashed line constructed using a global cubic spline.
Comparing the obtained curves with the course of the monthly average values of
$S\\!N$ observed up to June 2022 (small circles connected by thin lines), we
conclude that 2023.5 is the most preferable of the three years of the assumed
maximum from Table 2. However, in general, the monthly average values are
located, although close to the interpolation curve, but still somewhat below
her. Let’s assume a new maximum date — 2024.0. Note that the average annual
values can be calculated by averaging $S\\!N$ not only traditionally from
January to December but also from July to June of the next year, and. Let’s
calculate such annual averages, output regressions with determination
coefficients similar to those shown in Figure 5, and calculate the values for
the dates 2023.0–2028.0 in annual increments — see Figure 6a: black circles
connected by a thick spline curve. It can be seen that by assuming a maximum
year of 2024.0 we get a better agreement of the behavior of the predicted and
observed $S\\!N$.

Now about the values indicated in Figure 6a by crosses. Are there any general
trends in the last years of the cycle (in our case, the odd one)? Consider the
values of $S\\!N$ in the years of the final minima of the odd cycle $m$ and
the two years before it $m-1$ and $m-2$, depending on their distance from the
initial minimum $\Delta$. It turns out that this dependence can be described
by a parabola $S\\!N(\Delta)=(327\pm 87)-(50\pm 18)\Delta+(1.96\pm
0.87)\Delta^{2}$ with a correlation coefficient $\rho=0.811$, i.e.
sufficiently high, and we can use it to predict the descending phase of a new
cycle (as it seems, this is a new unexpected result). In Figure 6a crosses
represent data given by the date of the initial minimum of Cycle 25 as
$2019.5\pm\Delta$. Note here that for even cycles, a similar dependence also
occurs, although the correlation coefficient with the parabola is somewhat
less — $\rho=0.768$ . However, it can also be used for solar cycle prediction
in the future.

So far, for predictions, we have used confidence ranges of estimates for a
deviation equal to the standard deviation (i.e., with a probability of falling
into the interval equal to 68%) — as in brajsa22, as well as in many authors
of solar cycle predictions. In fact, if we want to achieve a reliability of at
least 95%, we need to set a prediction interval based on at least two standard
deviations. Therefore, we will repeat all the previous procedures for this new
requirement. The average annual values of $S\\!N$, interpreted for years with
traditional calculation from January to December, with confidence bands are
shown in Figure 6b.

Note that our methodology allows us to predict individual average annual
values of $S\\!N$ in a cycle, and not only the amplitude of its maximum, as
with most authors.

Having made a prediction of the Cycle 25, it is natural to check whether the
Gnevyshev-Ohl rule is fulfilled for a pair of Cycles 24–25. In Figure 1a, the
corresponding point is circled. As we can see, our prediction does not
contradict this well-known rule.

Figure 6:  (a) Various versions of predictions for the development of Cycle 25
— see the legend in the text; (b) The most probable prediction of Cycle 25.
Light gray areas are the prediction interval with a confidence of 68%, dark
gray areas are the same with a confidence of 95%, and a stepped line is the
most probable prediction of average annual values.

1

##  6 Results

The Gnevyshev-Ohl rule establishes a relation between adjacent 11-year cycles
in the 22-year “magnetic” Hale cycle as: “an even cycle determines the next
odd one,” and thus they form a pair. In this work, we found that the behavior
of an odd cycle determines the behavior of the subsequent even one. Namely, 3
years before the minimum, the value of $S\\!N$ in the odd cycle is associated
with the value of the maximum in the subsequent even cycle (correlation
coefficient $\rho=0.94$. Similar to GOR, in this sense, cycles are linked in
pairs, but opposite to the Rule. This is one of the important results of the
work.

Based on this result, we propose to use $S\\!N_{m-3}$ (Figure 3a) as a
precursor of the subsequent EVEN cycle during the descending phase of the odd
one — we call this method MI3E. For the prediction of an odd cycle or a
prediction without parity (as in the article by brajsa22), this method gives
less reliable results.

To predict the ODD cycle, we propose to use the precursor of the seventh year
to its maximum MA7O — $S\\!N_{M-7}$ during the descending phase of the
previous even cycle (Figure 3b). It turned out that in this case, we can
predict the years near the maximum with a high correlation coefficient
($\rho=0.90{-}0.94$). In addition, 7 years before the maximum, it is also
possible to predict an even cycle according to the MA7E precursor (Figure 5).

Also noteworthy is the similar behavior of $S\\!N$ for different cycles near
the final minimum of the cycle, depending on the distance from the initial
minimum. Thus, the proposed approaches make it possible to predict cycles of
different parity.

The current Cycle 25 in the year 2023 should reach a maximum of 154 average
annual values with a prediction interval of 25 with a confidence of 68% and 53
with a confidence of 95%. This is higher than the original official prediction
of NOAA/NASA/ISES but lower than the updated prediction of leamon21 according
to the methodology of macint20. In 2024, $S\\!N$ will be almost as high as in
2023 — 147 units, so with smaller time averaging scales, the maximum will fall
at the end of 2023. Here we note that the proposed approach makes it possible
to predict the values of $S\\!N$ in individual years of the cycle, and not
only the amplitude of its maximum.

We have not discussed in this article the relationship between the precursors
$S\\!N_{m-3}$ and $S\\!N_{M-7}$, their relative position on the descending
phase of cycles, as well as their physical meaning — this is the task of the
following studies.

Yury Nagovitsyn thanks the Ministry of Science and Higher Education of the
Russian Federation for financial support for project 075–15–2020–780.

## References

  * Brajša et al. (2022) Brajša, R.; Verbanac, G.; Bandić, M.; Hanslmeier, A.; Skokić, I.; Sudar, D.: 2022, Astron. Nachr., 343(3), article id. e13960. DOI: 10.1002/asna.202113960
  * Clette et al. (2014) Clette, F.; Svalgaard, L.; Vaquero, J.M.; Cliver, E.W.: 2014, Space Sci. Rev., 186(1-4), 35. DOI: 10.1007/s11214-014-0074-2
  * Clette et al. (2016) Clette, F.; Cliver, E.W.; Lefèvre, L.; Svalgaard, L.; Vaquero, J.M.; Leibacher, J.W.: 2016, Solar Phys., 291(9-10), 2479. DOI: 10.1007/s11207-016-1017-8
  * Gnevyshev and Ohl (1948) Gnevyshev, M. N.; Ohl, A. I.: 1948, Astron. Zh., 25, 18.
  * Leamon et al. (2021) Leamon, R.J., McIntosh, S.W., Chapman, S.C., Watkins, N.W.: 2021, Sol. Phys., v. 296(10), 151. DOI: 10.1007/s11207-021-01897-z
  * McIntosh et al. (2020) McIntosh, S.W.; Chapman, Sandra, Leamon, Robert J., Egeland, Ricky, Watkins Nicholas W.: 2020, Solar Phys, 295(12), 163. DOI: 10.1007/s11207-020-01723-y
  * Nagovitsyn (2009) Nagovitsyn Yu. A.: 2009, Astron. Lett., , 35(8), 564. DOI: 10.1134/S1063773709080064
  * Ohl (1966) Ohl, A. I.: 1966, Byulletin Solnechnye Dannye Akademie Nauk USSR, 12, 84
  * Pesnell (2014) Pesnell, W.D.: 2014, Solar Phys. 289(6), 2317. DOI: 10.1007/s11207-013-0470-x
  * Petrovay (2020) Petrovay, K.: 2020, Liv. Rev. Sol. Phys., 17(1), 2. DOI: 10.1007/s41116-020-0022-z
  * Slonim (1984) Slonim Yu.M. 1984: Byulletin Solnechnye Dannye Akademie Nauk USSR, 5, 78\. ADS: https://ui.adsabs.harvard.edu/abs/1984BSolD1984Q..78S/abstract
  * Svalgaard el al. (2005) Svalgaard, L.; Cliver, E.W.; Kamide, Y.: 2005, Geophysical Research Letters, 32, L01104.
  * Svalgaard and Schatten (2016) Svalgaard, L.; Schatten, K.H.: 2016, Solar Phys., 291(9-10), 2653. DOI: 10.1007/s11207-015-0815-8
  * Usoskin et al. (2016) Usoskin, I.G.; Kovaltsov, G.A.; Lockwood, M.; Mursula, K.; Owens, M.; Solanki, S.K.: 2016, Solar Phys., 291(9-10), 2685. DOI: 10.1007/s11207-015-0838-1
  * Vitinsky (1973) Vitinsky, Yu.I., Tsiklichnost’ i prognozy solnechnoi aktivnosti (Cyclicity and Forecasting of Solar Activity), Moscow: Nauka, 1973 [in Russian].
  * Wilson et al. (1998) Wilson, R. M.; Hathaway, D. H.; Reichmann, E. J.: 1998, Journ. Geophys. Res., 103(A4), 6595. DOI: 10.1029/97JA02777


# Superconvergent P1 honeycomb virtual elements and lifted P3 solutions

Yanping Lin  Department of Applied Mathematics, The Hong Kong Polytechnic
University, Hung Hom, Hong Kong<EMAIL_ADDRESS>,  Xuejun Xu  School
of Mathematical Science, Tongji University, Shanghai, 200092, China  Institute
of Computational Mathematics, AMSS, Chinese Academy of Sciences, Beijing,
100190, China<EMAIL_ADDRESS>and  Shangyou Zhang  Department of
Mathematical Sciences, University of Delaware, Newark, DE 19716, USA.
<EMAIL_ADDRESS>
###### Abstract.

When solving the Poisson equation on honeycomb hexagonal grids, we show that
the $P_{1}$ virtual element is three-order superconvergent in $H^{1}$-norm,
and two-order superconvergent in $L^{2}$ and $L^{\infty}$ norms. We define a
local post-process which lifts the superconvergent $P_{1}$ solution to a
$P_{3}$ solution of the optimal-order approximation. The theory is confirmed
by a numerical test.

###### Key words and phrases:

virtual element, honeycomb grid, superconvergence, Poisson equation.

###### 1991 Mathematics Subject Classification:

65N15, 65N30.

Yanping Lin is supported in part by HKSAR GRF 15302922 and polyu-CAS joint
Lab.

Xuejun Xu is supported by National Natural Science Foundation of China (Grant
Nos. 12071350), Shanghai Municipal Science and Technology Major Project No.
2021SHZDZX0100, and Science and Technology Commission of Shanghai
Municipality.

##  1\. Introduction

In this work, we study the $P_{1}$ virtual element on honeycomb meshes for
solving the following Poisson equation.

(1.1) $\displaystyle\begin{aligned} -\Delta u&=f&&\hbox{\rm{in }}\ \Omega,\\\
u&=0&&\hbox{\rm{on }}\ \partial\Omega,\end{aligned}$

where $f\in L^{2}(\Omega)$ and $\Omega\subset\mathbb{R}^{2}$ is a bounded
polygonal domain which can be meshed by one-size honeycomb hexagons and a few
additional polygons consisting of equilateral triangles only at domain
boundary, such as the one illustrated in Figure 1.

Figure 1.  The first three honeycomb grids on a regular hexagon domain
$\Omega$.

The virtual element method is proposed and studied in [2, 3, 9, 10, 12, 13,
14, 16, 17, 18, 19, 20, 31, 35]. In 2D, the virtual element functions are
continuous $P_{k}$ polynomials on edges of a polygonal mesh and extended to
inside polygons by the Poisson equation with $P_{k-2}$ moments. In particular,
for the 2D $P_{1}$ virtual elements studied in this paper, the continuous
$P_{1}$ polynomials on the edges of the mesh are harmonically extended into
polygons. Such non-polynomial virtual element functions are interpolated into
a computable polynomial space on each polygon. As the interpolation is usually
not a full-rank operator, a stabilizer has to be added to the discrete
equation. However, simply raising the polynomial degree of the local
interpolation space, in order to eliminate the stabilizer, does not work in
general, cf. [36].

The references [4, 5] propose to use some high-degree polynomial interpolation
space only for 2D $P_{1}$ virtual elements to get some stabilizer-free virtual
element methods. The reference [6] proposes to use only high-$P_{k}$ harmonic
polynomials as the interpolation functions in a stabilizer-free virtual
element method. The references [24, 36] use multi-piece polynomials as the
interpolation space on each polygon or polyhedron, so that the resulting
virtual element method works for any $P_{k}$ virtual element without using
stabilizer.

In this work, we apply the stabilizer-free $P_{1}$ virtual element method [24,
36] to the Poisson equation (1.1) on 2D honeycomb meshes, cf. Figure 1. Taking
the advance of extreme symmetry, we show that the $P_{1}$ virtual element
solution $u_{h}$ is three-order superconvergent (three orders above the
optimal order) in $H^{1}$-norm, i.e.,

(1.2) $\displaystyle\|\nabla(I_{h}u-u_{h})\|_{L^{2}}\leq
Ch^{4}|u|_{H^{4,\infty}},$

where $H^{4,\infty}$ is the Sobolev space of functions with all 4-th
derivatives bounded, and $I_{h}$ is the $P_{1}$ nodal interpolation at corners
of honeycombs. Additionally, the $P_{1}$ virtual element solution is two-order
superconvergent in $L^{2}$-norm and in $L^{\infty}$ norm, i.e.,

(1.3) $\displaystyle\|I_{h}u-u_{h}\|_{L^{2}}\leq
C\|I_{h}u-u_{h}\|_{L^{\infty}}\leq Ch^{4}|u|_{H^{4,\infty}}.$

We then define a local-lift post-process, which lifts the $P_{1}$ virtual
element solution $u_{h}$ to a piecewise $P_{3}$ solution $\tilde{u}_{h}$ of
the optimal order approximation, i.e.,

(1.4)
$\displaystyle\|u-\tilde{u}_{h}\|_{L^{2}}+h|u-\tilde{u}_{h}|_{H^{1},h}\leq
Ch^{4}|u|_{H^{4,\infty}},$

where $|\cdot|_{H^{1},h}$ is a piecewise $H^{1}$-semi-norm based on a patch
grid $\mathcal{R}_{h}$ (defined in section 3.) We remark that the $P_{1}$
virtual element solution itself can only converge at the optimal order, two
orders lower than that of (1.4),

$\displaystyle\|u-u_{h}\|_{L^{2}}+h|u-u_{h}|_{H^{1},h}\leq Ch^{2}|u|_{H^{2}}.$

The theory is confirmed by a numerical test.

Some early work and influential references on the superconvergence of $P_{1}$
finite elements can be found in [7, 8, 11, 15, 21, 23, 27, 29, 30, 32, 48, 49,
50]. Some superconvergence results for weak Galerkin and conforming
discontinuous Galerkin finite element methods can be found in [1, 33, 34, 37,
38, 39, 40, 41, 42, 43, 44, 45, 46, 47].

##  2\. The $P_{1}$ honeycomb virtual element

We define the $P_{1}$ honeycomb virtual element in this section and show that
the $P_{1}$ honeycomb virtual element equation has a unique solution.

${\bf x}_{0}$${\bf x}_{6}$${\bf x}_{1}$${\bf x}_{2}$${\bf x}_{3}$${\bf
x}_{5}$${\bf x}_{4}$${\bf x}_{3}$${\bf x}_{1}$${\bf x}_{2}$${\bf x}_{0}$${\bf
x}_{2}$${\bf x}_{3}$${\bf x}_{4}$${\bf x}_{1}$ Figure 2. Left: A regular
hexagon inside $\Omega$, where ${\bf x}_{0}$ is an auxiliary vertex of
triangular mesh $\overline{\mathcal{T}}_{h}$; Middle: A regular triangle which
can be only at a corner of $\Omega$; Right: A pentagon consisting of three
regular triangles which can be at an edge of $\Omega$ with ${\bf x}_{1}{\bf
x}_{4}\subset\partial\Omega$.

Let $\mathcal{T}_{h}=\\{T\\}$ be a uniform honeycomb mesh on the domain
$\Omega$, where $T$ is either a regular hexagon of edge-size $h$ inside
$\Omega$, or a regular triangle of edge-side $h$ at a corner of the polygonal
domain, or a pentagon consisting of three regular triangles of edge-side $h$
at an edge of domain, cf. Figure 2. We note that there is a vertex added at
the middle of straight edge ${\bf x}_{1}{\bf x}_{4}$, cf. Figure 2. For each
hexagon $T$, we add the central point ${\bf x}_{0}$ to it so that it is
subdivided into 6 regular triangles $\\{T_{i}\\}$, cf. the left diagram of
Figure 2. For a pentagon $T\in\mathcal{T}_{h}$, we connect its center point
${\bf x}_{0}$ with two vertices ${\bf x}_{2}$ and ${\bf x}_{3}$ so that it is
subdivided into 3 regular triangles, cf. the right diagram in Figure 2. In
addition to the regular triangles in $\mathcal{T}_{h}$ (the middle graph in
Figure 2), all theses regular triangles form a equilateral triangular mesh
$\overline{\mathcal{T}}_{h}$ on domain $\Omega$.

Let $\mathcal{E}_{h}$ denote the set of edges $e$ in $\mathcal{T}_{h}$. Let
$E=\cup_{e\in\mathcal{E}_{h}}e$ and $\mathbb{B}_{1}(E)$ be the space of
piecewise linear functions on domain $E$,

$\displaystyle\mathbb{B}_{1}(E)=\\{\tilde{v}\in C^{0}(E):\tilde{v}|_{e}\in
P_{1}(e)\ \hbox{\rm{ for all }}e\in\mathcal{E}_{h}\\}.$

The virtual element space on the mesh $\mathcal{T}_{h}$ is defined as

(2.1) $\displaystyle\begin{aligned} \tilde{V}_{h}=\\{\tilde{v}\in
H^{1}_{0}(\Omega)&:\tilde{v}\in\mathbb{B}_{k}(\mathcal{E}_{h}),\
\Delta\tilde{v}|_{K}=0\quad\forall T\in\mathcal{T}_{h}\\}.\end{aligned}$

In computation, the non-polynomial virtual element function is interpolated to
a piecewise polynomial. The piecewise polynomial space on each
$T\in\mathcal{T}$ is defined by, cf. Figure 2,

$\displaystyle\mathbb{V}_{1}(T)=\begin{cases}\\{v_{h}\in
H^{1}(T):v_{h}|_{T_{i}}\in P_{1}(T_{i})\ \hbox{\rm{if
$T=\bigcup_{i=1}^{6}T_{i}$ is a hexagon}}\\},\\\ \\{v_{h}\in
H^{1}(T):v_{h}|_{T_{1}}\in P_{1}(T_{1})\ \hbox{\rm{if $T=T_{1}$ is a
triangle}}\\},\\\ \\{v_{h}\in H^{1}(T):v_{h}|_{T_{i}}\in P_{1}(T_{i})\
\hbox{\rm{if $T=\bigcup_{i=1}^{4}T_{i}$ is a pentagon}}\\},\end{cases}$

where $T_{i}$ is an edge-size $h$ regular triangle.

The interpolated virtual element space on $\mathcal{T}_{h}$ is defined by

(2.2) $\displaystyle V_{h}=\\{v_{h}=\Pi_{h}^{\nabla}\tilde{v}\ :\
v_{h}|_{T}\in\mathbb{V}_{1}(T),\ T\in\mathcal{T}_{h};\
\tilde{v}\in\tilde{V}_{h}\\},$

where $v_{h}=\Pi_{h}^{\nabla}\tilde{v}$ is the local $H^{1}$-projection, i.e.,
$v_{h}$ is the local solution on $T\in\mathcal{T}_{h}$:

$\displaystyle\left\\{\begin{aligned} (\nabla(v_{h}-\tilde{v}),\nabla
w_{h})_{T}&=0\quad\forall
w_{h}=\Pi_{h}^{\nabla}\tilde{w}\in\mathbb{V}_{1}(T),\\\ \langle
v_{h}-\tilde{v},w_{h}\rangle_{\partial T}&=0\quad\forall
w_{h}\in\mathbb{V}_{1}(T).\end{aligned}\right.$

The stabilizer-free virtual element equation reads: Find
$u_{h}=\Pi_{h}^{\nabla}\tilde{u}\in V_{h}$ such that

(2.3) $\displaystyle(\nabla u_{h},\nabla
v_{h})=(f,v_{h})\quad\forall\tilde{v}\in\tilde{V}_{h},\
v_{h}=\Pi_{h}^{\nabla}\tilde{v}.$

######  Lemma 2.1.

The $P_{1}$ honeycomb virtual element equation (2.3) has a unique solution.

###### Proof.

(2.3) is a finite dimensional square system of linear equations. Its existence
of solution is implied by the uniqueness of the solution. Let $f=0$ and
$v_{h}=u_{h}$ in (2.3), where $u_{h}$ is a solution. (2.3) becomes $\|\nabla
u_{h}\|_{T}^{2}=0$ on every $T\in\mathcal{T}_{h}$. Thus $u_{h}\in P_{0}(T)$.
Because $u_{h}\in C^{0}(\Omega)$, $u_{h}=C$ is a global constant on $\Omega$.
By $u_{h}|_{\partial\Omega}=0$, $u_{h}=0$. The proof is complete. ∎

##  3\. The superconvergence in $H^{1}$-norm

We show in this section that the $P_{1}$ honeycomb virtual element is three-
order superconvergent in $H^{1}$-norm first, and is two-order superconvergent
in $L^{2}$-norm and in $L^{\infty}$-norm.

######  Theorem 3.1.

Let $u\in H^{4,\infty}(\Omega)\cap H^{1}_{0}(\Omega)$ and $u_{h}\in V_{h}$ of
(2.2) be the solutions of Poisson’s equation (1.1) and the finite element
equation (2.3), respectively. It holds that

$\displaystyle|I_{h}u-u_{h}|_{H^{1}}\leq Ch^{4}|u|_{H^{4,\infty}}.$

###### Proof.

For analysis, we introduce the auxiliary, conforming $P_{1}$ finite element on
a sub-triangular mesh of $\mathcal{T}_{h}$. Let $\mathcal{T}_{h}^{6}$ be the
set of all hexagons in $\mathcal{T}_{h}$. For each $T\in\mathcal{T}_{h}^{6}$,
we introduce a new vertex/node ${\bf x}_{0}$, the center of the regular
hexagon. Let $\mathcal{N}_{h}$ be the set of all vertices/nodes in
$\mathcal{T}_{h}$. Let

(3.1) $\displaystyle\overline{\mathcal{N}}_{h}=\mathcal{N}_{h}\bigcup\\{{\bf
x}_{0}:{\bf x}_{0}\in T\ \forall T\in\mathcal{T}_{h}^{6}\\}.$

Subdividing all hexagons and pentagons in $\mathcal{T}_{h}$ by connecting
${\bf x}_{0}$ with all ${\bf x}_{i}$, cf. Figure 2, we have an equilateral
triangular mesh on $\Omega$,

(3.2) $\displaystyle\overline{\mathcal{T}}_{h}=\\{T:T\ \hbox{\rm{ is a regular
triangle of edge-size $h$}}\\}.$

The auxiliary, conforming $P_{1}$ finite element space on triangular mesh
$\overline{\mathcal{T}}_{h}$ is defined by

(3.3) $\displaystyle\overline{V}_{h}=\\{v_{h}\in
H^{1}_{0}(\Omega):v_{h}|_{T}\in P_{1}(T)\ \forall
T\in\overline{\mathcal{T}}_{h}\\}.$

Let $e_{h}=I_{h}u-u_{h}$, where $I_{h}:H^{2}(\Omega)\cap H^{1}_{0}(\Omega)\to
V_{h}$ is the nodal interpolation operator based on nodal set
$\mathcal{N}_{h}$. Let $\overline{I}_{h}:H^{2}(\Omega)\cap
H^{1}_{0}(\Omega)\to\overline{V}_{h}$ (defined in (3.3)) be the nodal
interpolation operator based on nodal set $\overline{\mathcal{N}}_{h}$
(defined in (3.1)), i.e.,

$\displaystyle I_{u}u({\bf x}_{j})=u({\bf x}_{j})\quad\forall{\bf
x}_{j}\in\overline{\mathcal{N}}_{h}.$

Let

(3.4) $\displaystyle d_{h}=I_{h}u-\overline{I}_{h}u\in\overline{V}_{h}.$

By the harmonic solution in (2.2), $d_{h}$ is supported only inside all
hexagons $T\in\mathcal{T}^{6}_{h}$ and, cf Figure 2, due to a regular hexagon,

$\displaystyle d_{h}({\bf x}_{0})$
$\displaystyle=\frac{1}{6}\sum_{i=1}^{6}u({\bf x}_{i})-u({\bf x}_{0}),$
$\displaystyle d_{h}({\bf x}_{i})$ $\displaystyle=0,\quad i=1,\dots,6.$

Because

(3.5) $\displaystyle(\nabla u-\nabla u_{h},\nabla v_{h})=0\quad\forall
v_{h}\in V_{h},$

it follows that by the divergence theorem,

(3.6) $\displaystyle\begin{aligned}
|I_{h}u-u_{h}|_{H^{1}}^{2}&=(\nabla(I_{h}u-u),\nabla e_{h})\\\ &=(\nabla
d_{h},\nabla e_{h})+(\nabla(\overline{I}_{h}u-u),\nabla e_{h})\\\
&=\sum_{T\in\mathcal{T}_{h}^{6}}(\nabla d_{h},\nabla
e_{h})_{T}+\sum_{T\in\overline{\mathcal{T}}_{h}}(\nabla(\overline{I}_{h}u-u),\nabla
e_{h})_{T}\\\
&=0+\sum_{T\in\overline{\mathcal{T}}_{h}}\sum_{i=0}^{2}\int_{{\bf x}_{i}{\bf
x}_{i+1}}(\overline{I}_{h}u-u)\partial_{{\bf t}_{i+2}^{\perp}}e_{h}d{\bf
s},\end{aligned}$

where ${\bf t}_{i}^{\perp}$ is the $90$ degrees clockwise rotation of ${\bf
t}_{i}$, and ${\bf x}_{i}$ are three vertices of $T$ depicted in Figure 3.
Here $(\nabla d_{h},\nabla e_{h})_{T}=0$ because $e_{h}$ is discrete harmonic
inside a hexagon $T$ by (2.2), and $d_{h}$ is supported inside $T$.

${\bf x}_{1}$${\bf x}_{2}$${\bf x}_{0}={\bf x}_{3}$$T:$${\bf t}_{2}$${\bf
t}_{1}$${\bf t}_{3}={\bf t}_{0}$ Figure 3. A equilateral triangle with three
unit tangent vectors.

For linear polynomials, the normal derivative is a linear combination of two
tangential derivatives. We have, due to equilateral triangles, cf. Figure 3,

(3.7) $\displaystyle\begin{aligned} &\quad\
\sum_{T\in\overline{\mathcal{T}}_{h}}\sum_{i=0}^{2}\int_{{\bf x}_{i}{\bf
x}_{i+1}}(\overline{I}_{h}u-u)\partial_{{\bf t}_{i-1}^{\perp}}e_{h}d{\bf s}\\\
&=\sum_{T\in\overline{\mathcal{T}}_{h}}\sum_{i=0}^{2}\int_{{\bf x}_{i}{\bf
x}_{i+1}}(\overline{I}_{h}u-u)\frac{2}{\sqrt{3}}(\partial_{{\bf
t}_{i+1}}-\frac{1}{2}\partial_{{\bf t}_{i-1}})e_{h}d{\bf s}\\\
&=\sum_{T\in\overline{\mathcal{T}}_{h}}\sum_{i=0}^{2}\int_{{\bf x}_{i}{\bf
x}_{i+1}}(\overline{I}_{h}u-u)\frac{2}{\sqrt{3}}\partial_{{\bf
t}_{i+1}}e_{h}d{\bf s},\end{aligned}$

where the second term cancels itself from the other side of triangle because
${\bf t}_{i+1}$ has an opposite direction on the other side of edge ${\bf
x}_{i}{\bf x}_{i+1}$, or $e_{h}=0$ if ${\bf x}_{i}{\bf x}_{i+1}$ has no other
neighbor triangle, i.e., a boundary edge.

A version of the Euler-Maclaurin formula is, cf. [22],

$\displaystyle\int_{a}^{b}f(t)dt$
$\displaystyle=\sum_{n=a}^{b}f(n)-\frac{f(a)+f(b)}{2}-\sum_{i=1}^{k}\frac{b_{i}}{i!}(f^{(i-1)}(b)-f^{(i-1)}(a))$
$\displaystyle\quad\ -\int_{a}^{b}P_{k+1}(t)f^{(k+1)}(t)dt,$

where $b_{i}$ is a Bernoulli number and $P_{k+1}(t)$ is a scaled periodic
Bernoulli polynomial $B_{k+1}(\\{1-t\\})/((k+1)!)$. Other than $b_{1}$, all
odd Bernoulli numbers are zero. Letting $k=3$, $a=0$ and $b=1$, we have

$\displaystyle\int_{0}^{1}fdt$
$\displaystyle=h\frac{f(0)+f(1)}{2}+\frac{h^{2}}{4}\int_{0}^{1}f^{\prime\prime}dt-h^{4}\int_{0}^{1}P_{4}(t)f^{(4)}(t)dt.$

Noting that $\int_{{\bf x}_{1}{\bf x}_{2}}\overline{I}_{h}ud{\bf s}=(u({\bf
x}_{1})+u({\bf x}_{2}))h/2$, we get from (3.6) and (3.7),

(3.8) $\displaystyle\begin{aligned}
|e_{h}|_{H^{1}}^{2}&=-\sum_{T\in\overline{\mathcal{T}}_{h}}\sum_{i=0}^{2}\frac{h^{2}}{2\sqrt{3}}\int_{{\bf
x}_{i}{\bf x}_{i+1}}\partial_{{\bf t}_{i-1}}^{2}u\partial_{{\bf
t}_{i+1}}e_{h}d{\bf s}\\\ &\quad\
+\sum_{T\in\overline{\mathcal{T}}_{h}}\sum_{i=0}^{2}h^{4}\int_{{\bf x}_{i}{\bf
x}_{i+1}}P_{4}({\bf x})\partial^{4}_{{\bf t}_{0}}u\partial_{{\bf
t}_{i+1}}e_{h}d{\bf s},\\\
&=:\frac{h^{2}}{2\sqrt{3}}I_{1}+I_{2}.\end{aligned}$

By Cauchy-Schwarz inequality, we get, because $P_{4}$ has a zero-point inside
the edge ($\int_{{\bf x}_{i}{\bf x}_{i+1}}P_{4}({\bf x})d{\bf s}=0$),

(3.9) $\displaystyle\begin{aligned} |I_{2}|&\leq
Ch^{4}|u|_{H^{4,\infty}}\bigg{(}\sum_{T\in\overline{\mathcal{T}}_{h}}|e_{h}|_{H^{1}(T)}^{2}\bigg{)}^{1/2}=Ch^{4}|u|_{H^{4,\infty}}|e_{h}|_{H^{1}}.\end{aligned}$

For the first term in (3.8), we consider an area-integral, cf. Figure 3, by
integration by parts, noting that $\partial_{{\bf t}_{2}}e_{h}$ is a constant,

$\displaystyle\quad\ \int_{T}(\partial_{{\bf t}_{1}}\partial_{{\bf
t}_{0}}^{2}u)\partial_{{\bf t}_{2}}e_{h}d{\bf x}=\int_{{\bf x}_{1}{\bf
x}_{0}}\partial_{{\bf t}_{0}}^{2}u\partial_{{\bf t}_{2}}e_{h}d{\bf
s}-\int_{{\bf x}_{1}{\bf x}_{2}}\partial_{{\bf t}_{0}}^{2}u\partial_{{\bf
t}_{2}}e_{h}d{\bf s}.$

When summing the first term above, as ${\bf t}_{2}$ is tangential to ${\bf
x}_{1}{\bf x}_{0}$, the term is canceled with the term on the other side of
the triangle. Thus, we have, from (3.8),

$\displaystyle I_{1}$
$\displaystyle=\sum_{T\in\overline{\mathcal{T}}_{h}}\sum_{i=0}^{2}\int_{T}(\partial_{{\bf
t}_{i}}\partial_{{\bf t}_{i-1}}^{2}u)\partial_{{\bf t}_{i+1}}e_{h}d{\bf x}.$

Taking integration by parts, we get, noting that the angle between two ${\bf
t}_{i}$ is $120^{0}$,

(3.10) $\displaystyle\begin{aligned}
I_{1}&=\sum_{T\in\overline{\mathcal{T}}_{h}}\sum_{i=0}^{2}\bigg{(}-\int_{T}(\partial_{{\bf
t}_{i}}\partial_{{\bf t}_{i+1}}\partial_{{\bf t}_{i-1}}^{2}u)e_{h}d{\bf x}\\\
&\ \quad+\frac{1}{2}\int_{{\bf x}_{i}{\bf x}_{i+1}}(\partial_{{\bf
t}_{i}}\partial_{{\bf t}_{i-1}}^{2}u)e_{h}d{\bf s}+\frac{1}{2}\int_{{\bf
x}_{i+1}{\bf x}_{i-1}}(\partial_{{\bf t}_{i}}\partial_{{\bf
t}_{i-1}}^{2}u)e_{h}d{\bf s}\bigg{)}\\\
&=\sum_{T\in\overline{\mathcal{T}}_{h}}\sum_{i=0}^{2}-\int_{T}(\partial_{{\bf
t}_{i}}\partial_{{\bf t}_{i+1}}\partial_{{\bf
t}_{i-1}}\Big{(}\sum_{i=0}^{2}\partial_{{\bf t}_{i}}u\Big{)}e_{h}d{\bf x}\\\
&=0,\end{aligned}$

where $\sum\int_{{\bf x}_{i}{\bf x}_{i+1}}(\partial_{{\bf
t}_{i}}\partial_{{\bf t}_{i-1}}^{2}u)e_{h}d{\bf s}=0$ as ${\bf t}_{i}$ are
parallel but in opposite directions on the two sides of the edge ${\bf
x}_{i}{\bf x}_{i+1}$, $\sum\int_{{\bf x}_{i+1}{\bf x}_{i-1}}(\partial_{{\bf
t}_{i}}\partial_{{\bf t}_{i-1}}^{2}u)e_{h}d{\bf s}=0$ as ${\bf t}_{i}$ lies on
the edge but in opposite directions on the two sides of the edge ${\bf
x}_{i}{\bf x}_{i+1}$, and $\sum_{i=0}^{2}\partial_{{\bf t}_{i}}u=0$ as
$\sum_{i=0}^{2}{{\bf t}_{i}}={\bf 0}$. By (3.8), (3.9) and (3.10), the theorem
is proved. ∎

######  Theorem 3.2.

Let $u\in H^{4,\infty}(\Omega)\cap H^{1}_{0}(\Omega)$ and $u_{h}\in V_{h}$ of
(2.2) be the solutions of Poisson’s equation (1.1) and the finite element
equation (2.3), respectively. It holds that

$\displaystyle\|I_{h}u-u_{h}\|_{L^{2}}\leq C\|I_{h}u-u_{h}\|_{L^{\infty}}\leq
Ch^{4}|u|_{H^{4,\infty}}.$

###### Proof.

Let $g^{j}_{h}\in V_{h}$ be the discrete Green function at an interior vertex
${\bf x}_{j}\in\mathcal{N}_{h}$ satisfying

$\displaystyle(\nabla g^{j}_{h},\nabla v_{h})=v_{h}({\bf x}_{j})\quad\forall
v_{h}\in{V}_{h}.$

Because of (3.5), we have, similar to (3.6),

$\displaystyle\begin{aligned} (I_{h}u-u_{h})({\bf
x}_{j})&=(\nabla(I_{h}u-u),\nabla g^{j}_{h})\\\
&=\sum_{T\in\mathcal{T}_{h}^{6}}(\nabla d_{h},\nabla
g^{j}_{h})_{T}+\sum_{T\in\overline{\mathcal{T}}_{h}}(\nabla(\overline{I}_{h}u-u),\nabla
g^{j}_{h})_{T}\\\ &=\sum_{T\in\mathcal{T}_{h}}\sum_{i=0}^{2}\int_{{\bf
x}_{i}{\bf x}_{i+1}}(I_{h}u-u)\partial_{{\bf t}_{i+2}^{\perp}}g^{j}_{h}d{\bf
s},\end{aligned}$

where $d_{h}$ is defined in (3.4), and other notations are defined in (3.6).

Repeating the proof from (3.8) to (3.10), i.e., replacing $e_{h}$ there by
$g^{j}_{h}$, we get

(3.11) $\displaystyle|(I_{h}u-u_{h})({\bf x}_{j})|$
$\displaystyle=\bigg{|}\sum_{T\in\overline{\mathcal{T}}_{h}}h^{4}\int_{\partial
T}(P_{4})(\partial_{{\bf t}}^{4}u)(\partial_{{\bf t}}g^{j}_{h})d{\bf
s}\bigg{|}$ (3.12) $\displaystyle\leq
Ch^{4}|u|_{H^{4,\infty}}\bigg{(}\sum_{T\in\mathcal{T}_{h}}|g^{j}_{h}|_{H^{1}(T)}^{2}\bigg{)}^{1/2}$
(3.13) $\displaystyle\leq
Ch^{4}|u|_{H^{4,\infty}}|\overline{g}_{h}^{j}|_{H^{1}},$

where ${\bf t}$ is a unit tangential vector on each edge of $T$, and
$\overline{g}_{h}^{j}$ is the discrete function in $\overline{V}_{h}$, i.e.,
for a ${\bf x}_{j}\in\overline{\mathcal{N}}_{h}$,

$\displaystyle(\nabla\overline{g}_{h}^{j},\nabla v_{h})=v_{h}({\bf
x}_{j})\quad\forall v_{h}\in\overline{V}_{h}.$

Therefore, for a ${\bf x}_{j}\in\mathcal{N}_{h}$, as
$V_{h}\subset\overline{V}_{h}$,

$\displaystyle(\nabla{g}_{h}^{j},\nabla v_{h})=v_{h}({\bf
x}_{j})=(\nabla\overline{g}_{h}^{j},\nabla v_{h})\quad\forall
v_{h}\in{V}_{h}.$

It implies that ${g}_{h}^{j}$ is an $H^{1}$-orthogonal projection of
$\overline{g}_{h}^{j}$ in $V_{h}$. By (3.11) and the $H^{1,\infty}$ stability
of $\overline{g}_{h}^{j}$ from [26], we have

$\displaystyle|I_{h}u-u_{h}|_{L^{\infty}}$ $\displaystyle=\max_{{\bf
x}_{j}\in\mathcal{N}_{h}}|(I_{h}u-u_{h})({\bf x}_{j})|$ $\displaystyle\leq
Ch^{4}|u|_{H^{4,\infty}}|{g}_{h}^{j}|_{H^{1}}$ $\displaystyle\leq
Ch^{4}|u|_{H^{4,\infty}}|\overline{g}_{h}^{j}|_{H^{1}}$ $\displaystyle\leq
Ch^{4}|u|_{H^{4,\infty}}|\overline{g}_{h}^{j}|_{H^{1,\infty}}$
$\displaystyle\leq Ch^{4}|u|_{H^{4,\infty}}.$

Consequently,

$\displaystyle\|I_{h}u-u_{h}\|_{L^{2}}\leq\hbox{\rm{meas}}(\Omega)\|I_{h}u-u_{h}\|_{\infty}\leq
Ch^{4}|u|_{H^{4,\infty}}.$

The proof is complete. ∎

##  4\. Lift the $P_{1}$ solution to a $P_{3}$ solution

We can choose a larger patch, and/or several patterns of patches, for locally
lifting a solution. For a simple way of computation, we choose only one patch
$R$ of shape of one regular triangle of edge-size $2h$ on one honeycomb mesh
$\mathcal{T}_{h}$, cf. Figure 4.

$R:$$\mathcal{R}_{h}:$$e_{1}$$e_{2}$$e_{3}$ Figure 4. Left: A recovery patch
$R$ of an equilateral triangle; Right: a recovery patch grid
$\mathcal{R}_{h}=\\{\hbox{\rm{regular triangle }}\ R\\}$ on a regular hexagon
domain.

On each patch $R$, there are 11 vertices $\\{{\bf x}_{j},\
j=1,\dots,11\\}\subset\mathcal{N}_{h}$ of $\mathcal{T}_{h}$, cf. Figure 4. Let
$\tilde{u}_{h}\in P_{3}(R)$ be the least squares solution of following
$11\times 10$ linear system of equations,

(4.1) $\displaystyle\tilde{u}({\bf x}_{j})=u_{h}({\bf x}_{j}),\quad\
j=1,\dots,11,$

where $u_{h}\in V_{h}$ is the $P_{1}$ solution of (2.3).

######  Lemma 4.1.

The linear system (4.1) of equations has a unique solution.

###### Proof.

The coefficient matrix of (4.1) is of size $11\times 10$, as $\dim
P_{3}(R)=10$. To show it is a full-rank matrix, we can show $\tilde{u}_{h}=0$
if

$\displaystyle\tilde{u}_{h}({\bf x}_{j})=0,\quad\ j=1,\dots,11.$

Because the $P_{3}$ polynomial $\tilde{u}_{h}$ vanishes at four points on edge
$e_{1}$, cf. Figure 4, it vanishes on the whole edge $e_{1}$. Thus

$\displaystyle\tilde{u}_{h}=l_{1}p_{2},$

where $p_{2}\in P_{2}(R)$ and $l_{1}\in P_{1}(R)$ which vanishes only on
$e_{1}$.

Because $\tilde{u}_{h}$ vanishes but $l_{1}$ does not vanish at three top
points of edge $e_{2}$, cf. Figure 4, $p_{2}$ vanishes on the whole edge
$e_{2}$. Thus

$\displaystyle\tilde{u}_{h}=l_{1}l_{2}p_{1},$

where $p_{1}\in P_{1}(R)$ and $l_{2}\in P_{1}(R)$ which vanishes only on
$e_{2}$.

Because $\tilde{u}_{h}$ vanishes but $l_{1}$ and $l_{2}$ do not vanish at two
internal points of edge $e_{3}$, cf. Figure 4, $p_{1}$ vanishes on the whole
edge $e_{3}$. Thus

$\displaystyle\tilde{u}_{h}=Cl_{1}l_{2}l_{3},$

where $C\in P_{0}(R)$ and $l_{3}\in P_{1}(R)$ which vanishes only on $e_{3}$.

There are two internal ${\bf x}_{j}$ inside $R$. Because $\tilde{u}_{h}({\bf
x}_{j})=0$, $l_{1}({\bf x}_{j})\neq 0$, $l_{2}({\bf x}_{j})\neq 0$ and
$l_{3}({\bf x}_{j})\neq 0$, cf. Figure 4, $C=0$, i.e. $\tilde{u}_{h}=0$.

Thus the matrix in (4.1) is of full rank and the solution $\tilde{h}$ is well
defined uniquely by the least squares method. ∎

######  Theorem 4.2.

Let $u\in H^{1}_{0}(\Omega)\cap H^{4,\infty}(\Omega)$ and $\tilde{u}_{h}$ be
the solutions of Poisson’s equation (1.1) and the lifting equation (4.1),
respectively. It holds that

$\displaystyle\|u-\tilde{u}_{h}\|_{L^{2}}+h|u-\tilde{u}_{h}|_{H^{1},h}\leq
Ch^{4}|u|_{H^{4,\infty}},$

where $|u-\tilde{u}_{h}|_{H^{1},h}$ is defined piecewise on $\mathcal{R}_{h}$.

###### Proof.

Let $\tilde{I}_{h}u=\tilde{u}_{h}\in P_{3}(R)$ be the least squares solution
of (4.1). Because the interpolation $\tilde{I}_{h}u$ preserves $P_{3}$
polynomials and is stable (cf. [28]), we have

$\displaystyle\|u-\tilde{I}_{h}u\|_{L^{2}}+h|u-\tilde{I}_{h}u|_{H^{1},h}\leq
Ch^{4}|u|_{H^{4}}.$

As $\tilde{u}_{h}\in P_{3}(R)$, $\tilde{I}_{h}\tilde{u}_{h}=\tilde{u}_{h}$. By
definition (4.1), $\tilde{I}_{h}u=\tilde{I}_{h}I_{h}u$. Thus, as
$\tilde{I}_{h}$ is stable (cf. [28]),

$\displaystyle\quad\ \|u-\tilde{u}_{h}\|_{L^{2}}+h|u-\tilde{u}_{h}|_{H^{1},h}$
$\displaystyle\leq\|u-\tilde{I}_{h}u\|_{L^{2}}+\|\tilde{I}_{h}(I_{h}u-u_{h})\|_{L^{2}}$
$\displaystyle\qquad+h|u-\tilde{u}_{h}|_{H^{1},h}+h|\tilde{I}_{h}(I_{h}u-u_{h})\|_{H^{1},h}$
$\displaystyle\leq
Ch^{4}|u|_{H^{4,\infty}}+C\|I_{h}u-u_{h}\|_{L^{2}}+Ch^{4}|u|_{H^{4,\infty}}+Ch|I_{h}u-u_{h}|_{H^{1},h}$
$\displaystyle\leq Ch^{4}|u|_{H^{4,\infty}}.$

The theorem is proved. ∎

##  5\. Numerical test

We solve the Poisson equation (1.1), on a regular hexagon domain $\Omega$ of
size $1$ centered at the origin. The exact solution of (1.1) for the numerical
test is

(5.1) $\displaystyle\begin{aligned}
u(x,y)&=x^{2}\sin(\frac{\pi}{2}(\frac{y}{\sqrt{3}}+x+1))\sin\frac{\pi}{2}(\frac{y}{\sqrt{3}}-x+1))\\\
&\qquad\cdot\sin(\frac{\pi}{\sqrt{3}}(y+\frac{\sqrt{3}}{2})).\end{aligned}$

The first three grids are depicted in Figure 1. Here the refinement can not be
nested. Some honeycombs are split and combined with some neighbor honeycombs’
debris to form half-sized honeycombs. The higher level grids are obtained by
such half-sized refinements.

In Table 1, we list the error between the interpolation and the $P_{1}$
virtual element solution to show two-order or three-order superconvergence.
The optimal order of convergence is $2$ in $L^{2}$ and $L^{\infty}$ norms, but
we have an order $4$ for them in Table 1, i.e., order-two superconvergence.
The optimal order of convergence is $1$ for $P_{1}$ elements in $H^{1}$-norm,
but we have an order $4$ for it in Table 1, i.e., order-three
superconvergence. The results match the theory on the superconvergence. In
Table 1, on first level, there is no internal degree of freedom in $V_{h}$
that the errors are supposedly zero. They are round-off errors when evaluating
the zero boundary condition.

Table 1. Error profile of $u_{h}$ on Figure 1 meshes for (5.1).  Grid | $\|I_{h}u-u_{h}\|_{L^{2}}$ | $O(h^{r})$ | $|I_{h}u-u_{h}|_{H^{1}}$ | $O(h^{r})$ | $\|I_{h}u-u_{h}\|_{L^{\infty}}$ | $O(h^{r})$
---|---|---|---|---|---|---
1 | 0.1233E-31 | 0.00 | 0.2013E-31 | 0.00 | 0.1735E-31 | 0.00
2 | 0.1567E-02 | 0.00 | 0.4800E-02 | 0.00 | 0.1829E-02 | 0.00
3 | 0.1138E-03 | 3.78 | 0.3305E-03 | 3.86 | 0.1480E-03 | 3.63
4 | 0.7358E-05 | 3.95 | 0.2118E-04 | 3.96 | 0.9840E-05 | 3.91
5 | 0.4641E-06 | 3.99 | 0.1332E-05 | 3.99 | 0.6244E-06 | 3.98
6 | 0.2907E-07 | 4.00 | 0.8337E-07 | 4.00 | 0.3917E-07 | 3.99
7 | 0.1818E-08 | 4.00 | 0.5213E-08 | 4.00 | 0.2450E-08 | 4.00
8 | 0.1137E-09 | 4.00 | 0.3258E-09 | 4.00 | 0.1532E-09 | 4.00
9 | 0.7104E-11 | 4.00 | 0.2037E-10 | 4.00 | 0.9575E-11 | 4.00
10 | 0.4440E-12 | 4.00 | 0.1279E-11 | 3.99 | 0.5985E-12 | 4.00

In Table 2, we list the $L^{2}$-error between the solution and the $P_{1}$
virtual element solution, showing it can only converge at the optimal order in
$L^{2}$-norm, order 2. We also list the $L^{2}$ and the $H^{1}$ errors between
the exact solution and the $P_{3}$ post-processed solution. The results match
the theory on the optimal-order convergence of the lifted solution. In Table
2, we start to do lifting from level three.

Table 2. Error profile of $u_{h}$ and $\tilde{u}_{h}$ on Figure 1 meshes for (5.1).  Grid | $\|u-u_{h}\|_{L^{2}}$ | $O(h^{r})$ | $\|u-\tilde{u}_{h}\|_{L^{2}}$ | $O(h^{r})$ | $|u-\tilde{u}_{h}|_{H^{1},h}$ | $O(h^{r})$
---|---|---|---|---|---|---
1 | 0.8686E-01 | 0.00 | 0.0000E+00 | 0.00 | 0.0000E+00 | 0.00
2 | 0.3127E-01 | 1.47 | 0.0000E+00 | 0.00 | 0.0000E+00 | 0.00
3 | 0.9417E-02 | 1.73 | 0.3298E-02 | 0.00 | 0.2263E-01 | 0.00
4 | 0.2448E-02 | 1.94 | 0.2324E-03 | 3.83 | 0.2648E-02 | 3.10
5 | 0.6179E-03 | 1.99 | 0.1464E-04 | 3.99 | 0.3256E-03 | 3.02
6 | 0.1548E-03 | 2.00 | 0.9192E-06 | 3.99 | 0.4060E-04 | 3.00
7 | 0.3873E-04 | 2.00 | 0.5752E-07 | 4.00 | 0.5072E-05 | 3.00
8 | 0.9684E-05 | 2.00 | 0.3596E-08 | 4.00 | 0.6339E-06 | 3.00
9 | 0.2421E-05 | 2.00 | 0.2248E-09 | 4.00 | 0.7924E-07 | 3.00
10 | 0.6053E-06 | 2.00 | 0.1405E-10 | 4.00 | 0.9905E-08 | 3.00

##  6\. Ethical Statement

###  6.1. Compliance with Ethical Standards

The submitted work is original and is not published elsewhere in any form or
language.

###  6.2. Funding

Yanping Lin is supported in part by HKSAR GRF 15302922 and polyu-CAS joint
Lab.

Xuejun Xu is supported by National Natural Science Foundation of China (Grant
Nos. 12071350), Shanghai Municipal Science and Technology Major Project No.
2021SHZDZX0100, and Science and Technology Commission of Shanghai
Municipality.

###  6.3. Conflict of Interest

There is no potential conflict of interest.

###  6.4. Ethical approval

This article does not contain any studies involving animals. This article does
not contain any studies involving human participants.

###  6.5. Informed consent

This research does not have any human participant.

###  6.6. Availability of supporting data

This research does not use any external or author-collected data.

###  6.7. Authors’ contributions

All authors made equal contribution.

###  6.8. Acknowledgments

None.

## References

  * [1] A. Al-Twaeel, X. Wang, X. Ye and S. Zhang, A stabilizer free weak Galerkin element method with supercloseness of order two, Numer. Methods Partial Differential Equations 37 (2021), no. 2, 1012–1029.
  * [2] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini and A. Russo, Basic principles of virtual element methods, Math. Models Methods Appl. Sci. 23 (2013) 199–214.
  * [3] L. Beirão da Veiga, F. Brezzi, L. D. Marini and A. Russo, H(div) and H(curl)-conforming virtual element methods, Numer. Math. 133 (2016) 303–332.
  * [4] S. Berrone, A. Borio, F. Marcon and G. Teora, A first-order stabilization-free virtual element method, Appl. Math. Lett. 142 (2023), Paper No. 108641, 6 pp.
  * [5] S. Berrone, A. Borio and F. Marcon, Comparison of standard and stabilization free Virtual Elements on anisotropic elliptic problems, 2022, arXiv:2202.08571v1.
  * [6] S. Berrone, A. Borio and F. Marcon, Lowest order stabilization free Virtual Element Method for the Poisson equation, 2023, arXiv:2103.16896v1.
  * [7] H. Blum, Q. Lin and R. Rannacher, Asymptotic error expansion and Richardson extrapolation for linear finite elements, Numer. Math. 49 (1986), no. 1, 11–37.
  * [8] J. H. Bramble and A. H. Schatz, Higher order local accuracy by averaging in the finite element method, Math. Comp. 31 (1977), pp. 94–111.
  * [9] S. Cao, L. Chen and R. Guo, A virtual finite element method for two-dimensional Maxwell interface problems with a background unfitted mesh, Math. Models Methods Appl. Sci. 31 (2021), no. 14, 2907–2936.
  * [10] S. Cao, L. Chen and R. Guo, Immersed virtual element methods for electromagnetic interface problems in three dimensions, Math. Models Methods Appl. Sci. 33 (2023), no. 3, 455–503.
  * [11] H. Chen and J. Wang, An interior estimate of superconvergence for finite element solutions for second-order elliptic problems on quasi-uniform meshes by local projections, SIAM J. Numer.Anal. 41 (2003), no. 4, 1318–1338.
  * [12] L. Chen, H. Wei and M. Wen, An interface-fitted mesh generator and virtual element methods for elliptic interface problems, J. Comput. Phys. 334 (2017), 327–348.
  * [13] L. Chen and J. Huang, Some error analysis on virtual element methods, Calcolo 55 (2018), no. 1, Paper No. 5, 23 pp.
  * [14] L. Chen and X. Huang, Nonconforming virtual element method for 2mth order partial differential equations in Rn, Math. Comp. 89 (2020), no. 324, 1711–1744.
  * [15] J. Douglas and T. Dupont, Superconvergence for Galerkin methods for the two-point boundary problem via local projections, Numer. Math. 21 (1973), 270–278.
  * [16] F. Feng, W. Han and J. Huang, Virtual element methods for elliptic variational inequalities of the second kind, J. Sci. Comput. 80 (2019), no. 1, 60–80.
  * [17] F. Feng, J. Huang and Y. Yu, A non-consistent virtual element method for reaction diffusion equations, East Asian J. Appl. Math. 10 (2020), no. 4, 786–799.
  * [18] J. Huang and S. Lin, A C0P2 time-stepping virtual element method for linear wave equations on polygonal meshes, Electron. Res. Arch. 28 (2020), no. 2, 911–933.
  * [19] J. Huang and Y. Yu, Some estimates for virtual element methods in three dimensions, Comput. Methods Appl. Math. 23 (2023), no. 1, 177–187.
  * [20] J. Huang and Y. Yu, A medius error analysis for nonconforming virtual element methods for Poisson and biharmonic equations, J. Comput. Appl. Math. 386 (2021), Paper No. 113229, 20 pp.
  * [21] Y. Huang and S. Zhang, Supercloseness of the divergence-free finite element solutions on rectangular grids, Communications in Mathematics and Statistics 1 (2013), 143–162.
  * [22] V. Kač and P. Cheung, Quantum Calculus, Springer, New York, 2001.
  * [23] M. Krizek and P. Neittaanmaki, On superconvergence techniques, Acta Appl. Math. 9 (1987), 175–198.
  * [24] Y. Lin, M. Mu and S. Zhang, Stabilizer-free polygonal and polyhedral virtual elements, 2023, [archivenumber].
  * [25] L. Mu, X. Ye and S. Zhang, A stabilizer free, pressure robust, and superconvergence weak Galerkin finite element method for the Stokes Equations on polytopal mesh, SIAM J. Sci. Comput. 43 (2021), no. 4, A2614–A2637.
  * [26] R. Rannacher and R. Scott, Some optimal error estimates for piecewise linear finite element approximations, Math. Comp. 38 (1982), no. 158, 437–445.
  * [27] A. H. Schatz, I. H. Sloan and L. B. Wahlbin, Superconvergence in finite element methods and meshes that are symmetric with respect to a point, SIAM J. Numer. Anal. 33 (1996), 505–521.
  * [28] L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp. 54 (1990), no. 190, 483–493.
  * [29] L. B. Wahlbin, Superconvergence in Galerkin Finite Element Methods, Lecture Notes in Math. 1605, Springer-Verlag, New York, 1995.
  * [30] C. Wang, Superconvergence of Ritz-Galerkin finite element approximations for second order elliptic problems, Numer. Methods Partial Differential Equations 34 (2018), no. 3, 838–856.
  * [31] G. Wang, L. Mu, Y. Wang and Y. He, A pressure-robust virtual element method for the Stokes problem, Comput. Methods Appl. Mech. Engrg. 382 (2021), Paper No. 113879, 19 pp.
  * [32] J. Wang, A superconvergence analysis for finite element solutions by the least-squares surface fitting on irregular meshes for smooth problems, J. Math. Study, 33 (2000), 229–243.
  * [33] J. Wang, X. Wang, X. Ye, S. Zhang and Peng Zhu, Two-order superconvergence for a weak Galerkin method on rectangular and cuboid grids, Numer. Methods Partial Differential Equations 39 (2023), no. 1, 744–758.
  * [34] J. Wang, X. Wang, X. Ye, S. Zhang and Peng Zhu, On the superconvergence of a WG method for elliptic problem with variable coefficients, Science China, 2023, doi:10.1007/s11425-022-2097-8.
  * [35] Y. Wang, G. Wang and F. Wang, An adaptive virtual element method for incompressible flow, Comput. Math. Appl. 101 (2021), 63–73.
  * [36] X. Xu and S. Zhang, A family of stabilizer-free virtual elements on triangular meshes, 2023, arXiv:2309.09660.
  * [37] X. Ye and S. Zhang, A stabilizer free weak Galerkin finite element method on polytopal mesh: Part II, J. Comput. Appl. Math. 394 (2021), Paper No. 113525, 11 pp.
  * [38] X. Ye and S. Zhang, A stabilizer free weak Galerkin finite element method on polytopal mesh: Part III, J. Comput. Appl. Math. 394 (2021), Paper No. 113538, 9 pp.
  * [39] X. Ye and S. Zhang, A stabilizer free WG Method for the Stokes Equations with order two superconvergence on polytopal mesh, Electron. Res. Arch. 29 (2021), no. 6, 3609–3627.
  * [40] X. Ye and S. Zhang, A $P_{k+2}$ polynomial lifting operator on polygons and polyhedrons, Appl. Math. Lett. 116 (2021), Paper No. 107033, 6 pp.
  * [41] X. Ye and S. Zhang, Achieving superconvergence by one-dimensional discontinuous finite elements: the CDG method, East Asian J. Appl. Math. 12 (2022), no. 4, 781–790.
  * [42] X. Ye and S. Zhang, Order two superconvergence of the CDG method for the Stokes equations on triangle/tetrahedron, Journal of Applied Analysis and Computation, 12 (2022), no. 6, 2578–2592.
  * [43] X. Ye and S. Zhang, Achieving superconvergence by one-dimensional discontinuous finite elements: weak Galerkin method, East Asian J. Appl. Math. 12 (2022), no. 3, 590–598.
  * [44] X. Ye and S. Zhang, Constructing order two superconvergent WG finite elements on rectangular meshes, Numer. Math. Theory Methods Appl. 16 (2023), no. 1, 230–241.
  * [45] X. Ye and S. Zhang, Order two superconvergence of the CDG finite elements on triangular and tetrahedral meshes, CSIAM Trans. Appl. Math. 4 (2023), no. 2, 256–274.
  * [46] X. Ye and S. Zhang, Four-order superconvergent weak Galerkin methods for the biharmonic equation on triangular meshes, Communications on Applied Mathematics and Computation, Doi:10.1007/s42967-022-00201-5, 2023.
  * [47] X. Ye and S. Zhang, Four-order superconvergent CDG finite elements for the biharmonic equation on triangular meshes, J. Comput. Appl. Math., 2023, doi: 10.1016/j.cam.2023.115516.
  * [48] M. Zlamal, Superconvergence and reduced integration in the finite element method, Math Comp. 32 (1978), 663–685.
  * [49] Q. Zhu and Q. Lin, Superconvergence Theory of the Finite Element Methods, Hunan Science Press, Changsha, China, 1989.
  * [50] O. C. Zienkiewicz and J. Z. Zhu, The superconvergence patch recovery and a posteriori error estimates, Part 2, Error estimates and adaptivity, Internat. J. Numer. Methods Engrg. 33 (1992), 1365–1382.