# Integer partitions detect the primes

William Craig, Jan-Willem van Ittersum and Ken Ono  Department of Mathematics
and Computer Science, University of Cologne, Weyertal 86-90, 50931 Cologne,
Germany<EMAIL_ADDRESS><EMAIL_ADDRESS>Dept. of Mathematics,
University of Virginia, Charlottesville, VA 22904<EMAIL_ADDRESS>
###### Abstract.

We show that integer partitions, the fundamental building blocks in additive
number theory, detect prime numbers in an unexpected way. Answering a question
of Schneider, we show that the primes are the solutions to special equations
in partition functions. For example, an integer $n\geq 2$ is prime if and only
if

$(3n^{3}-13n^{2}+18n-8)M_{1}(n)+(12n^{2}-120n+212)M_{2}(n)-960M_{3}(n)=0,$

where the $M_{a}(n)$ are MacMahon’s well-studied partition functions. More
generally, for MacMahonesque partition functions $M_{\vec{a}}(n),$ we prove
that there are infinitely many such prime detecting equations with constant
coefficients, such as

$80M_{(1,1,1)}(n)-12M_{(2,0,1)}(n)+12M_{(2,1,0)}(n)+\dots-12M_{(1,3)}(n)-39M_{(3,1)}(n)=0.$

###### Key words and phrases:

partitions, primes, quasimodular forms

###### 2020 Mathematics Subject Classification:

Primary 11P81, 05A17; Secondary 11Fxx

##  1\. Introduction and Statement of Results

Diophantine sets $S$ are subsets of $\mathbb{Z}_{+}^{n}$ that arise from the
positive integer solutions of a Diophantine equation. For example, the set of
composite numbers is Diophantine, which is easily seen by considering the
polynomial equation

$C(x_{1},x_{2},x_{3}):=x_{1}-(x_{2}+1)(x_{3}+1)=0.$

Indeed, we have that $S_{C}:=\\{x_{1}\in\mathbb{Z}_{+}\ :\
C(x_{1},x_{2},x_{3})=0\ {\text{\rm with}}\ x_{2},x_{3}\in\mathbb{Z}_{+}\\}$ is
the set of composite numbers. In 1970, Matiyasevich [1] famously proved that a
set is Diophantine if and only if it is computably enumerable. Building on
work by Davis, Putnam, and Robinson, this result resolved Hilbert’s tenth
problem in the negative: There is no algorithm that determines whether a
generic Diophantine equation has integer solutions.

In contrast to the simple case of the composite numbers, Matiyasevich’s work
implies that the set of primes is Diophantine, and even more, is polynomially
representable. In 1976, Jones, Sato, Wada, and Wiens [2] made this explicit;
they produced a degree 25 polynomial in 26 variables whose positive values, as
the variables vary over nonnegative integers, is the set of primes.

In analogy, we consider whether the set of primes is partition theoretic, by
which we mean solutions to equations in partition functions (for background on
partitions see [3]). This is a question of Schneider, as an example of his
general philosophy and research program where classical number theoretic
topics (e.g. multiplicative number theory, prime number distribution, and
zeta-functions) are informed by partition theory (for example, see [4, 5] and
his works on arithmetic densities [6] and partition zeta-functions [7]).

We show that the primes are partition theoretic in infinitely many ways. Our
first glimpse of this phenomenon arises from MacMahon’s [8] $q$-series:

(1.1) $\displaystyle\mathcal{U}_{a}(q)=\sum_{n\geq
1}M_{a}(n)\,q^{n}:=\sum_{0<s_{1}<s_{2}<\cdots<s_{a}}\frac{q^{s_{1}+s_{2}+\cdots+s_{a}}}{(1-q^{s_{1}})^{2}(1-q^{s_{2}})^{2}\cdots(1-q^{s_{a}})^{2}}.$

For $a\geq 1$, the partition function $M_{a}(n)$ sums the products of the
_multiplicities_ of partitions of $n$ with $a$ different part sizes. To be
precise, we have

(1.2) $M_{a}(n)=\sum_{\begin{subarray}{c}0<s_{1}<s_{2}<\dots<s_{a}\\\
n=m_{1}s_{1}+m_{2}s_{2}+\dots+m_{a}s_{a}\end{subarray}}m_{1}m_{2}\cdots
m_{a}.$

We find linear equations in these functions, with polynomial coefficients,
that detect the primes.

######  Theorem 1.1.

The following are true:

(1) For positive integers $n,$ we have

$(n^{2}-3n+2)M_{1}(n)-8M_{2}(n)\geq 0,$

and for $n\geq 2,$ this expression vanishes if and only if $n$ is prime.

(2) For positive integers $n,$ we have

$(3n^{3}-13n^{2}+18n-8)M_{1}(n)+(12n^{2}-120n+212)M_{2}(n)-960M_{3}(n)\geq 0,$

and for $n\geq 2,$ this expression vanishes if and only if $n$ is prime.

######  Remark.

The expressions in Theorem 1.1 vanish trivially when $n=1.$

######  Remark.

By taking suitable linear combinations (with polynomial coefficients) of the
two expressions in Theorem 1.1, it is simple to produce infinitely many
further examples. For example, we have

$(3n^{4}-3n^{3}-17n^{2}+27n-10)M_{1}(n)-(240n-400)M_{2}(n)-1920M_{3}(n).$

There are more examples of such prime detecting expressions involving the
$M_{a}(n)$ functions.

######  Theorem 1.2.

Suppose that $E(n)$ is in Table 1 in the Appendix. For positive integers $n,$
we have $E(n)\geq 0,$ and for $n\geq 2,$ we have $E(n)=0$ if and only if $n$
is prime.

Based on a computer search for further such prime detecting partition
equations, we believe that this list generates all such expressions involving
MacMahon’s $M_{a}(n)$ partition functions.

######  Conjecture.

Let $\vec{P}(x):=(p_{1}(x),p_{2}(x),\dots,p_{a}(x))\in\mathbb{Z}[x]^{a}$ be a
vector of relatively prime integer polynomials. For integers $n\geq 2,$
suppose that

$E(n)=p_{1}(n)M_{1}(n)+p_{2}(n)M_{2}(n)+\dots+p_{a}(n)M_{a}(n)\geq 0,$

and vanishes precisely on the primes. Then $E(n)$ is a $\mathbb{Q}[n]$-linear
combination of entries in Table 1.

Despite this conjecture, we show that Theorems 1.1 and 1.2 are glimpses of a
rich infinite family of such expressions. To make this precise, we define
natural generalizations of MacMahon’s $M_{a}(n).$ For $a\geq 1$ and a vector
$\vec{a}=(v_{1},v_{2},\dots,v_{a})\in\mathbb{N}^{a}$, we define the
MacMahonesque partition function111These functions were also considered in
Bachmann’s master’s thesis (see Chapter 5 of [9]).

(1.3) $M_{\vec{a}}(n):=\sum_{\begin{subarray}{c}0<s_{1}<s_{2}<\dots<s_{a}\\\
m_{1},\ldots,m_{a}>0\\\
n=m_{1}s_{1}+m_{2}s_{2}+\dots+m_{a}s_{a}\end{subarray}}m_{1}^{v_{1}}m_{2}^{v_{2}}\cdots
m_{a}^{v_{a}},$

certain sums of degree $|\vec{a}|:=\nu_{1}+\cdots+\nu_{a}$ monomials in the
part multiplicities.

######  Remark.

We note that $M_{a}(n)\neq M_{(a)}(n),$ when $a>1.$ Instead, we have
$M_{a}(n)=M_{\vec{a}}(n)$, where $\vec{a}=(1,\dots,1)\in\mathbb{N}^{a}.$

For these generalized MacMahon functions, we establish the existence of
infinitely many linear equations with integer coefficients whose set of
solutions, where $n\geq 2$, is the set of primes.

######  Theorem 1.3.

If $d\geq 4$, then the following are true.

(1) There are relatively prime integers $c_{\vec{a}}$ such that for positive
integers $n$ we have

$\sum_{|\vec{a}|\leq d}c_{\vec{a}}\,M_{\vec{a}}(n)\geq 0,$

and for $n\geq 2,$ this expression vanishes if and only if $n$ is prime.

(2) As $d\rightarrow\infty$, there are $\gg d^{2}$ many linearly independent
expressions satisying (1).

######  Examples.

We offer two examples of Theorem 1.3. For positive integers $n$, we find that

$\begin{split}\Psi_{1}(n):=63M_{(2,2)}(n)-12&M_{(3,0)}(n)-39M_{(3,1)}(n)-12M_{(1,3)}(n)\\\
&+80M_{(1,1,1)}(n)-12M_{(2,0,1)}(n)+12M_{(2,1,0)}(n)+12M_{(3,0,0)}(n)\geq
0,\end{split}$

and

$\begin{split}\Psi_{2}(n):=14M_{(1)}&(n)-15M_{(2)}(n)-2M_{(3)}(n)+3M_{(4)}(n)+30M_{(2,0)}(n)-72M_{(3,0)}(n)\\\
&-36M_{(2,1)}(n)-6M_{(4,0)}(n)-12M_{(3,1)}(n)+72M_{(2,1,0)}(n)+72M_{(3,0,0)}(n)\geq
0.\end{split}$

For $n\geq 2$ these expressions vanish if and only if $n$ is prime. Finally,
we note the nontrivial identity $\Psi_{2}(n)=\frac{36}{11}\Psi_{1}(n).$

In addition to the equations covered by Theorem 1.3, there are many prime
detecting equations involving convolutions of the MacMahonesque functions. For
example, one can show that

$\begin{split}\Psi_{3}(n):=10M_{(1)}&(n)-17M_{(3)}(n)+7M_{(5)}(n)\\\
+&12\sum_{i+j=n}M_{(1)}(i)\left(M_{(1)}(j)-10M_{(3)}(j)\right)+96\\!\\!\sum_{i+j+\ell=n}\\!\\!M_{(1)}(i)M_{(1)}(j)M_{(1)}(\ell)\geq
0,\end{split}$

and similarly detects primes. However, it turns out that this expression can
be reformulated much more simply thanks to the identity
$\Psi_{3}(n)=\Psi_{2}(n).$

Such simplifications always exist. It turns out that arbitrary convolutions of
MacMahonesque functions are always linear combinations of such functions with
constant coefficients. Furthermore, linear combinations of MacMahonesque
functions with polynomial coefficients can be reduced similarly. To be
precise, we have the following result which is implicit in the proof of
Theorem 1.3.

######  Theorem 1.4.

The following are true.

(1) If $\vec{\alpha}\in\mathbb{Z}_{+}^{a}$ and
$\vec{\beta}\in\mathbb{Z}_{+}^{b}$, then there are rational numbers
$c_{\vec{a}}$ such that for every $n$ we have

$\sum_{i+j=n}M_{\vec{\alpha}}(i)M_{\vec{\beta}}(j)=\sum_{|\vec{a}|\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}|\vec{\alpha}|+|\vec{\beta}|+a+b-1}}c_{\vec{a}}M_{\vec{a}}(n).$

(2) If $\vec{\alpha}\in\mathbb{Z}_{+}^{a},$ then there are rational numbers
$c_{\vec{a}}$ such that for every $n$ we have

$nM_{\vec{\alpha}}(n)=\sum_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}|\vec{a}|\leq|\vec{\alpha}|+a+1}}c_{\vec{a}}M_{\vec{a}}(n).$

######  Remarks.

(1) By Theorem 1.4 (1) and (2), every convolution product of MacMahonesque
functions with polynomial coefficients is a linear combination of such
functions with rational coefficients.

(2) As the proof will show, Theorem 1.4 is a partition theoretic reformulation
of beautiful results of Bachmann and Kühn (see Theorem 1.3 and 1.7 of [10]).

######  Examples.

The identity $\Psi_{3}(n)=\Psi_{2}(n)$ offers examples of linear expressions
with constant coefficients obtained from Theorem 1.4 (1). As a simpler example
of Theorem 1.4 (1), we have

$\sum_{i+j=n}M_{(1)}(i)\,M_{(1)}(j)={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{1}{6}M_{(3)}(n)+2M_{(1,1)}(n)-\frac{1}{6}M_{(1)}(n).}$

As an example of Theorem 1.4 (2), we have

$\begin{split}nM_{(1,1)}(n)=\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{1}{22}\Bigl{(}-21M_{(3,1)}(n)&+72M_{(2,2)}(n)-9M_{(1,3)}(n)+24M_{(3,0)}(n)\\\
&\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-24M_{(3,0,0)}(n)-24M_{(2,1,0)}(n)+24M_{(2,0,1)}(n)-72M_{(1,1,1)}(n)\Bigr{)}.\end{split}$

Finally, by iterating Theorem 1.4 (2), we can reformulate the expression in
Theorem 1.1 (1) as

$\begin{split}(n^{2}-3n+2)M_{1}(n)-&8M_{2}(n)={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{3}{11}\Psi_{1}(n)},\end{split}$

a linear expression in MacMahonesque functions with constant coefficients.

To obtain these results, we make use of the theory of multiple $q$-zeta
values, largely as developed by Bachmann and Kühn [9, 10, 11], and the theory
of quasimodular forms. In Section 2, we recall facts about quasimodular forms,
and we establish conditions (see Theorem 2.3) under which their Fourier
expansions detect primes, a result which is of independent interest. In
Section 3, we use the fact that MacMahon’s $q$-series $\mathcal{U}_{a}(q)$
arise naturally as quasimodular forms (see [10, 12, 13]), combined with
Theorem 2.3, to prove Theorems 1.1 and 1.2.

In contrast, the generating series of the generalized MacMahonesque partition
functions

$\mathcal{U}_{\vec{a}}(q):=\sum_{n\geq 1}M_{\vec{a}}(n)q^{n}$

are generally not quasimodular. However, they occur naturally in the theory of
multiple $q$-zeta values. Bachmann and Kühn show [10] that the MacMahonesque
series generate an extraordinary differential algebra whose arithmetic
properties encode a rich web of identities satisfied by the MacMahonesque
generating functions. To make use of quasimodularity, for those $\vec{a}$ with
odd entries, we construct quasimodular forms
$\mathcal{U}_{\vec{a}}^{\operatorname{sym}}(q)$ by taking traces over the
natural action of the symmetric group $S_{a}$ on $\mathcal{U}_{\vec{a}}(q).$
In particular, we prove the following theorem (see Theorem 4.4 (2)).

######  Theorem 1.5.

All quasimodular forms are linear combinations of the MacMahonesque
$\mathcal{U}_{\vec{a}}^{\operatorname{sym}}(q)$.

Theorem 1.3 then follows by applying Theorem 2.3. These details are carried
out in Section 4.

## Acknowledgements

The authors would like to thank the organizers of the International Conference
on Modular Forms and $q$-Series at the University of Cologne, where this
project emerged. The authors would also like to thank Tewodros Amdeberhan,
George Andrews, Henrik Bachmann, Robert Schneider, and Ajit Singh for helpful
comments which improved the manuscript. The first and second authors thank the
support of the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (grant agreement No. 101001179)
and by the SFB/TRR 191 “Symplectic Structure in Geometry, Algebra and
Dynamics”, funded by the DFG (Projektnummer 281071066 TRR 191). The third
author thanks the Thomas Jefferson Fund and grants from the NSF (DMS-2002265
and DMS-2055118).

##  2\. Quasimodular forms that detect primes

Central to all of our results is a theorem (see Theorem 2.3) that
characterizes quasimodular forms in the space consisting of linear
combinations of derivatives of Eisenstein series whose coefficients are zero
at primes and positive at composite $n\geq 2$.

###  2.1. Preliminaries about quasimodular forms

To state and then prove this result, we require several important preliminary
facts. For $\tau$ in the upper-half of the complex plane, set $q:=e^{2\pi
i\tau}$. For integers $k\geq 1,$ we consider the Bernoulli numbers $B_{k}$ and
the Eisenstein series

(2.1) $G_{k}(\tau):=-\dfrac{B_{k}}{2k}+\sum_{n\geq 1}\sigma_{k-1}(n)q^{n},$

where $\sigma_{k-1}(n):=\sum_{d|n}d^{k-1}.$ It is well-known (for example, see
[14]) that $G_{k}$ is quasimodular if and only if $k$ is even, and modular if
$k\geq 4$ and even. Furthermore, the algebra of quasimodular forms, which we
denote by $\widetilde{M}$, is generated by $G_{2},G_{4}$, and $G_{6}$. We also
denote by $\widetilde{M}_{k}$ the subspace of $\widetilde{M}$ of forms with
weight $k$, and $\widetilde{M}_{\leq k}$ the space of all quasimodular forms
of mixed weight $\leq k$.

We require the important $q$-differential operator

(2.2) $\displaystyle D:=\frac{1}{2\pi i}\dfrac{d}{d\tau}=q\dfrac{d}{dq}.$

For $k\geq 0$, we observe that $D^{k}:q^{n}\mapsto n^{k}q^{n}$. Ramanujan
famously proved the identities (for example, see [14, Section 2.3])222These
identities are often stated for the normalized Eisenstein series
$E_{k}:=1+\dots$.

$\displaystyle DG_{2}=-2G_{2}^{2}+\dfrac{5}{6}G_{4},\ \ \
DG_{4}=-8G_{2}G_{4}+\dfrac{7}{10}G_{6},\ \ \
DG_{6}=-12G_{2}G_{6}+\dfrac{400}{7}G_{4}^{2}.$

These identities show that $D$ defines a map on the algebra of quasimodular
forms, which increases weights by $2$. Namely, we have that
$D:\widetilde{M}_{k}\to\widetilde{M}_{k+2}.$

###  2.2. Some prime-detecting quasimodular forms

Using elementary properties of $D,$ it is straightforward to produce
quasimodular forms whose Fourier coefficients detect the primes.

######  Lemma 2.1.

Given non-negative odd integers $k,\ell$ with $\ell>k$, consider the
quasimodular forms

$\displaystyle
f_{k,\ell}:=\bigl{(}D^{\ell}+1\bigr{)}G_{k+1}-\bigl{(}D^{k}+1\bigr{)}G_{\ell+1}.$

For $n\geq 2$, the $n$th Fourier coefficient of $f_{k,\ell}$ vanishes if and
only if $n$ is prime. Furthermore, all of the coefficients of $f_{k,\ell}$ are
non-negative.

######  Remarks.

(1) The assumption in Lemma 2.1 that $k$ and $\ell$ are odd is imposed to
ensure $f_{k,\ell}$ are quasimodular forms. In particular, the vanishing
property is true for any non-negative integers $k$ and $\ell$.

(2) Lemma 2.1 appears in an unpublished note by Lelièvre (see Lemma 2 of
[15]). For completeness we offer a streamlined proof here.

###### Proof.

By inspection, the coefficient of $q^{1}$ vanishes. Therefore, without loss of
generality we may assume that $n\geq 2$ is an integer. Thanks to (2.1), the
$n$th coefficient of $f_{k,\ell}$ is given by

$\displaystyle\sum_{\begin{subarray}{c}1<d<n\\\
d|n\end{subarray}}\left[\left(1+n^{\ell}\right)d^{k}-\left(1+n^{k}\right)d^{\ell}\right].$

This sum is empty and therefore vanishing if $n$ is prime, so we assume now
that $n$ is composite. Since $d\leq\frac{n}{2}$ must hold, we then have

$\displaystyle\left(1+n^{\ell}\right)d^{k}-\left(1+n^{k}\right)d^{\ell}$
$\displaystyle\geq
d^{k}\left(1+n^{\ell}-\left(\dfrac{n}{2}\right)^{\ell-k}-n^{k}\left(\dfrac{n}{2}\right)^{\ell-k}\right)$
$\displaystyle\geq
d^{k}\left(1+n^{\ell}\left(1-\dfrac{1}{2^{\ell-k-1}}\right)\right)\geq 1.$

Thus, the coefficients of $f_{k,\ell}$ are all non-negative. ∎

We now turn to the problem of determining those quasimodular forms whose
Fourier expansions detect primes. To this end, we let $\mathcal{E}$ denote the
_quasimodular Eisenstein space_ , the subspace of quasimodular forms generated
additively by the even weight Eisenstein series and their derivatives. We then
let $\Omega$ denote the prime-detecting quasimodular forms, where a
quasimodular form

$f=\sum_{n\geq 0}b_{n}(f)q^{n}$

is prime-detecting if for positive $n$ we have $b_{n}(f)\geq 0$, and for
$n\geq 2$ vanish if and only if $n$ is prime.

Based on numerical evidence, we expect that the property of vanishing of the
Fourier coefficients at primes is not satisfied by quasimodular forms which
are not in the quasimodular Eisenstein space.

######  Conjecture 2.2.

We have $\Omega=\mathcal{E}\cap\Omega.$

###  2.3. Characterizing $\mathcal{E}\cap\Omega$

We give an explicit description of $\mathcal{E}\cap\Omega$ (see Theorem 2.3).
To this end, for even $k\geq 6$, we define the distinguished quasimodular
forms of weight $k$

(2.3) $\displaystyle H_{k}:=\sum_{n\geq
0}b_{n}(H_{k})q^{n}:=\begin{cases}\frac{1}{6}\left((D^{2}-D+1)G_{2}-G_{4}\right)&\
\ \ {\text{\rm if}}\ k=6,\\\
\frac{1}{24}(-D^{2}G_{k-6}+(D^{2}+1)G_{k-4}-G_{k-2})&\ \ \ {\text{\rm if}}\
k\geq 8.\end{cases}$

The quasimodular forms $H_{k}$ give rise to all quasimodular forms which are
linear combinations of derivatives of Eisenstein series that are prime-
detecting.

######  Theorem 2.3.

For all $n\geq 0$ and $k\geq 6$, we have $D^{n}H_{k}\in\Omega$. Conversely, if
$f\in\mathcal{E}\cap\Omega$, then $f$ is a linear combination of the forms
$D^{n}H_{k}$ for $n\geq 0$ and $k\geq 6$.

Before we prove this theorem, we consider two additional results on properties
of these quasimodular forms $H_{k}$.

######  Lemma 2.4.

For $k\geq 6$ even, the following are true:
(1) We have that $b_{1}(H_{k})=0.$
(2) The Fourier coefficients $b_{n}(H_{k})$ for $n\geq 2$ are integers.
(3) The Fourier coefficients $b_{n}(H_{k}),$ where $n\geq 2,$ generate
$\mathbb{Z}$ as a ring.

###### Proof.

Claim (1) is clear by inspection. Therefore, we now assume that $n\geq 2$ is
an integer. Noting that $b_{4}(H_{6})=3$ and $b_{6}(H_{6})=20$ are relatively
prime, we see that (3) holds for $H_{6}$ assuming that (2) holds. To show (2)
for $H_{6}$, we must prove that

$(n^{2}-n+1)\sigma_{1}(n)-\sigma_{3}(n)\equiv 0\pmod{6}.$

Evenness is trivial, and so we consider this expression modulo 3. If
$n\not\equiv 2\pmod{3}$, then after multiplying by $\left(n+1\right)$, the
3-divisibility follows from

$\displaystyle\left(n^{3}+1\right)\sigma_{1}(n)-\left(n+1\right)\sigma_{3}(n)\equiv
0\pmod{3},$

which is an easy consequence of Fermat’s little theorem. If $n\equiv
2\pmod{3}$, then there is a prime $p\equiv 2\pmod{3}$ dividing $n$ with odd
multiplicity $2m+1$, and then 3-divisibility follows from multiplicativity of
divisor sums and $\sigma_{2j+1}\left(p^{2m+1}\right)\equiv 0\pmod{3}$. Thus,
both (2) and (3) hold for $H_{6}$.

We now consider even $k\geq 8.$ To prove (2), we must prove that

$\displaystyle-n^{2}\sigma_{k-7}(n)+\left(n^{2}+1\right)\sigma_{k-5}(n)-\sigma_{k-3}(n)\equiv
0\pmod{24}.$

It is straightforward to see this congruence modulo 3 from the fact that
$\sigma_{j+2}(n)\equiv\sigma_{j}(n)\pmod{3}$ for any $j$, and so we consider
the congruence above modulo 8. For any odd prime $p$, we have from
$p^{2}\equiv 1\pmod{8}$ that
$\sigma_{2j+1}\left(p^{m}\right)\equiv\sigma_{1}\left(p^{m}\right)\pmod{8}$
for any integer $j\geq 0$. Now, writing $n=2^{a}m$ for $m$ odd, it follows
from the multiplicativity of $\sigma_{j}$ that

(2.4) $\displaystyle-$ $\displaystyle
n^{2}\sigma_{k-7}(n)+\left(n^{2}+1\right)\sigma_{k-5}(n)-\sigma_{k-3}(n)$
$\displaystyle\ \
\equiv\left(-4^{a}m^{2}\sigma_{k-7}\left(2^{a}\right)+\left(4^{a}m^{2}+1\right)\sigma_{k-5}\left(2^{a}\right)-\sigma_{k-3}\left(2^{a}\right)\right)\sigma_{1}(m)\pmod{8}$
$\displaystyle\ \ \equiv
4^{a}m^{2}\left(\sigma_{k-5}\left(2^{a}\right)-\sigma_{k-7}\left(2^{a}\right)\right)\sigma_{1}(m)-\left(\sigma_{k-3}\left(2^{a}\right)-\sigma_{k-5}\left(2^{a}\right)\right)\sigma_{1}(m)\pmod{8}.$

Now, for $j\geq 1$ we observe that

$\displaystyle\sigma_{2j+1}\left(2^{a}\right)-\sigma_{2j-1}\left(2^{a}\right)=\sum_{b=0}^{a}\left(2^{(2j+1)b}-2^{(2j-1)b}\right)=\sum_{b=0}^{a}2^{(2j-1)b}\left(4^{b}-1\right).$

It is therefore clear that if $j\geq 2$,
$\sigma_{2j+1}\left(2^{a}\right)-\sigma_{2j-1}\left(2^{a}\right)\equiv
0\pmod{8}$. If then $j=1$, from the fact that
$\sigma_{3}\left(2^{a}\right)-\sigma_{1}\left(2^{a}\right)$ is even, we obtain
from (2.4) that

$\displaystyle-n^{2}\sigma_{k-7}(n)$
$\displaystyle+\left(n^{2}+1\right)\sigma_{k-5}(n)-\sigma_{k-3}(n)\equiv-\left(\sigma_{k-3}\left(2^{a}\right)-\sigma_{k-5}\left(2^{a}\right)\right)\sigma_{1}(m)\pmod{8}.$

But then for $2j+1=k-3$, the condition $k\geq 8$ corresponds to the condition
$j\geq 2$, and so we obtain (2) for all $H_{k}$.

We now prove (3). Let $g$ be the largest integer dividing

$24b_{n}(H_{k})=-n^{2}\sigma_{k-7}(n)+\left(n^{2}+1\right)\sigma_{k-5}(n)-\sigma_{k-3}(n)$

for all $n\geq 2$. We know $24\mid g$, and therefore it will be enough to show
that each of $16$, $9$, and $p$ for primes $p\geq 5$ do not divide $g$.
Beginning with any $p\geq 3$ prime and $m\geq 1$ coprime to $p$, it is easy to
see that

$\displaystyle
24b_{pm}(H_{k})\equiv\sigma_{k-5}(pm)-\sigma_{k-3}(pm)\equiv\sigma_{k-5}(m)-\sigma_{k-3}(m)\pmod{p}.$

We then observe that $24b_{2p}(H_{k})\not\equiv 0\pmod{p}$ for $p\geq 5$ and
$b_{5p}(H_{k})\not\equiv 0\pmod{9}$. Finally, for $m$ odd we see that

$\displaystyle
24b_{4m}(H_{k})\equiv\sigma_{k-5}(4)\sigma_{k-5}(m)-\sigma_{k-3}(4)\sigma_{k-3}(m)\pmod{16}.$

It is then easy to see that this expression is not a multiple of 16 for $m=1$
in the case $k=8$ and for $m=3$ if $k\geq 10$. Therefore, we have shown that
$g\mid 24$, so that $g=24$ as required. ∎

To prove Theorem 2.3, we also require the following about the linear
independence of the forms $D^{n}H_{\ell}$ and the spaces they generate.

######  Lemma 2.5.

The following are true.
(1) The quasimodular forms $D^{n}H_{\ell}$ (with $n\geq 0$ and $\ell\geq 6$)
are linearly independent.
(2) The subspace of elements of the $D^{n}H_{\ell}$ of weight $\leq k$ is of
dimension $\frac{1}{8}(k-2)(k-4)$.

###### Proof.

We first prove (1). We show the stronger statement that the quasimodular forms
$D^{n}E_{\ell}$, generating $\mathcal{E}$, are linearly independent. Observe
that any non-trivial relations between the quasimodular forms $D^{n}H_{\ell}$
would imply a relation between these derivatives of Eisenstein series.

To prove the claim, first of all, observe that one can restrict to
quasimodular forms of homogeneous weight. Namely, by the quasimodular
transformation behaviour one recovers the weight. For quasimodular forms of
homogeneous weight, the statement follows from [16, Proposition 20 (iii)]

$\widetilde{M}_{k}=\bigoplus_{r=0}^{k/2-1}D^{r}M_{k-2r}\oplus\mathbb{C}D^{k/2-1}E_{2}.$

For (2), recall that the operator $D$ increases the weight by $2$ and
$H_{\ell}$ is of weight $\ell$. We count the number of linear independent
generators of weight $\leq k$, as follows:

$\displaystyle\sum_{\ell\geq 3}\,\sum_{2\ell+2n\leq
k}\\!1\,=\,\frac{1}{8}(k-2)(k-4).$

###### Proof of Theorem 2.3.

If $f$ is a quasimodular form, then $f\in\Omega$ if and only if $Df\in\Omega.$
Therefore, it suffices to check that $H_{k}\in\Omega$ for $k\geq 6$. We first
observe that for $n\geq 2$, we have

$\displaystyle(n+1)b_{n}(H_{6})=\dfrac{n+1}{6}\left[\left(n^{2}-n+1\right)\sigma_{1}(n)-\sigma_{3}(n)\right]=\dfrac{1}{6}\left[\left(n^{3}+1\right)\sigma_{1}(n)-\left(n+1\right)\sigma_{3}(n)\right].$

This gives the identity $(D+1)H_{6}=\frac{1}{6}f_{1,3}$, and so Lemma 2.1
implies that $H_{6}\in\Omega$.

For $k\geq 8$, we calculate for $n\geq 2$ that

$\displaystyle b_{n}(H_{k})$
$\displaystyle=\sum_{d|n}\left[-n^{2}d^{k-7}+\left(n^{2}+1\right)d^{k-5}-d^{k-3}\right]=\sum_{\begin{subarray}{c}1<d<n\\\
d|n\end{subarray}}d^{k-7}\left(d^{2}-1\right)\left(n^{2}-d^{2}\right).$

By inspection, we see that $H_{k}\in\Omega$.

For the converse, let $f\in\mathcal{E}\cap\Omega$ be of weight $\leq k$. As
$D$ increases the weight by $2$, we can write

$f\,=\,\sum_{\ell=1}^{k/2}p_{2\ell}(D)\,G_{2\ell}$

with $p_{2\ell}$ a polynomial of degree at most $\frac{1}{2}k-\ell$. The
vanishing of the Fourier coefficients at all primes is equivalent to

(2.5) $\sum_{\ell=1}^{k/2}p_{2\ell}(x)(x^{2\ell-1}+1)=0.$

(Initially, this equality only holds for primes $x$, and then, since a
polynomial of bounded degree is uniquely determined by its values at primes,
this is an equality of polynomials in the variable $x$.) Observe the dimension
of the space of sequences of polynomials $p_{2},p_{4},\ldots,p_{k}$ with the
given degree bound equals

$\sum_{\ell=1}^{k/2}\left(\frac{1}{2}k-\ell+1\right)\,=\,\frac{1}{8}k(k+2).$

Note that the left-hand side of (2.5) is of degree $\leq k-1$. Hence,
condition (2.5) gives $\leq k$ linearly independent constraints between the
coefficients of the polynomials. For notational convenience, we let

$p_{2\ell}(x)=\sum_{j=0}^{k/2-\ell}a_{j,2\ell}\,x^{j}$

for real coefficients $a_{j,2\ell}$, and by $M$ denote the matrix of linear
relations of the $a_{j,2\ell}$ given by (2.5), where the rows are
parameterized by the degrees $0,\ldots,k-1$ and the columns by pairs
$(j,2\ell)$, ordered lexicographically.

We show that, after removing the $0$th row of $M$, the matrix is (up to
permuting the columns) in row echelon form. Namely, observe that in the $i$th
row the smallest pair $(j,2\ell)$ for which the corresponding entry in $M$ is
non-zero is given by

$\begin{cases}(0,2)&i=0,\\\ (i+1,0)&i>0\text{ odd},\\\ (i,1)&i>0\text{
even}.\end{cases}$

Hence, the matrix obtained from $M$ by removing the $0$th row is of full rank.
We conclude the rank of $M$ is at least $k-1$. As (2.5) is divisible by $x+1$
one can easily infer that the rank is exactly ${k-1}$. Therefore, the
dimension of the space of polynomials satisfying (2.5) is given by

$\sum_{\ell=1}^{k/2}\left(\frac{1}{2}k-\ell+1\right)-(k-1)=\frac{1}{8}(k-2)(k-4).$

Accordingly, this is also the dimension of the vector space of elements
$D^{n}H_{2\ell}$ of weight $\leq k$ by Lemma 2.5(2). ∎

##  3\. Quasimodularity of $\mathcal{U}_{a}(q)$ and the proofs Theorems 1.1
and 1.2

Here we use Theorem 2.3 to prove Theorems 1.1 and 1.2.

###  3.1. Proof of Theorem 1.1

MacMahon’s series $\mathcal{U}_{a}(q)$ have been the focus of much study (for
example, see [13, 17, 18, 19, 20]). The first glimpses of their beautiful
properties are captured by the identities (for example, see [8, 17])

$\mathcal{U}_{1}(q)=\sum_{n\geq 1}\sigma_{1}(n)q^{n}\ \ \ \ {\text{\rm and}}\
\ \ \ \mathcal{U}_{2}(q)=\frac{1}{8}\sum_{n\geq
1}\left((-2n+1)\sigma_{1}(n)+\sigma_{3}(n)\right)q^{n}.$

To prove (1), we note that the sequence we wish to study is the coefficients
of the series $F$ given by

$\displaystyle F:=\left(D^{2}-3D+2\right)\mathcal{U}_{1}(q)-8\
\mathcal{U}_{2}(q)=\sum_{n\geq
1}\left[\left(n^{2}-n+1\right)\sigma_{1}(n)-\sigma_{3}(n)\right]q^{n}.$

Noting the identities $F=6H_{6}$ and $\left(D+1\right)F=f_{1,3}$, part (1) of
Theorem 1.1 then follows by either Theorem 2.3 or Lemma 2.1, respectively.

To prove (2), we begin by noting that [8, p. 102]

$\displaystyle\mathcal{U}_{3}(q)=\dfrac{1}{1920}\sum_{n\geq
1}\left[\left(40n^{2}-100n+37\right)\sigma_{1}(n)-10(3n-5)\sigma_{3}(n)+3\sigma_{5}(n)\right]q^{n}.$

The sequence we wish to study form the coefficients of the series $G$ given by

$\displaystyle G$
$\displaystyle:=\left(3D^{3}-13D^{2}+18D-8\right)\mathcal{U}_{1}(q)-\left(12D^{2}+120D+212\right)\mathcal{U}_{2}(q)-960\
\mathcal{U}_{3}(q).$

It is easy to see using previous formulas for $\mathcal{U}_{j}(q)$ that
$G=36H_{8}$, and so (2) follows by Theorem 2.3.

###  3.2. Proof of Theorem 1.2

We recall the result of Andrews and Rose [12, Corollary 4] that
$\mathcal{U}_{a}(q)\in\widetilde{M}_{\leq 2a}$. By Theorem 2.3, any linear
combination of derivatives of various $\mathcal{U}_{a}(q)$ detects primes if
it can be represented by linear combinations of derivatives of the
quasimodular forms $H_{k}$. Since the spaces of quasimodular forms of (mixed)
weight $\leq 2a$ is finite-dimensional, the examples given in the Appendix can
be verified by a straightforward linear algebra calculation333In particular,
this can be done using a row-echelon basis for the spaces $\widetilde{M}_{\leq
a}$ with rational coefficients.. In particular, the examples in Table 1 are
exactly
$6H_{6}\mspace{1.0mu},36H_{8}\mspace{1.0mu},90H_{10}\mspace{1.0mu},90H_{12}\mspace{1.0mu}$,
and $30\left(D^{2}+1\right)H_{14}+30H_{16}$ respectively.

##  4\. MacMahonesque series and the proof of Theorem 1.3 and Theorem 1.4

###  4.1. Certain $q$-analogues of multiple zeta values

The proofs in the previous section rely on the fact that the $\mathcal{U}_{a}$
are quasimodular forms. Similarly, for odd $a$, we have that
$\mathcal{U}_{(a)}(q)$ is quasimodular, as we have

(4.1) $\mathcal{U}_{(a)}(q)=\sum_{n\geq
1}M_{(a)}(n)q^{n}=\frac{B_{a+1}}{2(a+1)}+G_{a+1}(\tau).$

The proof of Theorem 1.3 relies on a refinement of this phenomenon. Although
the generating functions

(4.2) $\mathcal{U}_{\vec{a}}(q):=\sum_{n=1}^{\infty}M_{\vec{a}}(n)q^{n}$

are not always a quasimodular forms, it turns out that they arise naturally in
the study of $q$-analogues of multiple zeta values (see [10, Section 6]),
which often turn out to be related to quasimodular forms.

For completeness, before we proceed with the tools for the proofs of Theorem
1.3 and Theorem 1.4, we consider the $q$-series generating functions for the
MacMahonesque $M_{\vec{a}}(n)$ that generalize (1.1).

######  Lemma 4.1 ([10, Lemma 2.5]).

If $\vec{a}\in\mathbb{N}^{a},$ then we have that

$\mathcal{U}_{\vec{a}}(q)=\sum_{n=1}^{\infty}M_{\vec{a}}(n)q^{n}=\sum_{0<s_{1}<s_{2}<\cdots<s_{a}}\frac{q^{s_{1}+s_{2}+\dots+s_{a}}P_{v_{1}}(q^{s_{1}})P_{v_{2}}(q^{s_{2}})\cdots
P_{v_{a}}(q^{s_{a}})}{(1-q^{s_{1}})^{v_{1}+1}(1-q^{s_{2}})^{v_{2}+1}\cdots(1-q^{s_{a}})^{v_{a}+1}},$

where $P_{k}(x)$ is the $k$th Eulerian polynomial, which are defined by

$\displaystyle\dfrac{xP_{n}(x)}{\left(1-x\right)^{n+1}}:=\sum_{k\geq
1}k^{n}x^{k}.$

######  Remark.

In analogy with (1.1), one might think that it makes sense to make use of the
more aesthetically pleasing $q$-series

$\mathcal{V}_{\vec{a}}(q)=\sum_{n\geq
1}N_{\vec{a}}(n)q^{n}:=\sum_{0<s_{1}<s_{2}<\cdots<s_{a}}\frac{q^{s_{1}+s_{2}+\cdots+s_{a}}}{(1-q^{s_{1}})^{v_{1}+1}(1-q^{s_{2}})^{v_{2}+1}\cdots(1-q^{s_{a}})^{v_{a}+1}}.$

However, these $q$-series are the generating functions for the less pleasing
partition functions

$N_{\vec{a}}(n)=\sum_{\begin{subarray}{c}0<s_{1}<s_{2}<\dots<s_{a}\\\
n=m_{1}s_{1}+m_{2}s_{2}+\dots+m_{a}s_{a}\end{subarray}}\binom{m_{1}+v_{1}-1}{v_{1}}\binom{m_{2}+v_{2}-1}{v_{2}}\cdots\binom{m_{a}+v_{a}-1}{v_{a}}.$

To derive this expression, one notes for non-negative integers $v$ (and
$|q|<1$) that

$\frac{q}{(1-q)^{v+1}}=\sum_{m=1}^{\infty}\binom{m+v-1}{v}q^{m}.$

All of our results regarding $M_{\vec{a}}(n)$ could be reformulated in terms
of $N_{\vec{a}}(n)$. Finally, we note that the natural definition of
$M_{\vec{a}}(n)$ in (1.3) necessitates the Eulerian polynomials in Lemma 4.1.
For convenience, we offer the first few polynomials

$P_{0}(x)=1,\ \ P_{1}(x)=1,\ \ P_{2}(x)=x+1,\ \ P_{3}(x)=x^{2}+4x+1,\ \
P_{4}(x)=x^{3}+11x^{2}+11x+1,\dots.$

For $\vec{a}=(1,1,\dots,1),$ we recover
$\mathcal{U}_{\vec{a}}(q)=\mathcal{V}_{\vec{a}}(q)$ as $P_{1}(x)=1.$

Write $\mathcal{Z}_{q}$ for the vector space generated additively by the
MacMahonesque series $\mathcal{U}_{\vec{a}}$. It turns out that
$\mathcal{Z}_{q}$ is strictly larger than the space of quasimodular forms.
Furthermore, this space enjoys two important properties; namely, (i) it is an
algebra, and (ii) it is closed under differentiation by $D=q\frac{d}{dq}$.
Using the theory of quasi-shuffle algebras, Bachmann and Kühn proved the
following theorem [10, Theorem 1.3 and 1.7].

######  Theorem 4.2.

The space $\mathcal{Z}_{q}$ is a differential subalgebra of
$\mathbb{Q}[\\![q]\\!]$ containing $\widetilde{M}$.

Observe that this theorem directly implies that any product of (in fact, any
polynomial in) the MacMahon series is a linear combination of MacMahonesque
series $\mathcal{U}_{\vec{a}}$. Similar, any derivative of a MacMahon series
is a linear combination of the MacMahonesque series $\mathcal{U}_{\vec{a}}$.
Expressing these observations in terms of the MacMahon functions, we obtain
Theorem 1.4.

The proof of this theorem allows one to compute the product of two
MacMahoneque series recursively, as we explain now. Let
$A=\\{z_{1},z_{2},\ldots\\}$ be the set of _letters_ $z_{j}$, one for each
positive integer. Furthermore, we let $\mathbb{Q}A$ be the $\mathbb{Q}$-vector
space generated by these letters. A _word_ is a non-commutative product of
letters, and we write $\mathbb{Q}\langle A\rangle$ for the (non-commutative)
polynomial algebra over $\mathbb{Q}$ generated by words with letters in $A$.
On $\mathbb{Q}A$, we define a commutative and associative product $\diamond$
by444 Since our $U_{\vec{a}}$ correspond to
$\left(\prod_{i}v_{i}!\right)[\nu_{1}+1,\ldots,\nu_{a}+1]$ in [10], we chose a
different normalization here.

(4.3) $z_{i}\diamond z_{j}:=\sum_{m=0}^{i+j}c_{i,j,m}\,z_{m+1}\,,$

where we let

$c_{i,j,{i+j}}:=\frac{i!j!}{(i+j+1)!},$

and for $m<i+j$ we let

$c_{i,j,m}:=\left[(-1)^{i}\binom{i}{m+1}+(-1)^{j}\binom{j}{m+1}\right]\frac{B_{i+j-m}}{i+j-m}.$

Then, for $x,y\in A$ and $w,v\in\mathbb{Q}\langle A\rangle$, we recursively
define the _quasi-shuffle product_ by $1*w=w*1=w$ and

$xw*yv:=x(w*yv)+y(xw*v)+(x\diamond y)(w*v).$

Identify $\mathbb{Q}\langle A\rangle$ with the vector space generated by
finite sequences of positive integers, i.e. $z_{v_{1}}\cdots z_{v_{a}}$
corresponds to $\vec{a}=(v_{1},\ldots,v_{a})$. Under this identification, the
quasi-shuffle product yields, e.g.,

$(1)*(1,1)=3(1,1,1)+\frac{1}{6}(3,1)+\frac{1}{6}(1,3)-\frac{1}{3}(1,1).$

We stress that this is a formal sum of vectors, which cannot be simplified by
coordinate-wise addition. Moreover, we extend the generalized MacMahon series
linearly, i.e., for a formal finite linear combinations of finite sequences
$\sum_{\vec{a}}c_{\vec{a}}\,\vec{a}$ we set

(4.4)
$\mathcal{U}_{\sum_{\vec{a}}c_{\vec{a}}\,\vec{a}}:=\sum\nolimits_{\vec{a}}c_{\vec{a}}\,\mathcal{U}_{\vec{a}}.$

Then by [10, Proposition 2.10], the power series multiplication of the
generalized MacMahon series can be computed recursively using the quasi-
shuffle product defined above.

######  Proposition 4.3.

For all $\vec{a}\in\mathbb{N}^{a}$ and $\vec{b}\in\mathbb{N}^{b}$, one has

$\mathcal{U}_{\vec{a}}(q)\,\mathcal{U}_{\vec{b}}(q)=\mathcal{U}_{{\vec{a}}*{\vec{b}}}(q).$

Moreover, assigning weight $|\vec{a}|+\ell(\vec{a})$ to $\vec{a}$, the algebra
$\mathcal{Z}_{q}$ becomes a filtered algebra.

######  Example.

As an example of Proposition 4.3, we have

$\mathcal{U}_{(1)}(q)\,\mathcal{U}_{(1,1)}(q)=3\mathcal{U}_{(1,1,1)}(q)+\frac{1}{6}\mathcal{U}_{(3,1)}(q)+\frac{1}{6}\mathcal{U}_{(1,3)}(q)-\frac{1}{3}\mathcal{U}_{(1,1)}(q)=\mathcal{U}_{(1)*(1,1)}(q),$

which is equivalent to

$\mathcal{U}_{1}(q)\,\mathcal{U}_{2}(q)=3\mathcal{U}_{3}(q)+\frac{1}{6}\mathcal{U}_{(3,1)}(q)+\frac{1}{6}\mathcal{U}_{(1,3)}(q)-\frac{1}{3}\mathcal{U}_{2}(q).$

This example illustrates the fact that the vector space generated by the
MacMahon series is not closed under multiplication, but $\mathcal{Z}_{q}$ is.

###  4.2. Symmetric group action

As we have emphasized above, a generic MacMahonesque $q$-series
$\mathcal{U}_{\vec{a}}(q)$ is not quasimodular. However, their “symmetrized
sums” are quasimodular, and provide canonical $q$-series that additively
generate the ring of quasimodular forms.

To make this precise, we recall that the generalized MacMahon series
$\mathcal{U}_{\vec{a}}(q)$ admit the natural action of the $a$-th symmetric
group $S_{a}$ by

(4.5)
$\displaystyle\sigma\mathcal{U}_{\vec{a}}(q):=\mathcal{U}_{\sigma\vec{a}}(q),$

where for $\sigma\in S_{a}$ we let
$\sigma\vec{a}:=\left(v_{\sigma(1)},\dots,v_{\sigma(a)}\right)$. For those
$\vec{a}$ consisting of only odd entries, we consider the symmetrized
MacMahonesque $q$-series

(4.6) $\mathcal{U}_{\vec{a}}^{\operatorname{sym}}(q):=\sum_{\sigma\in
S_{a}}\sigma\mathcal{U}_{\vec{a}}(q).$

######  Theorem 4.4.

The following are true.

(1) For all $\vec{a}\in(2\mathbb{Z}_{+}-1)^{a}$, we have that
$\mathcal{U}_{\vec{a}}^{\operatorname{sym}}(q)\in\widetilde{M}_{\leq|\vec{a}|+a}.$

(2) The $\mathcal{U}_{\vec{a}}^{\operatorname{sym}}(q)$ with $|\vec{a}|+a\leq
k$ additively generate $\widetilde{M}_{\leq k}$.

######  Remark.

Theorem 4.4 (1) is given as an exercise in [11]. Part (2) of the theorem is
the new observation.

###### Proof.

We now prove (1). Symmetric sums in the context of quasi-shuffle algebras have
been studied broadly. In particular, by [21, Theorem 13] we have

(4.7) $\mathcal{U}_{\vec{a}}^{\operatorname{sym}}(q)=\sum_{\sigma\in
S_{a}}\sigma\mathcal{U}_{\vec{a}}(q)=\sum_{B\in\Pi_{a}}c(B)\prod_{\beta\in
B}\mathcal{U}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\diamond_{\beta}}}(q),$

where $\Pi_{a}$ denotes all set partitions of the set $\\{1,\ldots,a\\}$,
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\diamond_{\beta}:=\diamond_{i\in\beta}(v_{i})}$
denotes the $\diamond$-product of the $z_{v_{i}}$ (written as a finite
sequence), and

$c(B)=(-1)^{a-|B|}\prod_{\beta\in B}(|\beta|-1)!$

is a constant. In particular, this expresses $\sum_{\sigma\in
S_{a}}\sigma\mathcal{U}_{\vec{a}}(q)$ as a polynomial in the modular
Eisenstein series (up to constant terms) $\mathcal{U}_{(v)}$ for certain
$v\in\mathbb{Z}_{+}$. We now show that, in fact, $v$ only takes odd values.

Recall that the entries $v_{1},\ldots,v_{a}$ of $\vec{a}$ are all odd. Now,
observe that for odd $i,j$, we have that $c_{i,j,m}=0$ if $m$ is odd. Namely,
for odd $i,j,m$ the Bernoulli number $B_{i+j-m}$ is zero, unless $i+j-m=1$. In
the latter case $m+1>i$ and $m+1>j$, and hence
$\binom{i}{m+1}=\binom{j}{m+1}=0$. Hence, in the $\diamond$-product of $z_{i}$
and $z_{j},$ where $i$ and $j$ are odd, we find that the coefficients of
$z_{m+1}$ vanish if $m$ is odd. We conclude that $\sum_{\sigma\in
S_{a}}\sigma\mathcal{U}_{\vec{a}}(q)$ can be expressed as a polynomial in
$\mathcal{U}_{(v)}$ for certain odd integers $v$. Now, (1) follows from (4.1).

We now turn to the proof of (2). By (4.1), the symmetrized MacMahonesque
$q$-series generate all the Eisenstein series. If $\vec{a}=(v_{1},v_{2})$ is
of length $2$, then by brute force using (4.7) we find that
$\mathcal{U}_{\vec{a}}^{\operatorname{sym}}(q)$ is—up to an Eisenstein
series—equal to the product of the Eisenstein series
$\mathcal{U}_{(v_{1})}(q)$ and $\mathcal{U}_{(v_{2})}(q)$. It follows that
$\mathcal{U}_{\vec{a}}^{\operatorname{sym}}(q)$ with $|\vec{a}|+a\leq k$
additively generate the algebra of modular forms $M_{\leq k}$, as it is well-
known that the space of modular forms is generated by polynomials in
Eisenstein series of degree at most $2$ (for example, see [22, Section 5]).
Finally, using (4.7) again, observe that for $\vec{a}=(v,1,1,\ldots,1)$ and
$\vec{a}=(v_{1},v_{2},1,1,\ldots,1)$ one obtains arbitrary powers of the
quasimodular $\mathcal{U}_{1}$ times Eisenstein series and products of
Eisenstein series respectively, modulo quasimodular forms of lower depth.
Here, we recall every quasimodular form can be written as a polynomial in the
quasimodular Eisenstein series $\mathcal{U}_{1}$ (which equals $G_{2}$ up to a
constant), and the depth is defined as the degree in $\mathcal{U}_{1}$. Claim
(2) now follows by induction on the depth. ∎

###  4.3. Proof of Theorem 1.3

For $d=4$, we have example $\Psi_{1}(n)$ in the Introduction. Multiplying by
$n^{d-4}$ and using Theorem 1.4 implies the existence for all $d\geq 4$,
thereby confirming (1).

We turn to (2). By Theorem 4.4 (2) the subspace of all
$\mathcal{U}_{\vec{a}}(q)$ of weight at most $k$ generates
$\widetilde{M}_{k}$. By Lemma 2.5 there are $\gg k^{2}$ linearly independent
quasimodular forms of weight at most $k$ for which the $n$th Fourier
coefficient ($n\geq 2$) is non-negative and vanish if and only if $n$ is
prime. We also note that the coefficient of $q^{1}$ in these forms also
vanish. This proves (2).

###  4.4. Proof of Theorem 1.4

The first part of the statement follows from Proposition 4.3. In particular,
the bound on $\vec{a}$ in the sum in the statement of the theorem, follows
from estimating $\ell(\vec{a})\geq 1$ and the fact that $\mathcal{Z}_{q}$ is a
filtered algebra. The second part of the statement follows from the fact that
the operator $D$ preserves $\mathcal{Z}_{q}$ and increases the weight by $2$
(see Theorem 4.2 and [10, Theorem 1.7]).

## Appendix

Table 1. Prime detecting expressions for Theorem 1.2. See (2.3) for the definition of the quasimodular forms $H_{k}$. Prime Detecting Expression | Quasimodular form
---|---
$\displaystyle(n^{2}-3n+2)M_{1}(n)-8M_{2}(n)$ | $6H_{6}$
$\displaystyle(3n^{3}-13n^{2}+18n-8)M_{1}(n)+(12n^{2}-120n+212)M_{2}(n)-960M_{3}(n)$ | $36H_{8}$
$\displaystyle(25n^{4}-171n^{3}+423n^{2}-447n+170)M_{1}(n)+(300n^{3}-3554n^{2}+12900n-14990)M_{2}(n)+$ |
$\displaystyle\quad(2400n^{2}-60480n+214080)M_{3}(n)-725760M_{4}(n)$ | $90H_{10}$
$\displaystyle(126n^{5}-1303n^{4}+5073n^{3}-9323n^{2}+8097n-2670)M_{1}(n)+$ |
$\displaystyle\quad(3024n^{4}-48900n^{3}+288014n^{2}-737100n+695490)M_{2}(n)+$ |
$\displaystyle\quad(60480n^{3}-1510080n^{2}+10644480n-23496480)M_{3}(n)+$ |
$\displaystyle\quad(725760n^{2}-36288000n+218453760)M_{4}(n)-580608000M_{5}(n)$ | $90H_{12}$
$\scriptstyle\quad(300n^{8}-1542n^{7}-33049n^{6}+377959n^{5}-1651959n^{4}+3726801n^{3}-4575760n^{2}+2903750n-746500)M_{1}(n)+$ |
$\scriptstyle\quad(12000n^{7}-91008n^{6}-2799900n^{5}+50637162n^{4}-351366300n^{3}+1239098170n^{2}-2210467000n+1585493500)M_{2}(n)+$ |
$\scriptstyle\quad(432000n^{6}-3548160n^{5}-236343840n^{4}+5133219840n^{3}-42370071840n^{2}+161101416000n-236150560800)M_{3}(n)+$ |
$\scriptstyle\quad(12096000n^{5}-72817920n^{4}-17599680000n^{3}+396192142080n^{2}-3123876672000n+8555162112000)M_{4}(n)+$ |
$\scriptstyle\quad(193536000n^{4}-1056513024000n^{2}+21310248960000n-112944125952000)M_{5}(n)+$ |
$\scriptstyle\quad(-46495088640000n+604436152320000)M_{6}(n)-1115882127360000M_{7}(n)$ | $\scriptstyle 30(D^{2}+1)H_{14}+30H_{16}$

## References

  * [1] Y. V. Matiyasevich, _Enumerable sets are Diophantine_ , Doklady Akademii Nauk SSSR (in Russian), 191 (1970), 279-282.
  * [2] J. P. Jones, D. Sato, H. Wada, and D. Wiens, _Diophantine representation of the set of prime numbers_ , American Mathematical Monthly 83 (1976), 449-464.
  * [3] G. E. Andrews, The theory of partitions, Reprint of the 1976 original, Cambridge Mathematical Library, Cambridge University Press, Cambridge, 1998.
  * [4] R. Schneider, _Eulerian series, zeta functions, and the arithmetic of partitions_ , Thesis (Ph.D.)- Emory University, 2018. 220 pages.
  * [5] M. Locus Dawsey, M. Just and R. Schneider, _A “supernormal” partition statistic_ , J. Number Theory 241 (2022), 120-141.
  * [6] R. Schneider and A. V. Sills, _Analysis and combinatorics of partition zeta functions_ , Int. J. Number Theory 17, No. 3 (2021), 805-814.
  * [7] R. Schneider, _Partition zeta functions_ , Research in Number Theory 2 (2016), Article 9.
  * [8] P. A. MacMahon, _Divisors of Numbers and their Continuations in the Theory of Partitions_ , Proc. London Math. Soc. (2) 19 (1920), no.1, 75-113 [also in Percy Alexander MacMahon Collected Papers, Vol.2, pp. 303–341 (ed. G.E. Andrews), MIT Press, Cambridge, 1986].
  * [9] H. Bachmann, _Multiple Zeta-Werte und die Verbindung zu Modulformen durch Multiple Eisensteinreihen_ , Masterarbeit in Mathematik, Universität Hamburg, 2012.
  * [10] H. Bachmann and U. Kühn, _The algebra of generating functions for multiple divisor sums and applications to multiple zeta values_ ,. Ramanujan J. 40 (2016), no.3, 605–648.
  * [11] H. Bachmann, _Lectures on Multiple zeta values and modular forms_ , Nagoya University, Spring 2020 https://www.henrikbachmann.com/uploads/7/7/6/3/77634444/mzv_mf_2020_v_5_4.pdf.
  * [12] G. E. Andrews and S. C. F. Rose, _MacMahon’s sum-of-divisors functions, Chebyshev polynomials, and quasimodular forms_ , J. Reine Angew. Math. 676 (2013), 97–103.
  * [13] S. C. F. Rose, _Quasimodularity of generalized sum-of-divisors functions_ , Research in Number Theory 1 (2015), Paper No. 18, 11 pp.
  * [14] K. Ono, _The web of modularity: arithmetic of the coefficient of modular forms and $q$-series_. CBMS Regional Conference Series in Mathematics, 102, Amer. Math. Soc., Providence, RI 2004.
  * [15] S. Lelièvre, _Quasimodular forms with Fourier coefficients zero at primes_ , unpublished note.
  * [16] Don Zagier, _Elliptic modular forms and their applications_ , In The 1-2-3 of modular forms, Universitext, pages 1–103. Springer, Berlin, 2008.
  * [17] T. Amdeberhan, K. Ono and A. Singh, _MacMahon’s sums-of-divisors and allied $q$-series_, https://arxiv.org/abs/2311.07496, preprint.
  * [18] T. Amdeberhan, G. E. Andrews, and R. Tauraso, _Extensions of MacMahon’s sums of divisors_ , Res. Math. Sci. 11, 8 (2024).
  * [19] H. Bachmann, MacMahon’s sums-of-divisors and their connection to multiple Eisenstein series, arXiv:2312.06132
  * [20] K. Ono and A. Singh, _Remarks on MacMahon’s $q$-series_, https://arxiv.org/abs/2402.08783, preprint.
  * [21] M. E. Hoffman, _Quasi-shuffle algebras and applications_ , Algebraic combinatorics, resurgence, moulds and applications (CARMA). Vol. 2, 327–348. IRMA Lect. Math. Theor. Phys., 32, EMS Publishing House, Berlin, 2020.
  * [22] D. Zagier, Modular forms whose Fourier coefficients involve zeta-functions of quadratic fields. In: Serre, J.-P., Zagier, D.B. (eds.), modular functions of one variable, VI (proceedings international conference, university Of Bonn, Bonn


# A new symmetric linearly implicit exponential integrator preserving
polynomial invariants or Lyapunov functions for conservative or dissipative
systems

Lu Li

###### Abstract

We present a new symmetric linearly implicit exponential integrator that
preserves the polynomial first integrals or the Lyapunov functions for the
conservative and dissipative stiff equations, respectively. The method is
tested by both oscillated ordinary differential equations and partial
differential equations, e.g., an averaged system in wind-induced oscillation,
the Fermi–Pasta–Ulam systems, and the polynomial pendulum oscillators. The
numerical simulations confirm the conservative properties of the proposed
method and demonstrate its good behavior in superior running speed when
compared with fully implicit schemes for long-time simulations.

###### keywords:

Linearly implicit , energy-preserving , exponential integrator , conservative
system , dissipative system

[inst1]organization=Machine Intelligence Department,addressline=Simula
Metropolitan Center for Digital Engineering, city=Oslo, postcode=0167,
country=Norway

##  1 Introduction

This paper focuses on the semilinear systems of the form

$\dot{y}(t)=Ay+f(y),\quad y(t_{0})=y_{0},$ (1)

where $A$ is a linear unbounded differential operator or a matrix which has
eigenvalues with large negative real part or with purely imaginary eigenvalues
of large modulus, and the non-linear term $f$ is supposed to be nonstiff
satisfying Lipschitz condition. The semilinear system (1) arises in many
applications, such as the charged-particle dynamic [1], the rapidly rotating
shallow water equations for semi-geostrophic particle motion [2], and the
semi-discretization of semilinear PDEs. Equation (1) is usually stiff and one
popular class of numerical integrators that are suitable for such problems is
the exponential integrators [3]. These methods normally permit larger step
sizes and provide higher accuracy than the non-exponential integrators. The
basic idea behind such methods is to solve the stiff part with an exact
solver. Based on the variation of constants formula, the exact solution of (1)
is given by

$y(t_{0}+h)=\text{exp}(hA)y_{0}+h\int_{0}^{1}\text{exp}((1-\tau)hA)f(y(t_{0}+\tau
h))d\tau,$ (2)

where the integration interval is $[t_{0},t_{0}+h]$. Most exponential
integrators can be obtained from an appropriate approximation of the integral
in exact solution (2). For example, the exponential Euler method is obtained
by interpolating the nonlinear term at $y_{0}$ with the form

$y_{1}=\text{exp}(hA)y_{0}+h\phi(hA)f(y_{0}),$

and the implicit exponential Euler method is given by interpolating the
nonlinear term at $y_{1}$ with the form

$y_{1}=\text{exp}(hA)y_{0}+h\phi(hA)f(y_{1}),$

where $\phi(z):=\frac{e^{z}-1}{z}$ [3]. More examples of exponential
integrators can be found in, e.g., [3, 4].

Equation (1) might possess important geometric structures. In particular, the
canonical Hamiltonian structure corresponds to equations of the form

$\dot{y}(t)=J\nabla H(y),\quad y(t_{0})=y_{0},$ (3)

where

$H(y)=\frac{1}{2}y^{T}My+U(y),$

and

$J=\left[\begin{matrix}0&I_{m}\\\ -I_{m}&0\end{matrix}\right].$

with $I_{m}$ an identity matrix, $M$ a symmetric matrix and $U(y)$ a scalar
function. Two prominent features of equation (3) are the conservation of the
energy function $H(y)$ and the preservation of the symplecticity. In this
work, we intend to focus on more general equations than the canonical
Hamiltonian systems, where the matrix $J$ in (3) is allowed to be a constant
skew-symmetric matrix or a negative semidefinite matrix. A skew-symmetric
matrix $J$ grantees the conservation of the energy, while a negative
semidefinite matrix $J$ will lead to a dissipative system with Lyapunov
function $H(y)$ monotonically decreasing. In view of these special structures
in equation (3), a type of candidate methods will be the structure-preserving
exponential integrators which have superior qualitative behavior over long-
time integration compared with the general-purpose designed higher-order
methods [5]. Examples include symmetric methods [6], symplectic methods [7]
and energy-preserving methods [8]. Here we would like to consider particularly
energy-preserving exponential integrators. Such methods in previous work are
fully implicit, e.g., [8, 9, 10, 11], except for the recent work [12], where a
linearly implicit method for the nonlinear Klein–Gordon equation was
considered using the scalar auxiliary variable (SAV) [13] approach. For
linearly implicit methods, only one linear system is solved in each iteration
and thus less computational cost is needed. For this consideration, we expect
to construct linearly implicit methods in this work.

There are two mostly used techniques in creating linearly implicit energy-
preserving methods for general conservative/dissipative systems with gradient
flow, according to the authors’ best knowledge. The first one follows from
Furihata and co-authors, where multiple-point methods are used so that the
nonlinearity can be portioned out [14]. Further studies of this method were
presented in [15] and [16] using the concept of polarized energy and polarized
discrete gradient. The second technique is to combine the linearly implicit
Crank-Nicolson method and the invariant energy quadratization (IEQ) [17] or
the scalar auxiliary variable (SAV) [13] approach. Both the IEQ-based and the
SAV-based methods are applicable for nonlinear systems, including
nonpolynomials; however, these methods need bounded free energy regarding the
nonlinear terms. Besides, the linearly implicit methods constructed based on
the second approach have no symmetric property. This paper will focus on the
technique using polarized energy for problems with polynomial energy
functions. There are several reasons for choosing this technique. First, there
are huge amounts of PDEs and ODEs with polynomial energy functions, e.g., the
nonlinear Schrödinger equations, the nonlinear wave equations, the KdV
equations, the Camassa-Holm equations, the wind-induced oscillator, the
polynomial pendulum oscillator, and so on. Second, there are no other
restrictions for the nonlinear terms except for being polynomials. Third, the
general scheme of the methods based on the first technique looks much more
straightforward than the second technique. Last but not least, the linearly
implicit methods constructed based on the first technique can be symmetric,
and it has been shown that symmetric methods applied to (near-)integrable
reversible systems share similar properties to symplectic methods: linear
error growth, long-time near-conservation of first integrals and existence of
invariant tori [18]. Thus methods with symmetric property usually provide
prominent long-time behavior.

The paper is organized as follows. First, we construct the symmetric linearly
implicit energy-preserving exponential integrators and discuss their
properties in Section 2. In Section 3, numerical examples are presented to
illustrate the performance of the proposed method. In the last section, we
conclude the paper with a summary of the properties and advantages shared by
the method.

##  2 Symmetric linearly implicit energy-preserving exponential scheme

In this section, we combine the idea of constructing linearly implicit methods
using polarized energy [15] and the idea of constructing energy-preserving
exponential methods using discrete gradient [8] to build the symmetric
linearly implicit energy-preserving exponential integrators. To present the
method more intuitively, we restrict the nonlinear term $U(y)$ (potential
energy) in equation (3) to be a cubic polynomial. However, the method is also
applicable to problems with any higher-order polynomials, for which the
results will be introduced briefly in this paper too.

The critical point of using polarized energy to construct a linearly implicit
method is to portion out the nonlinearity over consecutive time steps. This
can be carried out by constructing quadratic polarized energy and then
performing the polarized discrete gradient method, similarly as shown in [19]
for a cubic polynomial. A systematical way of constructing a quadratic
polarization for higher-order polynomial functions is presented in [15], for
example

  * 1.

$U(x)=x^{2}$ can be polarized by
$\bar{U}(x,y)=\theta\frac{x^{2}+y^{2}}{2}+(1-\theta)xy$, $\theta\in[0,1]$,

  * 2.

$U(x)=x^{3}$ can be polarized by $\bar{U}(x,y)=x\frac{x+y}{2}y,$

  * 3.

$U(x)=x^{4}$ can be polarized by $\bar{U}(x,y)=x^{2}y^{2},$

  * 4.

$U(x)=x^{5}$ can be polarized by $\bar{U}(x,y,z,w)=xyzw\frac{x+y+z+w}{4}$,

  * 5.

$U(x)=x^{6}$ can be polarized by $\bar{U}(x,y,z)=x^{2}y^{2}z^{2}$.

Following from [19], $\overline{\nabla}\bar{U}$ is said to be a polarized
discrete gradient for a polarized energy $\bar{U}$ if the following conditions
hold

$\begin{split}\bar{U}(y,z)-\bar{U}(x,y)&=\frac{1}{2}(z-x)^{T}\overline{\nabla}\bar{U}(x,y,z),\\\
\overline{\nabla}\bar{U}(x,x,x)&=\nabla U(x).\end{split}$ (4)

We will use the polarized discrete gradient to construct the linearly implicit
energy-preserving exponential integrators. Consider the variation of constants
formula on interval $[t_{0},t_{0}+2h]$ for problem (1) in the form

$y(t_{0}+2h)=\text{exp}(2hJM)y_{0}+2h\int_{0}^{1}\text{exp}(2(1-\xi)hJM)J\nabla
U(y(t_{0}+2h\xi))d\xi.$ (5)

Substituting $\nabla U(y(t_{0}+2h\xi))$ with the polarized discrete gradient
$\overline{\nabla}\bar{U}(y_{0},y_{1},y_{2})$ in (5), we obtain the energy-
preserving exponential integrator for problem (3) as follows

$y_{n+2}=\text{exp}(2hJM)y_{n}+2h\phi(2hJM)J\overline{\nabla}\bar{U}(y_{n},y_{n+1},y_{n+2})\quad\quad\text{(LIEEP)}.$
(6)

######  Remark 1.

Suppose that a quadratic polarization of a higher-order polynomial function
$U(y_{n})$ has the form $\bar{U}(y_{n},\cdots,y_{n+i},\cdots y_{n+p-1})$. Then
the generalization of the polarized discrete gradient in (4) can be given by

$\begin{split}\bar{U}(y_{n+1},\cdots y_{n+p})-\bar{U}(y_{n},\cdots
y_{n+p-1})&=\frac{1}{p}(y_{n+p}-y_{n})^{T}\overline{\nabla}\bar{U}(y_{n},\cdots,y_{n+p}),\\\
\overline{\nabla}\bar{U}(y_{n},\cdots,y_{n+p})&=\nabla U(y_{n}).\end{split}$
(7)

######  Remark 2.

Suppose that the polarized discrete gradient for a higher-order polynomial
energy function $U(y_{n})$ is given by
$\overline{\nabla}\bar{U}(y_{n},y_{n+1},\cdots,y_{n+p})$ following from
equation (7). Then the linearly implicit energy-preserving exponential
integrator for problem (3) can be given by

$y_{n+p}=\text{exp}(phJM)y_{n}+ph\phi(phJM)J\overline{\nabla}\bar{U}(y_{n},y_{n+1},\cdots,y_{n+p}).$
(8)

######  Remark 3.

Scheme (6) is a special case of scheme (8) with $p=2$.

Before presenting the main theorems about the proposed method’s conservative
properties, we begin with a lemma.

######  Lemma 1.

For any symmetric matrix $M$, positive integer $p$ and scalar $h>0$, the
matrix

$B=\text{exp}(phJM)^{T}M\text{exp}(phJM)-M$

is zero when $J$ is skew-symmetric and negative semidefinite when $J$ is
negative semidefinite.

The proof of this lemma follows directly from Lemma 2.2 in [8] since the
result does not change when $h$ is replaced by $ph$.

######  Theorem 1.

Scheme (6) preserves the following polarized energy

$\bar{H}(y_{n},y_{n+1})=\frac{1}{4}(y_{n}^{T}My_{n}+y_{n+1}^{T}My_{n+1})+\bar{U}(y_{n},y_{n+1}).$
(9)

###### Proof.

We first assume that the matrix $M$ is not singular and denote by $V=2hJM$,
$M^{-1}\overline{\nabla}\bar{U}=\overline{\nabla}\tilde{U}$. The energy error
has the form

$\begin{split}&\bar{H}(y_{n+1},y_{n+2})-\bar{H}(y_{n},y_{n+1})\\\
&=\frac{1}{4}(y_{n+2}^{T}My_{n+2}+y_{n+1}^{T}My_{n+1})-\frac{1}{4}(y_{n}^{T}My_{n}+y_{n+1}^{T}My_{n+1})\\\
&+\bar{U}(y_{n+1},y_{n+2})-\bar{U}(y_{n},y_{n+1}).\\\ \end{split}$ (10)

By replacing
$y_{n+2}=\text{exp}(V)y_{n}+2h\phi(V)J\overline{\nabla}\bar{U}(y_{n},y_{n+1},y_{n+2})$
and using $\phi(V)V=\text{exp}(V)-I$, we get the following equations

$\begin{split}&\frac{1}{4}(y_{n+2}^{T}My_{n+2}+y_{n+1}^{T}My_{n+1})-\frac{1}{4}(y_{n}^{T}My_{n}+y_{n+1}^{T}My_{n+1})\\\
&=\frac{1}{4}y_{n}^{T}(\text{exp}(V)^{T}M\text{exp}(V)-M)y_{n}+hy_{n}^{T}\text{exp}(V)^{T}M\phi(V)J\overline{\nabla}\bar{U}\\\
&+h^{2}(\overline{\nabla}\bar{U})^{T}J^{T}\phi(V)^{T}M\phi(V)J\overline{\nabla}\bar{U}\\\
&=\frac{1}{4}y_{n}^{T}(\text{exp}(V)^{T}M\text{exp}(V)-M)y_{n}+\frac{1}{2}y_{n}^{T}\text{exp}(V)^{T}M(\text{exp}(V)-I)\overline{\nabla}\tilde{U}\\\
&+\frac{1}{4}(\overline{\nabla}\tilde{U})^{T}(\text{exp}(V)-I)^{T}M(\text{exp}(V)-I)\overline{\nabla}\tilde{U},\end{split}$
(11)

and

$\begin{split}&\bar{U}(y_{n+1},y_{n+2})-\bar{U}(y_{n},y_{n+1})\\\
&=\frac{(y_{n+2}-y_{n})^{T}}{2}\overline{\nabla}\bar{U}(y_{n},y_{n+1},y_{n+2})\\\
&=\frac{1}{2}y_{n}^{T}(\text{exp}(V)^{T}-I)\overline{\nabla}\bar{U}+h(\overline{\nabla}\bar{U})^{T}J^{T}\phi(V)^{T}\overline{\nabla}\bar{U}\\\
&=\frac{1}{2}y_{n}^{T}(\text{exp}(V)^{T}M-M)\overline{\nabla}\tilde{U}+\frac{1}{2}\overline{\nabla}\tilde{U}^{T}V^{T}\phi(V)^{T}M\overline{\nabla}\tilde{U}\\\
&=\frac{1}{2}y_{n}^{T}(\text{exp}(V)^{T}M-M)\overline{\nabla}\tilde{U}+\frac{1}{2}\overline{\nabla}\tilde{U}^{T}(\text{exp}(V)^{T}M-M)\overline{\nabla}\tilde{U}.\end{split}$
(12)

Inserting equations (11) and (12) to equation (10), we obtain the following
results

$\begin{split}&\bar{H}(y_{n+1},y_{n+2})-\bar{H}(y_{n},y_{n+1})\\\
&=\frac{1}{4}y_{n}^{T}(\text{exp}(V)^{T}M\text{exp}(V)-M)y_{n}+\frac{1}{2}y_{n}^{T}(\text{exp}(V)^{T}M\text{exp}(V)-M)\overline{\nabla}\tilde{U}\\\
&+\frac{1}{4}\overline{\nabla}\tilde{U}^{T}(\text{exp}(V)^{T}M\text{exp}(V)-M)\overline{\nabla}\tilde{U}+\frac{1}{4}\overline{\nabla}\tilde{U}^{T}(\text{exp}(V)^{T}M-M\text{exp}(V))\overline{\nabla}\tilde{U}\\\
&=\frac{1}{4}(y_{n}+\overline{\nabla}\tilde{U})^{T}(\text{exp}(V)^{T}M\text{exp}(V)-M)(y_{n}+\overline{\nabla}\tilde{U})\\\
&+\frac{1}{4}\overline{\nabla}\tilde{U}^{T}(\text{exp}(V)^{T}M-M\text{exp}(V))\overline{\nabla}\tilde{U}\\\
&=0,\end{split}$ (13)

where the last step follows from the fact that
$\text{exp}(V)^{T}M-M\text{exp}(V)$ is skew-symmetric, and
$\text{exp}(V)^{T}M\text{exp}(V)-M$ is also skew-symmetric from Lemma 1.

For a singular $M$, one can find a series of non-singular and symmetric
matrices $M_{\epsilon}$ such that $M_{\epsilon}=M$ when $\epsilon\rightarrow
0$. For any $M_{\epsilon}$, we can follow the proof above and show that the
polarized energy function in the form

$\bar{H_{\epsilon}}(y^{\epsilon}_{n},y^{\epsilon}_{n+1})=\frac{1}{4}({y^{\epsilon}_{n}}^{T}M_{\epsilon}y^{\epsilon}_{n}+{y^{\epsilon}_{n+1}}^{T}M_{\epsilon}y^{\epsilon}_{n+1})+\bar{U}(y^{\epsilon}_{n},y^{\epsilon}_{n+1})$

is preserved by the approximation given by

$y^{\epsilon}_{n+2}=\text{exp}(2hJM_{\epsilon})y_{n}^{\epsilon}+2h\phi(2hJM_{\epsilon})J\overline{\nabla}\bar{U}(y_{n}^{\epsilon},y_{n+1}^{\epsilon},y_{n+2}^{\epsilon})$

for the following problem

$\dot{y}(t)=JM_{\epsilon}y+J\nabla U(y(t)),\quad y(t_{0})=y_{0}.$

Therefore,
$\bar{H_{\epsilon}}(y^{\epsilon}_{n},y^{\epsilon}_{n+1})=\bar{H}(y_{n},y_{n+1})$
is preserved by method (6) when $\epsilon\rightarrow 0$.

∎

For problems (3) with higher-order polynomial energy functions, similar
conservation property can be obtained by scheme (8), as shown in the following
corollary.

######  Corollary 1.

Scheme (8) preserves the following polarized energy

$\bar{H}(y_{n},\cdots
y_{n+p-1})=\frac{1}{2p}\sum_{i=0}^{p-1}y_{n+i}^{T}My_{n+i}+\bar{U}(y_{n},\cdots,y_{n+p-1}).$
(14)

The proof is similar to Theorem 1 and thus is omitted here.

######  Theorem 2.

If $J$ is a constant negative semidefinite matrix, scheme (6) preserves the
polarized Lyapunov function $\bar{H}$ for problem (3):

$\bar{H}(y_{n+1},y_{n+2})\leq\bar{H}(y_{n+1},y_{n}),$

where $\bar{H}(y_{n},y_{n+1})$ is defined by equation (9).

###### Proof.

Let us suppose $M$ to be non-singular; otherwise we follow the similar
procedure in the proof for Theorem 1, i.e., constructing a series of
$M_{\epsilon}$ convergent to $M$ to achieve the result.

For a constant negative semidefinite matrix $J$, the error of Lyapnov function
has the same form as the energy error in (13):

$\begin{split}&\bar{H}(y_{n+1},y_{n+2})-\bar{H}(y_{n},y_{n+1})\\\
&=\frac{1}{4}(y_{n+2}^{T}My_{n+2}+y_{n+1}^{T}My_{n+1})-\frac{1}{4}(y_{n}^{T}My_{n}+y_{n+1}^{T}My_{n+1})\\\
&+\bar{U}(y_{n+1},y_{n+2})-\bar{U}(y_{n},y_{n+1})\\\
&=\frac{1}{4}(y_{n}+\overline{\nabla}\tilde{U})^{T}(\text{exp}(V)^{T}M\text{exp}(V)-M)(y_{n}+\overline{\nabla}\tilde{U})\\\
&+\frac{1}{4}\overline{\nabla}\tilde{U}^{T}(\text{exp}(V)^{T}M-M\text{exp}(V))\overline{\nabla}\tilde{U}\\\
&\leq 0,\end{split}$

where the last step follows from the fact that
$\text{exp}(V)^{T}M-M\text{exp}(V)$ is skew-symmetric, and
$\text{exp}(V)^{T}M\text{exp}(V)-M$ is negative semidefinite according to
Lemma 1.

∎

Similarly, for problems (3) with a higher-order polynomial energy, we have the
following corollary.

######  Corollary 2.

If $J$ is a constant negative semidefinite matrix, scheme (8) preserves the
polarized Lyapunov function:

$\bar{H}(y_{n+1},\cdots,y_{n+p})\leq\bar{H}(y_{n},\cdots,y_{n+p-1}),$

where $\bar{H}(y_{n},\cdots y_{n+p-1})$ has the same form as equation (14).

The proof is similar to Theorem 2 and thus is omitted here.

######  Theorem 3.

Scheme (6) is symmetric.

###### Proof.

Exchanging $y_{n}$, $y_{n+1}$ $\leftrightarrow$ $y_{n+2}$, $y_{n+1}$ and
replacing $2h$ by $-2h$ in (6), we obtain

$y_{n}=\text{exp}(-V)y_{n+2}-2h\phi(-V)J\overline{\nabla}\bar{U}(y_{n+2},y_{n+1},y_{n}).$

Following from the definition in (4) and the cyclic permutation free property
of the polarized energy [15], we get

$\overline{\nabla}\bar{U}(y_{n+2},y_{n+1},y_{n})=\overline{\nabla}\bar{U}(y_{n},y_{n+1},y_{n+2}).$
(15)

Using $\text{exp}(V)\phi(-V)=\phi(V)$ and the relation in (15), we then obtain

$y_{n+2}=\text{exp}(V)y_{n}+2h\phi(V)J\overline{\nabla}\bar{U}(y_{n},y_{n+1},y_{n+2}).$

∎

Scheme (8) for problems with higher-order polynomial $U(y)$ also holds the
symmetric property if the polarization of the function $U(y)$ is invariant
when the order of its arguments is reversed. This can be obtained by
symmetrizing over dihedral group [15]. Although only cyclic permutation free
is required in the definition of the polarized energy, we can actually always
get a permutation free quadratic polarization for any higher order polynomial
$\bar{U}(y)$. In fact, the polarization examples shown above are all
permutation free. In this paper, we always consider permutation free
polarization, i.e., quadratic polarization satisfying the following condition

$\bar{U}(y_{n},\cdots y_{n+p-1})=\bar{U}(y_{n+i_{1}-1},\cdots
y_{n+i_{p}-1}),\quad(i_{1},\cdots,i_{p})\in S_{p},$

where $S_{p}$ is a symmetric group.

######  Corollary 3.

Scheme (8) is symmetric if the polarization of function $U(y)$ is permutation
free.

The proof is similarly to the proof of Theorem 3 except that the cyclic
permutation free property should be replaced by the permutation free property.

##  3 Numerical experiment

The proposed method is suitable for conservative or dissipative differential
equations of the form (3) with $J$ a constant skew-symmetric or negative
semidefinite matrix and $U(y)$ a scalar polynomial function of any order.
These equations include the highly oscillatory conservative or dissipative
ODEs and also the semi-discrete systems arising from PDEs. In this section, we
test our method by three differential equations. The first two examples are
used to demonstrate the efficient behavior of the method compared to the fully
implicit method, e.g., the energy-preserving exponential integrators based on
the averaged vector field method, denoted by EAVF. The third example is chosen
to show the application of the method for problems with higher-order
polynomial energy functions.

EAVF method was put forward in [8], which has the form

$y_{n+1}=\text{exp}(hJM)y_{n}+h\phi(hJM)J\overline{\nabla}U(y_{n},y_{n+1}),\quad\quad$
(16)

where $\overline{\nabla}U(y_{n},y_{n+1})=\int_{0}^{1}\nabla
U((1-\tau)y_{n}+\tau y_{n+1})d\tau$. Besides, scheme (16) has been shown to
preserve the discrete energy of the form

$H(y_{n})=\frac{1}{2}y_{n}^{T}My_{n}+U(y_{n}).$ (17)

The integral in EAVF method is evaluated by the 2-point GL quadrature formula,
which gives an exact approximation of the integration. In most cases, the
terms $\text{exp}(phJM)$ and $\phi(phJM)$ ($p$ is a positive integer number)
can not be calculated explicitly, and we use the MATLAB package proposed in
[20] to compute them, where Pade approximations are used.

In the experiments, the global error is defined by

$\underset{n\geq 0}{\max}\|y_{n}-y(t_{n})\|,$

where $t_{n}=t_{0}+nh$ with $h$ the time step size, and $y(t_{n})$ is the
reference exact solution. In this work, we compute the reference solution by
the 6-order continuous Runge–Kutta (CRK) method [21] with the form

$\begin{cases}y_{n+1/3}=y_{n}+hJ\int_{0}^{1}\big{(}\frac{37}{27}-\frac{32}{9}\sigma+\frac{20}{9}\sigma^{2}\big{)}\nabla
H(Y_{\sigma})d\sigma\\\
\vspace{5pt}y_{n+2/3}=y_{n}+hJ\int_{0}^{1}\big{(}\frac{26}{27}+\frac{8}{9}\sigma-\frac{20}{9}\sigma^{2}\big{)}\nabla
H(Y_{\sigma})d\sigma\\\ y_{n+1}=y_{n}+hJ\int_{0}^{1}\nabla
H(Y_{\sigma})d\sigma\end{cases},$

where

$\begin{split}Y_{\sigma}=&-\frac{(3\sigma-1)(3\sigma-2)(\sigma-1)}{2}y_{n}+\frac{3\sigma(3\sigma-2)(3\sigma-3)}{2}y_{n+1/3}\\\
&-\frac{3\sigma(3\sigma-1)(3\sigma-3)}{2}y_{n+2/3}+\frac{\sigma(3\sigma-1)(3\sigma-2)}{2}y_{n+1},\end{split}$

and the integrals are evaluated exactly by the 5-point GL quadrature. For all
fully implicit schemes, we solve the nonlinear system by the fixed point
iteration with tolerance as $10^{-14}$. All the numerical results presented
are obtained from schemes implemented in MATLAB (2020a release), running on a
MacBook Pro with a dual-core 2.6 GHz Intel 6-Core i7 processor and 16 GB of
2667 MHz DDR4 RAM.

Test problem one. We consider an averaged system in wind-induced oscillation
[22]

$\begin{split}\dot{x}_{1}&=-\zeta x_{1}-\lambda x_{2}+x_{1}x_{2},\\\
\dot{x}_{2}&=\lambda x_{1}-\zeta
x_{2}+\frac{1}{2}(x_{1}^{2}-x_{2}^{2}),\end{split}$ (18)

where $\zeta\geq 0$ is a damping factor and $\lambda$ is a detuning parameter
with $\zeta=r\text{cos}(\theta)$, $\lambda=r\text{sin}(\theta)$, $r\geq 0$,
$0\leq\theta\leq\pi/2$. Equation (18) can be rewritten into the form (3) with

$J=\begin{bmatrix}-cos(\theta)&-sin(\theta)\\\
sin(\theta)&-cos(\theta)\end{bmatrix}\quad M=\begin{bmatrix}r&0\\\
0&r\end{bmatrix},$
$U=-\frac{1}{2}sin(\theta)(x_{1}x_{2}^{2}-\frac{1}{3}x_{1}^{3})+\frac{1}{2}cos(\theta)(\frac{1}{3}x_{2}^{3}-x_{1}^{2}x_{2}).$

Its energy function (when $\theta=\pi/2$) or Lyapunov function (dissipative
case, when $\theta\leq\pi/2$) is

$H=\frac{1}{2}r(x_{1}^{2}+x_{2}^{2})-\frac{1}{2}sin(\theta)(x_{1}x_{2}^{2}-\frac{1}{3}x_{1}^{3})+\frac{1}{2}cos(\theta)(\frac{1}{3}x_{2}^{3}-x_{1}^{2}x_{2}).$

The matrix exponential in scheme (6) for problem (18) can be calculated
explicitly as follows

$exp(V)=\begin{bmatrix}exp(-2hcr)cos(2hsr)&-exp(-2hcr)sin(2hsr)\\\
exp(-2hcr)sin(2hsr)&exp(-2hcr)cos(2hsr)\end{bmatrix},$

with $c=cos(\theta)$ and $s=sin(\theta)$. We can obtain a polarized discrete
gradient $\nabla\bar{U}(x^{n},x^{n+1},x^{n+2})$ based on a polarization of $U$
given by

$\begin{split}\bar{U}(x^{n},x^{n+1})=&-\frac{1}{2}sin(\theta)\big{(}a\frac{x_{1}^{n}+x_{1}^{n+1}}{2}x_{2}^{n}x_{2}^{n+1}+(1-a)\frac{x_{1}^{n}(x_{2}^{n+1})^{2}+x_{1}^{n+1}(x_{2}^{n})^{2}}{2}\\\
&-\frac{1}{3}x_{1}^{n}\frac{x_{1}^{n}+x_{1}^{n+1}}{2}x_{1}^{n+1}\big{)}+\frac{1}{2}cos(\theta)\big{(}\frac{1}{3}x_{2}^{n}\frac{x_{2}^{n}+x_{2}^{n+1}}{2}x_{2}^{n+1}\\\
&-ax_{1}^{n}x_{1}^{n+1}\frac{x_{2}^{n}+x_{2}^{n+1}}{2}-(1-a)\frac{x_{2}^{n}(x_{1}^{n+1})^{2}+x_{2}^{n+1}(x_{1}^{n})^{2}}{2}\big{)}.\end{split}$
(19)

Then we get the linearly implicit energy-preserving scheme in the form of (6)
and the polarized energy in the form of (9).

Consider the initial vector $x_{1}(0)=0,\quad x_{2}(0)=1$, $r=20$, step size
$h=1/20$, and the parameters $\theta=\pi/2$ or $\theta=\pi/2-10^{-4}$. For
LIEEP method, the starting point $[x_{1}^{1},x_{2}^{1}]$ is computed by the
Matlab function _ode45_. $\theta=\pi/2$ provides a conservative system, and
Figure 1(a) confirms that EAVF method preserves the discrete energy (17) and
LIEEP method preserves the polarized energy (9). While $\theta=\pi/2-10^{-4}$
leads to a dissipative system, and Figure 1(b) shows that EAVF method and
LIEEP method preserves the dissipation of the Lyapunov function in (17) and
(9), respectively.

(a) $\theta=\pi/2$

(b) $\theta=\pi/2-10^{-4}$

Figure 1:  The energy behaviour of both EAVF and LIEEP method with time step
size $h=1/20$.

In Figure 2 and 3, we consider the global errors and the computational cost
using step sizes $h=1/10/2^{i}$, with $i=0,1,\cdots,5$. Surprisingly, Figure
2(a) shows that LIEEP method is superconvergent for the conservative system
($\theta=\pi/2$). We find that this behavior is closely related to the
parameter $a$ in the polarized potential energy in (19). We have tried
$a=0,1/2,1/4,1$, but only $a=1/2$ gives a three-order behaviour. Figure 2(b)
shows that the proposed method is more efficient than the fully implicit EAVF
method. When $\theta=\pi/2-10^{-4}$, i.e., the system is dissipative, the
superconvergent behavior disappears for LIEEP method, see Figure 3(a). From
this figure, we observe that LIEEP method has a convergent issue when the step
size is $h=1/10$, but with the decrease of the time step size, LIEEP method
gets convergent and behaves even better than EAVF method. Figure 3(b)
indicates that the proposed method is much more efficient than EAVF method for
the dissipative system.

(a) order plot

(b) efficiency plot

Figure 2: $T=1000$, $\theta=\pi/2$, the time step sizes are $h=1/10/2^{i}$,
for $i=0,1,\cdots,5$.

(a) order plot

(b) efficiency plot

Figure 3: $T=1000$, $\theta=\pi/2-10^{-4}$, the time step sizes are
$h=1/10/2^{i}$, for $i=0,1,\cdots,5$.

Test problem two. We consider a continuous generalization of an $\alpha$-FPU
(Fermi-Pasta-Ulam) system [23]:

$\frac{\partial^{2}u}{\partial t^{2}}=\beta\frac{\partial^{3}u}{\partial
t\partial x^{2}}+\frac{\partial^{2}u}{\partial
x^{2}}(1+\epsilon(\frac{\partial u}{\partial x})^{p})-\gamma\frac{\partial
u}{\partial t}-m^{2}u,$ (20)

where $\epsilon>0$, $\beta\geq 0$ is the coefficient of the internal damping,
$\gamma\geq 0$ is the coefficient of the external damping, and
$(x,t)\in[0,L]\times[0,T]$. Taking $\frac{\partial u}{\partial t}=v$, equation
(20) can be rewritten as

$\begin{split}\partial_{t}u&=v\\\
\partial_{t}v&=\beta\frac{\partial^{2}v}{\partial
x^{2}}+\frac{\partial^{2}u}{\partial x^{2}}(1+\epsilon(\frac{\partial
u}{\partial x})^{p})-\gamma v-m^{2}u.\end{split}$ (21)

Denoting by $y=[u,v]^{T}$, equation (21) can be reformulated as the following
Hamiltonian form

$\frac{\partial y}{\partial t}=\mathcal{Q}\frac{\delta\mathcal{H}}{\delta y},$

where

$\mathcal{Q}=\begin{bmatrix}0&1\\\
-1&\beta\partial^{2}_{x}-\gamma\end{bmatrix},\quad\mathcal{H}=\int_{0}^{L}E(t,u,v,u_{x})dx,$

with

$E(t,u,v,u_{x})=\frac{1}{2}u_{x}^{2}+\frac{m^{2}}{2}u^{2}+\frac{1}{2}v^{2}+\epsilon\frac{u_{x}^{p+2}}{(p+2)(p+1)}.$

The function $E(t,u,v,u_{x})$ physically represents the local energy density
of system (20) at any time $t$.

Consider $p=1$ and the homogeneous Dirichlet boundary conditions
$u(0,t)=u(L,t)=0$. Discretizing $\partial_{x}^{2}$ with the central difference
operator and $\partial_{x}$ with the forward difference operator, we obtain
the following semi-discrete ODE system

$\dot{y}=Q(My+\nabla U(y)),$

where

$Q=\begin{bmatrix}0&I\\\ -I&\beta D-\gamma I\end{bmatrix},\quad
M=\begin{bmatrix}m^{2}I-D&0\\\ 0&I\end{bmatrix},\quad
U(y)=\sum_{j=0}^{N-1}\frac{\epsilon}{6}(\frac{u_{j+1}-u_{j}}{\Delta x})^{3}.$

Setting $w_{j}^{n}=\frac{u_{j+1}^{n}-u_{j}^{n}}{\Delta x}$, and defining the
polarized energy

$\bar{U}(w_{j}^{n},w_{j}^{n+1})=\sum_{j=0}^{N-1}\frac{\epsilon}{6}w_{j}^{n}\frac{w_{j}^{n}+w_{j}^{n+1}}{2}w_{j}^{n+1},$

we can obtain the polarized discrete gradient

$\begin{split}\bar{\nabla}\bar{U}(w_{j}^{n},w_{j}^{n+1},w_{j}^{n+2})&=\frac{\epsilon}{6\Delta
x}w_{j-1}^{n+1}(w_{j-1}^{n}+w_{j-1}^{n+1}+w_{j-1}^{n+2})\\\
&-\frac{\epsilon}{6\Delta
x}w_{j}^{n+1}(w_{j}^{n}+w_{j}^{n+1}+w_{j}^{n+2}),\end{split}$

and the discrete gradient

$\bar{\nabla}U(w_{j}^{n})=\frac{\epsilon}{2\Delta
x}(w_{j-1}^{n})^{2}-\frac{\epsilon}{2\Delta x}(w_{j}^{n})^{2}.$

We consider $m=0$, $\epsilon=\frac{3}{4}$, $L=128,\quad T=100$ and spatial
step size $\Delta x=1$. The initial conditions are set to be
$u_{j}(0)=q_{j}(0)$, $v_{j}(0)=\dot{q}_{j}(0)$ and

$q_{j}(t)=5\text{ln}\frac{1+\text{exp}\big{(}2(\alpha(j-97)+t\text{sinh}(\alpha))\big{)}}{1+\text{exp}\big{(}2(\alpha(j-96)+t\text{sinh}(\alpha))\big{)}}+5\text{ln}\frac{1+\text{exp}\big{(}2(\alpha(j-32)+t\text{sinh}(\alpha))\big{)}}{1+\text{exp}\big{(}2(\alpha(j-33)+t\text{sinh}(\alpha))\big{)}},$

where $\alpha=0.1$. For LIEEP method, the starting point $y_{1}$ is computed
by the 6-order CRK method.

In Figure 4(a), we fix the external damping coefficient to be zero
($\gamma=0$) and present the energy behavior of LIEEP method for systems with
different internal damping coefficients and a long simulation time $T=500$. We
observe that the numerical method preserves the energy when there is no
damping ($\beta=0$) and also preserves the dissipation property when the
internal damping coefficient is greater than zero, consistent with what is
observed in [23], where a fully implicit four-step method is considered.
Similar behavior is observed in Figure 4(b), where the internal damping
coefficient is set to be zero ($\beta=0$). Figure 5 and 6 confirm that both
EAVF and LIEEP method are of order 2 in time, and the comparison of the
computational cost between these two methods gives a clear evidence that the
proposed method is more efficient than EAVF method. In this experiment, we
also present the numerical solutions given by LIEEP method for systems with
different settings of $\beta$ and $\gamma$, see Figure 7. These figures
clearly demonstrate the dissipative nature of the external damping
coefficient, see the change of the colors between Figure 7(a) and Figure 7(b),
and the internal damping coefficient, see the change of the shapes between
Figure 7(a) and Figure 7(c). These observations in the numerical solutions are
in accordance with the results shown by the fully implicit four-step method in
[23].

(a) $\gamma=0$

(b) $\beta=0$.

Figure 4: Polarized energy by LIEEP method for $\alpha$-FPU system with
different settings of internal and external damping coefficients. $T=500$ and
time step size $h=0.025$.

Figure 5:  $T=100$, $\gamma=0.005$, $\beta=0$ and time step size
$h=\frac{1}{2^{i}}$, $i=1,\cdots,5$, Left: efficiency comparison; right: order
plot.

Figure 6: $T=100$, $\gamma=0$, $\beta=2$ and time step size
$h=\frac{1}{2^{i}}$, $i=1,\cdots,5$. Left: efficiency comparison; right: order
plot.

(a) $\gamma=0$, $\beta=0$.

(b) $\gamma=0.005$, $\beta=0$.

(c) $\gamma=0$, $\beta=2$.

Figure 7:  The numerical solution of the $\alpha$-FPU system with different
settings of internal and external damping coefficients on time interval
$[0,500]$ and with time step size $h=0.025$.

Test problem three. We consider the polynomial pendulum oscillator, and the
main focus of this example is to illustrate the energy conservation property
shown in Corollary 1 for the proposed method in Remark 2. Consider the
nonlinear pendulum problem with the Hamiltonian

$H(p,q)=\frac{1}{2}p^{2}+1-\text{cos}q,$

and a truncated Taylor expansion of the cosine function:

$H(p,q)=\frac{1}{2}p^{2}+\frac{1}{2}q^{2}-\frac{1}{24}q^{4}+\frac{1}{720}q^{6}.$
(22)

The approximation in (22) to the original problem will be more accurate if
even higher-order polynomial is used and $\lvert q\rvert$ is sufficiently
small, e.g., $\lvert q\rvert<\frac{1}{2}$ [24]. Denoting by $y=[q;p]$, the
polynomial pendulum oscillator with energy function (22) can be rewritten into
form (3) with $m=2$, $J$ the canonical skew-symmetric matrix, M the identity
matrix and

$U(y)=-\frac{1}{24}q^{4}+\frac{1}{720}q^{6}.$ (23)

Consider a polarization of (23) as follows

$\bar{U}(y_{n},y_{n+1},y_{n+2})=-\frac{1}{24}q_{n}q_{n+1}q_{n+2}\frac{q_{n}+q_{n+1}+q_{n+2}}{3}+\frac{1}{720}{q_{n}}^{2}{q_{n+1}}^{2}{q_{n+2}}^{2}.$
(24)

We can obtain a polarized discrete gradient of the form

$\begin{split}\bar{\nabla}\bar{U}(y_{n},y_{n+1},y_{n+2},y_{n+3})=&\frac{1}{240}{q_{n+1}}^{2}{q_{n+2}}^{2}(q_{n}+q_{n+3})\\\
&-\frac{1}{24}q_{n+1}q_{n+2}(q_{n}+q_{n+1}+q_{n+2}+q_{n+3}).\end{split}$

Take the initial value as $q_{0}=0.5$, $p_{0}=1$ and the integration interval
as $[0,1000]$. We compute the first two starting points $y_{1}$ and $y_{2}$ by
Matlab function _ode15s_. The polarized energy is reported in Figure 8(a), and
we observe that it is exactly preserved by LIEEP method defined by equation
(8). In this figure, we also present the original discrete energy by LIEEP
method , i.e.,

$H(p_{n},q_{n})=\frac{1}{2}p_{n}^{2}+1-\text{cos}q_{n}.$ (25)

Although LIEEP method does not preserve the exact original energy, Figure 8(a)
shows that the discrete energy given by LIEEP method in the form of (25) stays
oscillated and bounded over a long-time integration. Besides, we observe that
the numerical solution by LIEEP method applied to the truncated equation
provides an approximation with a similar behaviour as the exact solution of
the nonlinear pendulum oscillator if a small time step size is considered,
e.g., h=0.3, i.e., the phase space is a cylinder, as illustrated in Figure
8(b).

(a) Energy

(b) solution

Figure 8:  In 8(a), the time step size is $h=1$, the polarized energy is
defined by equation (14) with $p=3$ and $\bar{U}(y_{n},y_{n+1},y_{n+2})$
defined by (24); the discrete energy is defined by equation (25). In 8(b), the
time step size is $h=0.3$.

##  4 Conclusion

This paper constructs a novel symmetric linearly implicit exponential
integrator that holds the conservative properties for semi-linear problems
with polynomial energy functions. The method is developed based on combining
the idea of using polarized discrete gradient to build linearly implicit
methods and the idea of using discrete gradient to create energy-preserving
exponential integrators. Besides conservative properties, the method is shown
to be symmetric, which guarantees excellent long-time behavior.

We test our methods on three types of differential equations, including an
oscillated ODE, i.e., the averaged wind-induced oscillator, an oscillated PDE,
i.e., the damped FPU problem, and also an ODE with higher-order polynomial
energy fucntion, i.e., the polynomial pendulum oscillator. The numerical
experiments confirm that the proposed method preserves the polarized energy or
the Lyapunov function, and the method is of order two. Moreover, it has been
shown that the proposed method has superconvergent behavior for some
particular systems when a proper polarized energy is considered. Compared with
the fully implicit method (EVAF), our method shows a significantly lower
computational cost. In view of the nice properties and the good behavior, we
recommend the proposed method for problems with polynomial energy function.

## Acknowledgement

The author would like to thank Isaac Newton Institute for Mathematical
Sciences, Cambridge, for support and hospitality during the programme
_Geometry, compatibility and structure preservation in computational
differential equations (2019)_ under Grant number EP/R014604/1, where work on
this paper was partly carried out.

The author would also like to thank the European Union Horizon 2020 research
and innovation programme under the Marie Skłodowska-Curie grant agreement No.
691070 CHiPS.

## References

  * Cary and Brizard [2009] J. R. Cary, A. J. Brizard,  Hamiltonian theory of guiding-center motion,  Rev. Modern Phys. 81 (2009) 693–738. URL: https://doi.org/10.1103/RevModPhys.81.693. doi:10.1103/RevModPhys.81.693.
  * Cotter and Reich [2006] C. J. Cotter, S. Reich,  Semigeostrophic particle motion and exponentially accurate normal forms,  Multiscale Model. Simul. 5 (2006) 476–496. URL: https://doi.org/10.1137/05064326X. doi:10.1137/05064326X.
  * Hochbruck and Ostermann [2010] M. Hochbruck, A. Ostermann,  Exponential integrators,  Acta Numer. 19 (2010) 209–286. URL: https://doi.org/10.1017/S0962492910000048. doi:10.1017/S0962492910000048.
  * Hochbruck et al. [0809] M. Hochbruck, A. Ostermann, J. Schweitzer,  Exponential Rosenbrock-type methods,  SIAM J. Numer. Anal. 47 (2008/09) 786–803. URL: https://doi.org/10.1137/080717717. doi:10.1137/080717717.
  * Hairer et al. [2006] E. Hairer, C. Lubich, G. Wanner, Geometric numerical integration, volume 31 of Springer Series in Computational Mathematics, second ed., Springer-Verlag, Berlin, 2006\. Structure-preserving algorithms for ordinary differential equations.
  * Celledoni et al. [2008] E. Celledoni, D. Cohen, B. Owren,  Symmetric exponential integrators with an application to the cubic Schrödinger equation,  Found. Comput. Math. 8 (2008) 303–317. URL: https://doi.org/10.1007/s10208-007-9016-7. doi:10.1007/s10208-007-9016-7.
  * Wu and Wang [2018] X. Wu, B. Wang, Recent developments in structure-preserving algorithms for oscillatory differential equations, Science Press Beijing, Beijing; Springer, Singapore, 2018. URL: https://doi.org/10.1007/978-981-10-9004-2. doi:10.1007/978-981-10-9004-2.
  * Li and Wu [2016] Y.-W. Li, X. Wu,  Exponential integrators preserving first integrals or Lyapunov functions for conservative or dissipative systems,  SIAM J. Sci. Comput. 38 (2016) A1876–A1895. URL: https://doi.org/10.1137/15M1023257. doi:10.1137/15M1023257.
  * Miyatake [2014] Y. Miyatake,  An energy-preserving exponentially-fitted continuous stage Runge–Kutta method for Hamiltonian systems,  BIT Numerical Mathematics 54 (2014) 777–799.
  * Cui et al. [2021] J. Cui, Z. Xu, Y. Wang, C. Jiang,  Mass-and energy-preserving exponential Runge–Kutta methods for the nonlinear Schrödinger equation,  Applied Mathematics Letters 112 (2021) 106770.
  * Shen and Leok [2019] X. Shen, M. Leok,  Geometric exponential integrators,  Journal of Computational Physics 382 (2019) 27–42.
  * Jiang et al. [2020] C. Jiang, Y. Wang, W. Cai,  A linearly implicit energy-preserving exponential integrator for the nonlinear Klein-Gordon equation,  Journal of Computational Physics 419 (2020) 109690.
  * Shen et al. [2019] J. Shen, J. Xu, J. Yang,  A new class of efficient and robust energy stable schemes for gradient flows,  SIAM Review 61 (2019) 474–506.
  * Furihata and Matsuo [2011] D. Furihata, T. Matsuo, Discrete variational derivative method, Chapman & Hall/CRC Numerical Analysis and Scientific Computing, CRC Press, Boca Raton, FL, 2011. A structure-preserving numerical method for partial differential equations.
  * Dahlby and Owren [2011] M. Dahlby, B. Owren,  A general framework for deriving integral preserving numerical methods for PDEs,  SIAM J. Sci. Comput. 33 (2011) 2318–2340. URL: https://doi.org/10.1137/100810174. doi:10.1137/100810174.
  * Eidnes and Li [2020] S. Eidnes, L. Li,  Linearly implicit local and global energy-preserving methods for PDEs with a cubic Hamiltonian,  SIAM Journal on Scientific Computing 42 (2020) A2865–A2888.
  * Zhao et al. [2017] J. Zhao, Q. Wang, X. Yang,  Numerical approximations for a phase field dendritic crystal growth model based on the invariant energy quadratization approach,  International Journal for Numerical Methods in Engineering 110 (2017) 279–300.
  * Hairer et al. [2006] E. Hairer, C. Lubich, G. Wanner, Geometric numerical integration: structure-preserving algorithms for ordinary differential equations, volume 31, Springer Science & Business Media, 2006.
  * Eidnes et al. [2019] S. Eidnes, L. Li, S. Sato,  Linearly implicit structure-preserving schemes for Hamiltonian systems,  Journal of Computational and Applied Mathematics (2019) 112489.
  * Berland et al. [2007] H. Berland, B. Skaflestad, W. M. Wright,  Expint—a Matlab package for exponential integrators,  ACM Transactions on Mathematical Software (TOMS) 33 (2007) 4–es.
  * Hairer [2010] E. Hairer,  Energy-preserving variant of collocation methods,  Journal of Numerical Analysis, Industrial and Applied Mathematics 5 (2010) 73–84.
  * McLachlan et al. [1998] R. I. McLachlan, G. Quispel, N. Robidoux,  Unified approach to Hamiltonian systems, Poisson systems, gradient systems, and systems with Lyapunov functions or first integrals,  Physical Review Letters 81 (1998) 2399.
  * Macías-Díaz and Medina-Ramírez [2009] J. Macías-Díaz, I. Medina-Ramírez,  An implicit four-step computational method in the study on the effects of damping in a modified $\alpha$-Fermi–Pasta–Ulam medium,  Communications in Nonlinear Science and Numerical Simulation 14 (2009) 3200–3212.
  * Iavernaro and Trigiante [2009] F. Iavernaro, D. Trigiante,  High-order symmetric schemes for the energy conservation of polynomial Hamiltonian problems,  J. Numer. Anal. Ind. Appl. Math 4 (2009) 87–101.