# Beating resonance patterns and extreme power flux skewing in anisotropic
elastic plates

Daniel A. Kiefer,∗ Sylvain Mezil, Claire Prada

Institut Langevin, ESPCI Paris, Université PSL, CNRS, 75005 Paris, France
∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
###### Abstract

Elastic waves in anisotropic media can exhibit a power flux that is not
collinear with the wave vector. This has notable consequences for waves guided
in a plate. Through laser-ultrasonic experiments, we evidence remarkable
phenomena due to slow waves in a single crystal silicon wafer. Waves
exhibiting power flux orthogonal to their wave vector are identified. A pulsed
line source that excites these waves reveals a wave packet radiated parallel
to the line. Furthermore, there exist precisely eight plane waves with zero
power flux. These so-called zero-group-velocity modes are oriented along the
crystal’s principal axes. Time acts as a filter in the wave vector domain that
selects these modes. Thus, a point source leads to beating resonance patterns
with moving nodal curves on the surface of the infinite plate. We observe this
pattern as it emerges naturally after a pulsed excitation.
Published version: Science Advances 9 (51), p. eadk6846, doi
10.1126/sciadv.adk6846.

##  1 Introduction

Engineering materials often exhibit anisotropic stiffness, in particular,
single crystals and composites. Among them, monocrystalline silicon is the
single most important material for fabrication of integrated circuits,
microelectromechanical systems (MEMS) and photovoltaic cells [1, 2]. These are
generally produced from a thin wafer of this material. Understanding the
intricate mechanics of elastic wave propagation in these plates is not only of
importance for their quality evaluation [3, 4, 5], but also for the functional
design of MEMS devices such as surface acoustic wave (SAW) and Lamb wave
filters and sensors [6, 7, 8]. The latter usually involve layers of
piezoelectric materials, which constitute another kind of medium where
anisotropy plays a major role [9].

As the structures in the mentioned applications are usually thin, guided
elastic waves are of great relevance. These waves propagate dispersively [10,
11, 12], i.e., their angular frequency $\omega$ is non-linearly related to
their wavenumber $k$ through the dispersion relation $\omega(k)$. Of
particular interest are solutions where the group velocity component
${\partial\omega/\partial k}$ vanishes while the wavenumber remains finite
[13, 14, 15, 16, 17, 18]. In isotropic media, these zero-group-velocity (ZGV)
points represent local resonances that are due to the finite thickness of the
semi-infinite structure. At sufficiently high frequencies they usually
dominate the overall mechanical response of the structure. They are rather
simple to excite and measure with contactless laser-ultrasonic techniques. ZGV
resonances are used in nondestructive evaluation to determine various
structural properties such as material parameters, thickness or bonding state
with very high precision in a spatially resolved manner [19, 20, 21, 22, 23].

Figure 1: Angular dependence of the ZGV resonance. The surface displacement
was acquired at the center of a line source and obtained on two different time
windows. The resonance frequency clearly depends on the line orientation with
a 90∘-periodicity. $\theta$ is the angle of the excited wave vector and is
counted from the [110] axis. For long times, the resonance remains at
$\theta=n\times 45^{\circ}$, $n\in\mathbb{Z}$, but vanishes elsewhere. The
measurement data was taken from the experiment detailed in Ref. [24], where
(A) has already been reported. Figure 2: Sample and measurement setup. (A)
[001]-cut silicon wafer sample. The
${\mathbf{e}_{X}\mathbf{e}_{Y}\mathbf{e}_{Z}}$ coordinate system is fixed as
indicated, while the ${\mathbf{e}_{r}\mathbf{e}_{\theta}\mathbf{e}_{z}}$
system is aligned with the wave propagation direction. The wave vector
$\mathbf{k}$ is inclined at angle $\theta$, while the group velocity vector
$\mathbf{c}_{\mathup{g}}$ is at angle $\theta+\alpha$. (B) schematic
representation of the laser-ultrasonic measurement setup. BE: Beam Expander,
L: Lense.

ZGV resonances in anisotropic plates were the object of some theoretical [25,
26, 27] as well as experimental studies [24, 28]. Prada et al. [24] revealed
the directional dependence of the ZGV resonance frequency in a silicon wafer
and we reproduce this result in Fig. 1A (we will show further on that not all
of the components really correspond to ZGV resonances). Note that the
directional dependence is distinctive of ZGV resonances as this is impossible
for common thickness resonances characterized by ${k=0}$. It is striking that
the extremal frequencies in Fig. 1 play a special role, in particular, when
inspecting later times as depicted in Fig. 1B. Moreover, in Ref. [24] the
response to the point-source excitation appears to be dominated by the two ZGV
frequencies associated to principal directions of the crystal. This fact was
confirmed by numerical computations in Ref. [26]. However, the described
effects remain unexplained in the literature.

The aim of the present work is twofold. Firstly, the special role of ZGV
resonances along principal directions of the material is explained. The
contributions along other directions are explained by slow waves that
propagate parallel to the line source, which means that they exhibit a group-
velocity vector that is orthogonal to the wave vector. Secondly, we predict
and confirm the formation of complex time-dependent resonance patterns on the
surface of the plate after it has been excited by an impulsive point load. We
show that this beating pattern can be explained by the interference of eight
ZGV modes associated to the material’s principal directions.

##  2 Results

We study transversely propagating waves and ZGV resonances in anisotropic
plates. Both effects are confirmed experimentally on a [001]-cut
monocrystalline silicon wafer of 524.6 µm thickness, see Fig. 2A. The
material’s stiffness is of cubic anisotropy (Voigt-notated stiffness
$C_{11}=\mathup{165.6\,GPa}$, $C_{12}=\mathup{63.9\,GPa}$,
$C_{44}=\mathup{79.5\,GPa}$, density $\rho=\mathup{2330\,kg/m^{3}}$). The
laser-ultrasonic system depicted in Fig. 2B is used to observe the waves. It
consists of a pulsed laser source that is defocused to excite the wavenumber
range of interest of up to $\approx$7 rad/mm and of a laser interferometer
that measures the normal surface displacements of the wafer. A two-dimensional
scan of the surface is performed by moving the laser source. To speed up the
measurement, we only scan one quarter of the wave field and reconstruct the
full field by exploiting the sample’s material symmetries. Some of the results
are bandpass filtered between 7.6 MHz to 7.8 MHz in a postprocessing step
(indicated explicitly in the figures), as this frequency band contains the
phenomena of interest. Further specificities on the measurement setup and the
postprocessing can be found in Sec. 4.2. Before presenting the experimental
results, the underlying theory is developed in the following.

Figure 3: Dispersion surface of the S1/S2b mode close to the minima in
frequency. (A) iso-frequency contours in the wave vector plane and (B) surface
in the frequency-wave vector space. ZGV: zero-group velocity. TGV: transverse-
group velocity. The ZGV points are located at the marked minima and saddle-
points (red points). For visualization purposes, (B) has been rotated by 22.5∘
and cropped to the shown quarter plane. Figure 4: Dispersion of wave vectors
and group velocity. (A) frequency vs. wavenumbers. (B) some group velocity
vectors sketched on a selected iso-frequency contour. Dispersion curves are
cuts along constant $\theta$ of this surface. (C) components $c_{\mathup{g}r}$
and $c_{\mathup{g}\theta}$ of the group velocity vector vs. frequency, compare
also to (B). The transverse power flux due to $c_{\mathup{g}\theta}$ does not
vanish for TGV waves, while it does vanish at ZGV points. (D) magnitude of the
group velocity vector vs. frequency. (E) skew angles $\alpha$ between the
group velocity vectors $\mathbf{c}{g}$ and the wave vectors $\mathbf{k}$
plotted vs. frequency.

###  2.1 Guided waves in a silicon wafer

We study guided elastic waves propagating in the wafer. The fixed
$\mathbf{e}_{X}\mathbf{e}_{Y}\mathbf{e}_{Z}$ coordinate system is aligned with
the [110],[$\bar{1}$10] and [001] crystallographic axes, respectively.
Moreover, the $\mathbf{e}_{r}\mathbf{e}_{\theta}\mathbf{e}_{z}$ system is
oriented with the wave vector, as depicted in Fig. 2A. Taking the point of
view of a plane wave, $\mathbf{e}_{r}$ denotes the _axial direction_ , while
$\mathbf{e}_{\theta}$ refers to the _transverse direction_. The plate has
stress-free surfaces. Furthermore, it is considered to be of infinite lateral
dimensions, so that waves reflected from the border can be disregarded.

Guided waves are characterized by the angular frequency $\omega$, the wave
vector
${\mathbf{k}=k\mathbf{e}_{r}(\theta)}=k_{X}\mathbf{e}_{X}+k_{Y}\mathbf{e}_{Y}$,
as well as the through-thickness displacement distribution $\mathbf{u}(z)$.
Only certain combinations of frequency and wave vector can propagate, which is
described by the _dispersion relation_ [10, 11, 12, 29]. The solutions form
surfaces $\omega(k_{X},k_{Y})$, or equivalently $\omega(k,\theta)$. We use a
semi-analytical method to obtain these solutions, i.e., the associated
eigenvalue problem is solved numerically. Our implementation is made available
as GEWtool [30]. Therewith, we are also able to compute ZGV points (and the
transversely propagating waves that will be discussed later on) by employing
the numerical methods that we have recently developed for this purpose [25,
31]. For details on the model and the numerical methods, see Section 4.1.

The computed dispersion surface of the first modes that exhibit ZGV points is
shown in Fig. 3. Note that due to the cubic anisotropy of silicon, the
dispersion surface has an angular periodicity of 90∘. Furthermore, we observe
reflection symmetry across the directions $\theta=\leavevmode\nobreak\ $0∘ and
$\theta=\leavevmode\nobreak\ $45∘, which correspond to the [110] and [010]
crystallographic axes, respectively. Furthermore, the surface exhibits four
minima on the $\langle$010$\rangle$ axes and four saddle-points along the
$\langle$110$\rangle$ directions.

Often, dispersion curves are plotted along a fixed propagation direction
$\theta$, corresponding to cuts across the dispersion surface
${\omega(k,\theta)}$. The curves in three different propagation directions are
depicted in Fig. 4A and they differ most in the directions 0∘ and 45∘.
Although pure Lamb/SH waves do not generally exist in the silicon plate
because they are coupled, we label the mode branches consistently to the
commonly used notation. In particular, following the notation in Ref. [17],
the positive-slope branch shown in Fig. 4A will be denoted as S1, while the
negative-slope section is called S2b.

The power flux of the waves is proportional to the group velocity [12, 32, 10,
11]. The latter is defined as the gradient of $\omega$ with respect to the
components of the wave vector, i.e.,

$\mathbf{c}{g}:=\nabla_{\mathbf{k}}\omega=\frac{\partial\omega}{\partial
k}\mathbf{e}_{r}+\frac{1}{k}\frac{\partial\omega}{\partial\theta}\mathbf{e}_{\theta}\,.$
(1)

It is worth remarking that these vectors are–per definition–orthogonal to the
iso-frequency lines drawn in Fig. 3A. This is illustrated in Fig. 4B, where
some group velocity vectors have been sketched for a chosen iso-frequency
contour.

The power flux is discussed by comparing propagation at the principal
direction $\theta=\leavevmode\nobreak\ $0∘ and the nonprincipal direction
$\theta=\leavevmode\nobreak\ $26.6∘. For both directions we compute the group
velocity vectors as a function of frequency. The resulting vector components
are depicted in Fig. 4C, while the magnitude is shown in Fig. 4D. Due to the
material’s reflection symmetry across principal axes, the transverse power
flux $c_{\mathup{g}\theta}$ is identically zero on principal directions. This
is not the case for other propagation directions, as Fig. 4C illustrates. Due
to this, the corresponding group velocity magnitude in Fig. 4D does no longer
vanish at the minimal frequency. Furthermore, depending on whether the axial
component $c_{\mathup{g}r}$ is positive or negative, we speak of _forward
waves_ or _backward waves_ [17, 33, 34, 35], respectively. With respect to a
source, a backward wave exhibits outgoing flux but incoming phase fronts. This
is the case for the S2b wave. Note that the phenomena investigated hereinafter
are strongly related to the existence of such waves.

The transverse group velocity component $c_{\mathup{g}\theta}$ is due to the
anisotropy $\partial\omega/\partial\theta$ and causes a transverse power flux.
We remark that it is the coupling of Lamb- and SH-polarizations that implies
this (usually) nonzero transverse power flux [25, Appendix B]. Due to this
transverse component, the overall power flux (or group velocity) is not
collinear with the wave vector. The angle between the two is denoted as
$\alpha$ and called _steering_ or _skew angle_ [12, 36, 37]. The frequency-
dependent skew angles of the S1/S2b modes are shown in Fig. 4E. It is
remarkable that they cover almost 180∘ in a rather narrow frequency range.

###  2.2 Zero- and transverse-group-velocity waves: ZGV and TGV

Resonances appear where the power flux vanishes and such points are denoted as
ZGV points. This requires the vanishing of the two components of the group
velocity vector given in (1), i.e., these resonances are associated to
stationary points of the dispersion surface of Fig. 3. It has been common
practice in the literature to plot the dispersion curves for a given
propagation direction as in Fig. 4A to identify ZGV modes, see Ref [24, 26,
27]. However, these curves only reveal the axial component
$c_{\mathup{g}r}=\partial\omega/\partial k$ of the group velocity vector.
Thus, identifying ZGV points as extrema of these curves is only valid for
isotropic media or for propagation along reflection symmetry planes of any
anistroptic material, because in these cases the group-velocity vector is
collinear to the wave vector, i.e.,
$c_{\mathup{g}\theta}\sim\partial\omega/\partial\theta\equiv 0$. Hence, while
the marked points on the curves of Fig. 4A appear to be ZGV points in the
conventional sense that $\partial\omega/\partial k=0$, their transverse power
flux might actually be nonzero due to anisotropy. For this reason we
distinguish between

  1. 1.

ZGV points: waves with overall zero power flow, i.e., $\partial\omega/\partial
k=0$ and $\partial\omega/\partial\theta=0$, and

  2. 2.

transverse-group-velocity (TGV) waves: waves with zero _axial_ power flux,
i.e., ${\partial\omega/\partial k=0}$.

A _ZGV point_ is an extremum or saddle point of the dispersion surface
depicted in Fig. 3. The cubic material exhibits eight isolated ZGV points on
the principal directions. The four minima in the $\langle$010$\rangle$
directions are denoted as ZGV1 and occur at $\omega/2\pi=\leavevmode\nobreak\
$7.7079 MHz and $k=\leavevmode\nobreak\ $3.421 rad/mm. The ZGV2 are saddle
points at $\langle$110$\rangle$ with $\omega/2\pi=\leavevmode\nobreak\ $7.7551
MHz and $k=\leavevmode\nobreak\ $3.142 rad/mm.

The existence of _TGV waves_ is evident from the closed iso-frequency contours
that enclose a minimum in Fig. 4B. These always exhibit points where the wave
vector is tangent to the contour and, consequently, orthogonal to the group
velocity. The wavenumber and frequency of TGV waves depend on the wave vector
orientation $\theta$ and form the continuous red curve drawn in Fig. 3A and
Fig. 3B. The closed TGV-curve separates the forward wave region (outside) from
the backward wave region (inside).

####  2.2.1 TGV wave radiated by a line source

The existence of TGV waves is quite remarkable as their power flux is
_orthogonal_ to the wave vector. While this fact is obvious from their
definition, the evolution of power flux when approaching the TGV-frequency can
be assessed in Fig. 4E: the skew angles of the S1 and S2b modes both converge
to $-$90∘ when reaching the TGV point. Note that the weaker the coupling
between Lamb and SH modes, the sharper we expect the transition towards $-$90∘
to take place. Ultimately, when these wave families decouple (i.e., on a
principal direction or in an isotropic plate), we obtain a discontinuity as
observed in Fig. 4E for $\theta=0^{\circ}$. The combination of orthogonal
propagation and wide skew-angle spectrum lead to extraordinary diffraction
effects.

Figure 5: Measured TGV wave field radiated by a synthetic line source inclined
26.6∘. (A) wavenumber spectrum at the TGV-frequency compared to theory prior
to filtering in the wavenumber domain; (B) instantaneous intensity
distributions after having applied a Gaussian filter in the wavenumber domain
as indicated in (A). The two wave packets propagate along the line source
while the wave vectors are orthogonal. Note that the phases propagate counter-
wise in the two pulses, leading to interference where the two overlap.
Frequency band: 7.6 MHz to 7.8 MHz. Full field reconstructed from one
quadrant. (cf. Supplementary Video 1.)

To observe these effects we synthesize a line-source response from the point-
source measurements, for details see Sec. 4.2. Waves usually radiate in normal
direction from the line source and diffract around its edges. However, this is
not the case for the TGV wave. This can be observed in Fig. 5, which presents
the response to a 6.4 mm line source preferentially exciting wave vectors at
26.6∘ and $-$153.4∘. The associated wavenumber spectrum at the TGV-frequency
is compared to the theoretical wavenumbers in Fig. 5A. To better observe the
TGV wave, we apply a Gaussian filter in the wavenumber plane. Its full width
at half maximum is delineated in Fig. 5A. A subsequent inverse Fourier
transform into the spatial domain yields the instantaneous intensity
distributions depicted in Fig. 5B. Two wave packets propagating along the line
source can clearly be discerned, i.e., their skew angle is indeed $\pm$90∘ as
expected. From another point of view, this wave “diffracts” around the line’s
edge while maintaining the orientation of the phase fronts aligned with the
line source.

A further validation is done by comparing the group velocity of the wave
packets to theory. After 40 µs the maxima of the envelops are 6.81 mm from the
source’s center, yielding a transverse velocity of $\pm$0.170 mm/µs. This is
in very good agreement with the computed transverse group velocity of $-$0.177
mm/µs for the TGV wave at $\theta=\leavevmode\nobreak\ $26.6∘, as shown in
Fig. 4C. Note that the wave vectors oriented at $\theta=-$153.4∘ lead to the
wave packet traveling in opposite direction.

Propagation of energy purely orthogonal to the wave vector is a strictly
monochromatic process, both in frequency and wavenumber. The neighboring
spectral components exhibit a nonzero axial power flux and these waves are
responsible for the packets’ lateral extend. Wave packets and their power flux
can only be observed when considering a spectrum of finite width. This is
ensured by the applied Gaussian filter in Fig. 5A, which has a half-width of
0.7 rad/mm. While all spectral components of the wave packets possess very
similar wave vectors, their power flux cover a wide range of directions.
Indeed, as the power flux is orthogonal to the iso-frequency lines, Fig. 5A
shows that the skew-angle spectrum spans almost 180∘.

Lastly, it is worth observing that the two wave packet’s phases propagate in
opposite directions, as indicated in Fig. 5B. The counter-propagating wave
packets lead to quickly varying interference in the regions where the two
superpose. It is remarkable that the phase fronts of one packet travel all in
the same direction. This means that the phases propagate _towards the source_
on one side, and _away from the source_ on the other side, as indicated by the
wave vector arrows in Fig. 5B. Accordingly, each wave packet clearly consists
of the S1 forward wave and the S2b backward wave, separated spatially by the
TGV component. This is consistent with the smaller wavelengths that we observe
in Fig. 5B for the S1 components compared to the S2b components. The spatial
allocation of modes is also in accordance with Fig. 4B as well as Fig. 4E.

####  2.2.2 Decay of the TGV-wave contributions

The decay observed for some spectral components in the introductory Fig. 1 can
now be explained. While the frequency extrema correspond to the ZGV
resonances, the remaining spectral components are due to TGV waves. The line
source excites the latter all along its length. The TGV waves then propagate
along the source, where they are measured at the center point by the
interferometer. Due to their propagative nature, the waves leave the detection
point after a certain time. As seen in Fig. 1B, the TGV waves disappear
completely after about 40 µs, while the expected eight ZGV resonances at
$\theta=n\times\text{45}^{\circ},n\in\mathbb{Z}$, perdure.

Figure 6: Formation of the S1/S2b-resonance pattern in an infinite silicon
plate. (A,C,E) show the wave vectors (not to scale); while (B,D,F) depict the
corresponding wave fields. The ZGV modes in (A,B) and (C,D) are characterized
by close but different wavenumbers $|\mathbf{k}|$ and frequencies $\omega$.
These two standing wave fields interfere to form a non-stationary, time-
dependent beating pattern as shown for example in (F) when the phase shift
between the two is zero. (cf. Supplementary Video 2.)

The time for the TGV waves to escape the measurement point depends both on the
length of the line source and the wave’s group velocity. Compared to other
wave components, the TGV wave exhibits a very low group velocity, as can well
be appreciated in Fig. 4D. Therefore, the TGV waves are detected for a
relatively long time, resulting in a well pronounced spectral peak in Fig. 1A.
This explains why they play an important role in the presented measurements.
Note that the group velocity tends to zero as the line orientation tends to
one of the material’s principal axes. Hence, the closer the line orientation
is to a ZGV point, the more pronounced the corresponding spectral peak will
be.

###  2.3 Resonance pattern of ZGV modes

The previous results demonstrate that time acts as a filter in the wave
vector-frequency domain. After sufficient time, only the ZGV resonances remain
in the finite spatial region of observation. As a consequence, a resonance
pattern develops in the silicon crystal plate. We first lay out the theory
explaining this resonance and subsequently we present direct measurements
thereof.

####  2.3.1 Theory of ZGV resonances

For the frequency region of interest, the resonance pattern is explained by
the interference of the eight ZGV resonances sketched in Fig. 6A and Fig. 6C.
The interference of the four wave vectors pertaining to a given ZGV mode leads
to the periodic square wave pattern depicted in Fig. 6B and Fig. 6D. Note that
for each wave vector $\mathbf{k}$, there exists a counter-propagating wave
$-\mathbf{k}$. As a consequence, the square patterns are standing wave fields
[17]. It is notable that this standing wave field is not due to the edges of
the plate but actually emerges in the _infinite plate_.

It is pertinent to discuss differences between the isotropic and anisotropic
cases. While the isotropic ZGV resonance consist of a continuous wave vector
spectrum, i.e., it exhibits wave vectors in all possible directions, the ZGV
resonances in an anisotropic medium consist of a finite set of wave vectors.
Accordingly, the ZGV field due to a point load on an isotropic plate is a
cylindrical standing wave [18]. In particular, this means that the nodal
curves seen on the plate’s surface are closed circles around the source that
never cross. In contrast to this, the field of the ZGV1 or ZGV2 resonances are
“checkerboard standing waves” as shown in Fig. 6B and 6D. Hence, anisotropy
leads to straight, open and crossing nodal lines.

The ZGV1 and ZGV2 waves interfere to form a resonance pattern as depicted in
Fig. 6F. The nodal curves form “bubbles” of different shapes that are no
longer square. As each of the wave components is flux free, no energy is
propagated. Nonetheless, this is no longer a standing wave field because the
ZGV1 and ZGV2 components are at different frequencies. Instead, the “bubbles”
move radially out or into the center of the pattern while changing their shape
(see Supplementary Video 2). Note that they shortly stop moving when the two
ZGV components are in phase or in paraphase. For the field synthesis shown in
Fig. 6F we arbitrarily assumed identical amplitudes and a zero phase shift
between ZGV1 and ZGV2. Due to the similarity of the displacement
eigenfunctions $\mathbf{u}{ZGV1}(z)$ and $\mathbf{u}{ZGV2}(z)$, both modes are
excited similarly by the point source. For the same reason, we expect both
ZGVs to be in phase at the moment of generation by the pulse.

The phase shift between ZGV1 and ZGV2 is modulated in time due to their
slightly different frequencies. In other words, temporal beating is expected
as a consequence of their superposition. Hence, the instantaneously observed
phase shift loops $2\pi$ in a period of

$\Delta T=\frac{2\pi}{\omega{ZGV2}-\omega{ZGV1}}=\text{{21.2\,\textmu}s}\,.$
(2)

####  2.3.2 Measurement of ZGV resonances

Figure 7: Measured wave fields at selected times. A resonance pattern forms
due to the superposition of the eight plane waves corresponding to ZGV
resonances in the single-crystal silicon wafer. The pattern forms naturally
over time as only the ZGV resonances remain in the spatial window close to the
point source. The inner region of 4 mm $\times$ 4 mm is marked for reference.
Frequency band: 7.6 MHz to 7.8 MHz. Full field reconstructed from one
quadrant. (cf. Supplementary Video 3 and 4.)

The described resonance pattern is directly observable by measuring the wave
field excited by a point-like source. As the S1/S2b resonance is the one that
is most strongly excited by our system, the theoretically expected moving
resonance pattern can be discerned in the raw data without postprocessing, see
Supplementary Video 4. To obtain a clean result, other resonances as well as
reflections from the border of the plate should be avoided. We achieve this
with the mentioned frequency bandpass filter in combination with the large
source that excites low wavenumbers, see Sec. 4.2 for details. The thus
observed frequency-wavenumber range contains only slow modes in the vicinity
of the ZGV resonances. The filtered displacement fields are depicted in Fig. 7
for selected time instants.

Propagating waves dominate the picture at first and their phase fronts (nodal
curves) encircle the source. But a nodal curve pattern with the mentioned
“bubbles” forms over time, see Fig. 7 and Supplementary Video 3. The pattern
develops slowly from the source (in the center) outwards as the propagative
waves leave the region of excitation. When the pattern becomes first visible,
the two ZGV resonances are not in phase and the “bubbles” move in and out on
the $\langle$110$\rangle$ and $\langle$100$\rangle$ axes as expected from the
theory.

After one beating period $\Delta T$, the two ZGVs should be in-phase for the
first time and, hence, we expect a pattern similar to Fig. 6F. At 21.25 µs the
pattern in the inner area marked by the square in Fig. 7 resembles the
corresponding region in Fig. 6F. Outside this region, the contribution of
propagating waves is still too high. Nonetheless, after two beating periods,
at 42.4µs, the non-ZGV waves propagated and we can now recognize a fully
developed ZGV resonance pattern in the region of observation. Note, for
example, the “eight-pointed star” centered at the center and at the edges of
the marked inner region. Discrepancies to Fig. 6F can mostly be attributed to
temporal sampling of the measurement.

Figure 8: Measured displacement fields at the two ZGV frequencies. (A,B)
harmonic field at arbitrary phase in the physical $X$-$Y$-space obtained by a
temporal Fourier transform. (C,D) spectral magnitude in the reciprocal
$k_{X}$-$k_{Y}$-space obtained by a spatio-temporal Fourier transform. The
numerically computed ZGV wave vectors are indicated by arrows. Full field
reconstructed from one quadrant.

A quantitative validation is performed by comparing measurement and theory in
the spectral domain. First, a time-Fourier transform yields the measured wave
fields at the ZGV frequencies. The result is depicted in Fig. 8A and Fig. 8B
and both confirm the wave fields expected in Figs. 6B and 6D. Second, we
additionally perform a spatial 2D Fourier transform into the wave vector-
frequency domain. The obtained spectral amplitudes are displayed in Fig. 8C
and Fig. 8D together with the ZGV wave vectors expected from theory. We
observe that the energy confines closely to the regions predicted by the
computed wave vectors. While the wave field at $\omega{ZGV1}$ consists almost
purely of the ZGV resonances, propagating modes exist at $\omega{ZGV2}$ and
are also excited. This is because $\omega{ZGV2}$ corresponds to the saddle-
points and the corresponding iso-frequency contours are two curves that
cross/touch at the ZGV2-points. In contrast to this, the iso-frequency
contours at $\omega{ZGV1}$ reduce to the isolated ZGV1-points. Note that
spatio-temporal gating can be used to remove the propagating waves of Fig. 8D,
as was done in Fig. 1. Overall, this confirms that the resonances observed in
an infinite plate of anisotropic elasticity are formed by the superposition of
a discrete set of ZGV modes.

##  3 Discussion

The existence of isolated critical points in the S1/S2b mode dispersion
surface entails two distinct and extraordinary effects: transverse wave
propagation (TGV waves) and beating ZGV resonances. We presented the
underlying theory as well as experimental evidence of both effects. TGV waves
are characterized by a power flux that is orthogonal to their wave vector.
They propagate along a line source that excites them. For this reason they
manifest as a peak in the spectrum measured at a chosen point of the line
source.

ZGV resonances are already used very successfully for ultrasonic
nondestructive testing and material characterization. With the developed
physical understanding, we expect similar procedures to become attractive for
novel microelectromechanical systems, in particular sensors. The present work
overcomes one of the major difficulties in designing such systems: properly
accounting for anisotropy, such as is encountered for silicon.

It is evident form the theory we presented that the ZGV resonances directly
reveal information about the principal axes of the material, independent of
the actual orientation of the sample. Measuring the ZGV1 and ZGV2 resonance
frequencies as well as the frequency of a regular thickness resonance is
sufficient to fully characterize a material of cubic symmetry. These
frequencies are readily obtained from the recorded normal surface displacement
after a single pulsed excitation because they all manifest as sharp peaks in
the spectrum [24]. This avoids the technically complicated scan of the wave
field, as needs to be done for widespread techniques based on guided waves.
Overall, this could enable nondestructive, contactless and real-time
characterization of complex materials.

The present contribution restricted to ZGV resonances of the S1/S2b modes.
However, it is well known that other ZGV modes exist at higher-frequencies.
They are substantially weaker and were frequency-filtered in this work. Note
that anisotropic plates can also support multiple ZGV resonances pertaining to
the same modes [25]. In this case, not only local minima, but also local
maxima can appear in the dispersion surface. The analysis and mechanics
presented here apply analogously.

The discussed phenomena are very particular to two-dimensional (and partly
also three-dimensional) wave propagation in anisotropic media. Analogous
phenomena should be expected when the material pertains to a different
anisotropy class than the one studied here. Moreover, similar effects are
expected in phononic crystals or meta-materials, as dispersion surfaces with
similar features haven been reported for these materials [38, 39, 40, 41]. But
note that the critical points (ZGV points) that we studied here exist even
without periodicity, i.e., in a homogeneous anisotropic plate. Lastly,
hyperbolic polaritons share many of the wave propagation features described
here, including power flux orthogonal to the wave vector [42].

##  4 Methods

The procedure that we use to compute guided elastic waves is sketched in
Subsec. 4.1. The experiment is detailed in Subsec. 4.2.

###  4.1 Computing guided waves and power flux

The waves propagate in a plate that is assumed to be infinite in the
$X$-$Y$-plane, see Fig. 2A. Its material is of anisotropic stiffness
$\mathbf{c}$ and mass density $\rho$. The plane waves of wavenumber $k$ shall
be harmonic with angular frequency $\omega$. Hence, their wave field is of the
form

$\mathbf{\tilde{u}}(r,\theta,z,t)=\mathbf{u}(z)\operatorname{e}^{\mathup{i}\mkern
1.0mukr}\operatorname{e}^{-\mathup{i}\mkern 1.0mu\omega t}\,.$ (3)

We proceed by arbitrarily fixing the propagation vector
$\mathbf{k}=k\mathbf{e}_{r}$ and computing the corresponding frequency
$\omega$ and modal displacements $\mathbf{u}(z)$. These are obtained as the
solutions of a differential eigenvalue problem [43, 29, 25]. $\omega$ and
$\mathbf{u}(z)$ represent the eigenvalues and the eigenfunctions,
respectively. For a concise derivation of the concrete problem formulation see
Ref. [25].

We compute solutions with a _semi-analytical procedure_ that consists in two
steps: (i) discretize the differential eigenvalue problem and (ii) solve the
resulting algebraic eigenvalue problem numerically. Many variants of this
procedure have been discussed in the literature [43, 44, 45]. A concise
implementation based on spectral collocation is GEW dispersion script [46].
For the current work, we perform the discretization using the _spectral
element method_ , i.e., one finite element of high polynomial order [47]. This
method leads to a regular Hermitian eigenvalue problem [48], which allows us
to reliably compute the ZGV and TGV points in general elastic waveguides [25,
31]. We make our implementations available under the name of GEWtool [30] and
it includes all required post-processing methods to reproduce the results of
this contribution. In particular, the script dispersionSurface_silicon_ZGV.m
in the examples directory produces Fig. 3A.

The coupling and decoupling of wave families is of importance [29, 25, 49,
48]. Symmetric and anti-symmetric waves decouple due to the plate’s symmetry
across its mid-plane. We only compute the symmetric waves by modeling the top
half of the plate and fixing the $u_{z}$ displacement component at the mid-
plane. Moreover, the shear-horizontal (SH) polarization ($u_{\theta}$) only
decouples form the Lamb polarization ($u_{r},u_{z}$) for
$\theta=n\times\text{45}^{\circ},n\in\mathbb{Z}$. For this reason, we always
compute the _fully-coupled waves_ , meaning that all three displacement
components are accounted for in the displacement vector $\mathbf{u}(z)$.

After computing guided wave solutions, their power flux and group velocity can
be computed in a post-processing step. To this end, we exploit the fact that
in nondissipative waveguides (more precisely: real-valued $\mathbf{k}$ and
$\omega$) the group velocity is equal to the _energy velocity_ [32, 10, 50]
and is given by

$\mathbf{c}{g}=\mathbf{c}{e}=\frac{\int\mathbf{p}\,\mathup{d}z}{\int\mathcal{H}\,\mathup{d}z}\,.$
(4)

It is defined through the power flux density vector $\mathbf{p}$ and the
average total stored energy density $\mathcal{H}$. Using the particle velocity
$\mathbf{v}=-\mathup{i}\mkern 1.0mu\omega\mathbf{u}$ and stress
$\mathbf{T}={\mathbf{c}:\nabla\mathbf{u}}={\mathbf{c}:(\mathup{i}\mkern
1.0muk\mathbf{e}_{r}+\mathbf{e}_{z}\partial_{z})\mathbf{u}}$, the power flux
density can be computed as

$\displaystyle\mathbf{p}$
$\displaystyle=-\frac{1}{2}\operatorname{Re}\\{\mathbf{v}^{*}\cdot\mathbf{T}\\}$
$\displaystyle=-\frac{1}{2}\operatorname{Re}\\{\mathup{i}\mkern
1.0mu\omega\mathbf{u}^{*}\cdot\mathbf{c}:(\mathup{i}\mkern
1.0muk\mathbf{e}_{r}+\mathbf{e}_{z}\partial_{z})\mathbf{u}\\}\,,$ (5)

where “∗” denotes complex conjugation, “$\operatorname{Re}$” the real-part
operator and $\partial_{z}=\partial/\partial z$. Furthermore, due to
equipartition of energy [12], we can compute the average total stored energy
density by

$\mathcal{H}=\frac{1}{2}\rho\omega^{2}\mathbf{u}^{*}\cdot\mathbf{u}\,.$ (6)

Our previously outlined procedure to solve for guided waves yields all
quantities required to compute equation (5) and equation (6). The
differentiation in equation (5) and the integration in equation (4) can be
performed numerically. In this way, we can compute the energy velocity of each
point of the dispersion curves independently. Most importantly, proceeding in
this way explicitly provides both the axial and the transverse components of
the energy velocity vector. The $z$-component is identically zero due to the
power-flux-free surfaces of the plate.

Lastly, as the wave vector is $\mathbf{k}=k\mathbf{e}_{r}$, the _skew angle_
$\alpha$ can be computed from the energy velocity vector as

$\alpha=\begin{cases}\arctan(c_{\mathup{e}\theta}/c_{\mathup{e}r})&\text{for
}c_{\mathup{e}r}\geq 0\,,\\\
\arctan(c_{\mathup{e}\theta}/c_{\mathup{e}r})-\pi&\text{otherwise}\,.\end{cases}$
(7)

###  4.2 Measurement setup

Measurements are achieved all-optically and are presented schematically in
Fig. 2B. Guided waves are generated with a Q-switched Nd:YAG (yttrium aluminum
garnet) laser (Quantel Laser, France, Centurion, 1064 nm optical wavelength,
100 Hz repetition rate) that delivers 10 ns pulses of 9 mJ. The laser output
beam is expanded and then focused onto the silicon plate with a lens (100 mm
focal length). The focal spot is kept rather large (beam width $\approx$1 mm),
thus exciting wavenumbers up to $\approx$7 rad/mm. This favors the generation
of the first ZGV resonance while limiting the generation of the fundamental
guided modes in the frequency range of interest [51]. The wafer (525 µm
nominal thickness, 524.6 µm measured thickness, 125 mm diameter) has a 100 nm
aluminum coating on the excitation side, which reduces the optical penetration
depth. However, similar results were obtained when exciting on the side
without the coating. The coating is thin enough not to affect the elastic
waves in the silicon plate.

The normal surface displacement is detected on the opposite side with a
heterodyne interferometer of 532 nm optical wavelength with a power of 100 mW
and a focal spot of $\approx$50 µm. A hardware high-pass filter with 1.25 MHz
cut-off frequency is used to avoid saturation of the interferometer due to the
large low-frequency displacements of the A0 mode. Signals are recorded with
100 MHz sampling rate by an oscilloscope connected to a computer. Each signal
is averaged 128 times. The 25.05 mm $\times$ 25.05 mm scan is performed by
moving the excitation unit with a two-axes translation stage along
$\mathbf{e}_{X}$ and $\mathbf{e}_{Y}$ with a 150 µm pitch. Note that scanning
the field is time consuming. For this reason we exploited the cubic symmetry
of the material and scanned only one quadrant around the source, as indicated
in the inset of Fig. 2B. The full fields have been reconstructed by an
appropriate symmetrization.

While the measurements directly provide the point-source response, the line-
source response shown in Fig. 5 is synthesized a posteriori. To this end, we
superpose 20 shifted point-source responses. The synthetic line source
exciting wave vectors at $\theta=\leavevmode\nobreak\ $26.6∘ is obtained by
shifting two pitches vertically for every horizontal one.

Fast waves outside the ZGV region are also excited and their reflections from
the border of the plate disturb the long-time observation of the slow TGV
waves in Fig. 5 as well as the resonance patterns in Fig. 7. To avoid this, we
use a band-pass filtered between 7.6 MHz and 7.8 MHz for these figures. Note
that due to the large source size, large wavenumbers are not observed. The
remaining frequency-wavenumber range contains only the slow modes close to the
ZGV points.

## References

  * [1] Matthew A. Hopcroft, William D. Nix and Thomas W. Kenny  “What is the Young’s Modulus of Silicon?”  In _J. Microelectromech. Syst._ 19.2, 2010, pp. 229–238  DOI: 10.1109/JMEMS.2009.2039697
  * [2] A. Belyaev et al.  “Crack detection and analyses using resonance ultrasonic vibrations in full-size crystalline silicon wafers”  In _Appl. Phys. Lett._ 88.11, 2006, pp. 111907  DOI: 10.1063/1.2186393
  * [3] Sunil Kishore Chakrapani, M. Padiyar and Krishnan Balasubramaniam  “Crack Detection in Full Size Cz-Silicon Wafers Using Lamb Wave Air Coupled Ultrasonic Testing (LAC-UT)”  In _J. Nondestr. Eval._ 31.1, 2012, pp. 46–55  DOI: 10.1007/s10921-011-0119-3
  * [4] Min-Kyoo Song and Kyung-Young Jhang  “Crack Detection in Single-Crystalline Silicon Wafer Using Laser Generated Lamb Wave”  In _Adv. Mater. Sci. Eng._ 2013, 2013, pp. e950791  DOI: 10.1155/2013/950791
  * [5] Bernard Masserey, Mathieu Simon, Jean-Luc Robyr and Paul Fromme  “Defect detection in monocrystalline silicon wafers using high frequency guided waves”  In _AIP Conference Proceedings_ 2102.1, 2019, pp. 020013  DOI: 10.1063/1.5099717
  * [6] L. Hackett et al.  “Aluminum scandium nitride films for piezoelectric transduction into silicon at gigahertz frequencies”  In _Appl. Phys. Lett._ 123.7, 2023, pp. 073502  DOI: 10.1063/5.0151434
  * [7] Ventsislav Yantchev, Lilia Arapan, Ilia Katardjiev and Victor Plessky  “Thin-film zero-group-velocity Lamb wave resonator”  In _Appl. Phys. Lett._ 99.3, 2011, pp. 033505  DOI: 10.1063/1.3614559
  * [8] C. Caliendo and M. Hamidullah  “Zero-group-velocity acoustic waveguides for high-frequency resonators”  In _J. Phys. D: Appl. Phys._ 50.47, 2017, pp. 474002  DOI: 10.1088/1361-6463/aa900f
  * [9] Stefan J. Rupitsch  “Piezoelectric Sensors and Actuators - Fundamentals and Applications”, Topics in Mining, Metallurgy and Materials Engineering  Berlin: Springer, 2018
  * [10] B.. Auld  “Acoustic Fields and Waves in Solids”  Malabar: Krieger Publishing Company, 1990
  * [11] Daniel Royer and Tony Valier-Brasier  “Elastic Waves in Solids 1: Propagation”  London: ISTEJohn Wiley & Sons, 2022
  * [12] Karl-Jörg Langenberg, René Marklein and Klaus Mayer  “Theoretische Grundlagen der zerstörungsfreien Materialprüfung mit Ultraschall”  Munich: De Gruyter Oldenbourg, 2009
  * [13] C. Prada, O. Balogun and T.. Murray  “Laser-based ultrasonic generation and detection of zero-group velocity Lamb waves in thin plates”  In _Appl. Phys. Lett._ 87.19, 2005, pp. 194109  DOI: 10.1063/1.2128063
  * [14] Dominique Clorennec, Claire Prada, Daniel Royer and Todd W. Murray  “Laser impulse generation and interferometer detection of zero group velocity Lamb mode resonance”  In _Appl. Phys. Lett._ 89.2, 2006, pp. 024101  DOI: 10.1063/1.2220010
  * [15] Dominique Clorennec, Claire Prada and Daniel Royer  “Local and noncontact measurements of bulk acoustic wave velocities in thin isotropic plates and shells using zero group velocity Lamb modes”  In _J. Appl. Phys._ 101.3, 2007, pp. 034908  DOI: 10.1063/1.2434824
  * [16] O. Balogun, T.. Murray and C. Prada  “Simulation and measurement of the optical excitation of the S1 zero group velocity Lamb wave resonance in plates”  In _J. Appl. Phys._ 102.6, 2007, pp. 064914  DOI: 10.1063/1.2784031
  * [17] Claire Prada, Dominique Clorennec and Daniel Royer  “Local vibration of an elastic plate and zero-group velocity Lamb modes”  In _J. Acoust. Soc. Am._ 124.1, 2008, pp. 203–212  DOI: 10.1121/1.2918543
  * [18] Qingnan Xie et al.  “Imaging gigahertz zero-group-velocity Lamb waves”  In _Nat. Commun._ 10.1, 2019, pp. 2228  DOI: 10.1038/s41467-019-10085-4
  * [19] M. Cès, D. Royer and C. Prada  “Characterization of mechanical properties of a hollow cylinder with zero group velocity Lamb modes”  In _J. Acoust. Soc. Am._ 132.1, 2012, pp. 180–185  DOI: 10.1121/1.4726033
  * [20] Clemens Grünsteidl, Todd W. Murray, Thomas Berer and István A. Veres  “Inverse characterization of plates using zero group velocity Lamb modes”  In _Ultrasonics_ 65, 2016, pp. 1–4  DOI: 10.1016/j.ultras.2015.10.015
  * [21] Clemens Grünsteidl, Thomas Berer, Mike Hettich and István Veres  “Determination of thickness and bulk sound velocities of isotropic plates using zero-group-velocity Lamb waves”  In _Appl. Phys. Lett._ 112.25, 2018, pp. 251905  DOI: 10.1063/1.5034313
  * [22] Sylvain Mezil, Jérôme Laurent, Daniel Royer and Claire Prada  “Non contact probing of interfacial stiffnesses between two plates by zero-group velocity Lamb modes”  In _Appl. Phys. Lett._ 105.2, 2014, pp. 021605  DOI: 10.1063/1.4890110
  * [23] Martin Ryzy et al.  “Determining longitudinal and transverse elastic wave attenuation from zero-group-velocity Lamb waves in a pair of plates”  In _J. Acoust. Soc. Am._ 153.4, 2023, pp. 2090  DOI: 10.1121/10.0017652
  * [24] Claire Prada, Dominique Clorennec, Todd W. Murray and Daniel Royer  “Influence of the anisotropy on zero-group velocity Lamb modes”  In _J. Acoust. Soc. Am._ 126.2, 2009, pp. 620–625  DOI: 10.1121/1.3167277
  * [25] Daniel A. Kiefer, Bor Plestenjak, Hauke Gravenkamp and Claire Prada  “Computing zero-group-velocity points in anisotropic elastic waveguides: Globally and locally convergent methods”  In _J. Acoust. Soc. Am._ 153.2, 2023, pp. 1386–1398  DOI: 10.1121/10.0017252
  * [26] E.. Glushkov and N.. Glushkova  “Multiple zero-group velocity resonances in elastic layered structures”  In _J. Sound Vib._ 500, 2021, pp. 116023  DOI: 10.1016/j.jsv.2021.116023
  * [27] Takasar Hussain and Faiz Ahmad  “Lamb modes with multiple zero-group velocity points in an orthotropic plate”  In _J. Acoust. Soc. Am._ 132.2, 2012, pp. 641–645  DOI: 10.1121/1.4730891
  * [28] Marc Thelen et al.  “Laser-excited elastic guided waves reveal the complex mechanics of nanoporous silicon”  In _Nat. Commun._ 12.1, 2021, pp. 3597  DOI: 10.1038/s41467-021-23398-0
  * [29] Daniel A. Kiefer  “Elastodynamic quasi-guided waves for transit-time ultrasonic flow metering”, FAU Forschungen, Reihe B, Medizin, Naturwissenschaft, Technik 42  Erlangen: FAU University Press, 2022  DOI: 10.25593/978-3-96147-550-6
  * [30] Daniel A. Kiefer  “GEWtool” DOI: 10.5281/zenodo.10114244, 2023  DOI: 10.5281/zenodo.10114244
  * [31] Bor Plestenjak and Daniel A. Kiefer  “GEW ZGV computation” DOI: 10.5281/zenodo.7537441, 2023  DOI: 10.5281/zenodo.7537441
  * [32] Karl J. Langenberg, Rene Marklein and Klaus Mayer  “Energy vs. group velocity for elastic waves in homogeneous anisotropic solid media”  In _2010 URSI International Symposium on Electromagnetic Theory_ Berlin, Germany: IEEE, 2010, pp. 733–736  DOI: 10.1109/URSI-EMTS.2010.5637253
  * [33] A.. Shuvalov and O. Poncelet  “On the backward Lamb waves near thickness resonances in anisotropic plates”  In _Int. J. Solids Struct._ 45.11, 2008, pp. 3430–3448  DOI: 10.1016/j.ijsolstr.2008.02.004
  * [34] Suraj Bramhavar et al.  “Negative refraction and focusing of elastic Lamb waves at an interface”  In _Phys. Rev. B_ 83.1, 2011, pp. 014106  DOI: 10.1103/PhysRevB.83.014106
  * [35] François Legrand et al.  “Cloaking, trapping and superlensing of lamb waves with negative refraction”  In _Sci. Rep._ 11.1, 2021, pp. 23901  DOI: 10.1038/s41598-021-03146-6
  * [36] Bastien Chapuis, Nicolas Terrien and Daniel Royer  “Excitation and focusing of Lamb waves in a multilayered anisotropic plate”  In _J. Acoust. Soc. Am._ 127.1, 2010, pp. 198–203  DOI: 10.1121/1.3263607
  * [37] Alexander Karmazin, Evgenia Kirillova, Wolfgang Seemann and Pavel Syromyatnikov  “A study of time harmonic guided Lamb waves and their caustics in composite plates”  In _Ultrasonics_ 53.1, 2013, pp. 283–293  DOI: 10.1016/j.ultras.2012.06.012
  * [38] Aleksi Bossart and Romain Fleury  “Extreme Spatial Dispersion in Nonlocally Resonant Elastic Metamaterials”  In _Phys. Rev. Lett._ 130.20, 2023, pp. 207201  DOI: 10.1103/PhysRevLett.130.207201
  * [39] A A Maznev, O B Wright and O Matsuda  “Mapping the band structure of a surface phononic crystal”  In _New J. Phys._ 13.1, 2011, pp. 013037  DOI: 10.1088/1367-2630/13/1/013037
  * [40] Nadège Kaina, Fabrice Lemoult, Mathias Fink and Geoffroy Lerosey  “Negative refractive index and acoustic superlens from multiple scattering in single negative metamaterials”  In _Nature_ 525.7567, 2015, pp. 77–81  DOI: 10.1038/nature14678
  * [41] Simon Yves, Fabrice Lemoult, Mathias Fink and Geoffroy Lerosey  “Crystalline Soda Can Metamaterial exhibiting Graphene-like Dispersion at subwavelength scale”  In _Sci. Rep._ 7.1, 2017, pp. 15359  DOI: 10.1038/s41598-017-15335-3
  * [42] Xin Zhang et al.  “Ultrafast anisotropic dynamics of hyperbolic nanolight pulse propagation”  In _Sci. Adv._ 9.34, 2023, pp. eadi4407  DOI: 10.1126/sciadv.adi4407
  * [43] A… Adamou and R.. Craster  “Spectral methods for modelling guided waves in elastic media”  In _J. Acoust. Soc. Am._ 116.3, 2004, pp. 1524–1535  DOI: 10.1121/1.1777871
  * [44] Hauke Gravenkamp, Chongmin Song and Jens Prager  “A numerical approach for the computation of dispersion relations for plate structures using the Scaled Boundary Finite Element Method”  In _J. Sound Vib._ 331.11, 2012, pp. 2543–2557  DOI: 10.1016/j.jsv.2012.01.029
  * [45] S. Finnveden  “Evaluation of modal density and group velocity by a finite element method”  In _J. Sound Vib._ 273.1, 2004, pp. 51–75  DOI: 10.1016/j.jsv.2003.04.004
  * [46] Daniel A. Kiefer  “GEW dispersion script” DOI: 10.5281/zenodo.7010603, 2022  DOI: 10.5281/zenodo.7010603
  * [47] Hauke Gravenkamp, Albert A. Saputra and Sascha Duczek  “High-Order Shape Functions in the Scaled Boundary Finite Element Method Revisited”  In _Arch. Comput. Methods Eng._ 28.2, 2021, pp. 473–494  DOI: 10.1007/s11831-019-09385-1
  * [48] Hauke Gravenkamp, Bor Plestenjak and Daniel A. Kiefer  “Notes on osculations and mode tracing in semi-analytical waveguide modeling”  In _Ultrasonics_ , 2023, pp. 107112  DOI: 10.1016/j.ultras.2023.107112
  * [49] F. Hernando Quintanilla, M… Lowe and R.. Craster  “The symmetry and coupling properties of solutions in general anisotropic multilayer waveguides”  In _J. Acoust. Soc. Am._ 141.1, 2017, pp. 406–418  DOI: 10.1121/1.4973543
  * [50] M.. Biot  “General Theorems on the Equivalence of Group Velocity and Energy Transport”  In _Phys. Rev._ 105.4, 1957, pp. 1129–1137  DOI: 10.1103/PhysRev.105.1129
  * [51] François Bruno et al.  “Laser beam shaping for enhanced Zero-Group Velocity Lamb modes generation”  In _J. Acoust. Soc. Am._ 140.4, 2016, pp. 2829–2838  DOI: 10.1121/1.4965291

## Acknowledgements

Funding: This work has received support under the program «Investissements
d’Avenir» launched by the French Government under Reference No.
ANR-10-LABX-24.

Author contributions: DAK and CP conceived and designed the study. SM
conceived and performed the experiment and the data curation. DAK performed
the modeling and implemented the numerical model. All authors analyzed and
interpreted the experimental results. DAK and SM wrote the manuscript. All
authors proofread and revised the manuscript.

Competing interests: The authors declare no competing interests.

Data and material availability: All data needed to evaluate the conclusions in
the paper are present in the paper and/or the Supplementary Materials. Our
code is provided at https://dx.doi.org/10.5281/zenodo.10114244.

## Supplementary Materials

We provide four supplementary videos showing the spatio-temporal dynamics of
the studied waves.


On a problem involving the squares of odd harmonic numbers

John M. Campbell111Corresponding author<EMAIL_ADDRESS>
Department of Mathematics and Statistics, York University, 4700 Keele St,
Toronto, Ontario, Canada M3J 1P3

Paul Levrie

Faculty of Applied Engineering, UAntwerpen, Groenenborgerlaan 171, 2020
Antwerp, Belgium, and Department of Computer Science, KU Leuven,
Celestijnenlaan 200A, 3001 Heverlee (Leuven), Belgium, P. O. Box 2402

Ce Xu

School of Mathematics and Statistics, Anhui Normal University, Wuhu 241002,
People’s Republic of China

Jianqiang Zhao

Department of Mathematics, The Bishop’s School, La Jolla, CA 92037, USA

###### Abstract

We solve a problem recently considered by Wang and Chu concerning series
involving the squares of finite sums of the form
$1+\frac{1}{3}+\cdots+\frac{1}{2n-1}$; many closely related and similarly
nontrivial series and integral evaluations are also introduced.  Our main
proof involves the use of an identity derived from Wilf–Zeilberger theory and
techniques concerning the theory of colored multiple zeta values.

_MSC:_ primary 11Y60 secondary 33B15

_Keywords:_ Harmonic sum, odd harmonic number, closed form, infinite series,
colored multiple zeta value

##  1 Introduction

As Adamchik has expressed [1] (cf. [7]), there is a great history concerning
the study of series of the following form:

$S(r)=\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{k+r}.$
(1)

As in [1, 7], we recall that Ramanujan was interested in the $S$-function
shown above [3, p. 39], and we record Ramanujan’s identity whereby

$S(r)=\frac{16^{r}}{\pi
r^{2}\binom{2r}{r}^{2}}\sum_{k=0}^{r-1}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}$
(2)

for $r\in\mathbb{N}_{0}$. Ramanujan’s discoveries concerning the $S$-function
in (1), together with the large amount of research, over the years, concerning
this function [1], have inspired the development of _Ramanujan-like_ series
given by introducing summand factors, such as harmonic or harmonic-type
numbers, in the series shown in (1) [4, 15]. In this article, we solve a
problem recently considered by Wang and Chu [16] concerning a family of series
resembling (1) and involving harmonic-type numbers, and we introduce many
related evaluations.

We write $O_{n}=1+\frac{1}{3}+\cdots+\frac{1}{2n-1}$ to denote the
$n^{\text{th}}$ odd harmonic number, by analogy with the usual notation
$H_{n}=1+\frac{1}{2}+\cdots+\frac{1}{n}$ for the $n^{\text{th}}$ entry in the
classical sequence of harmonic numbers. In the recent article [16], Wang and
Chu left the problem of evaluating series of the form

$\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{(1+2n-2\lambda)^{2}}$
(3)

as an open problem, letting $\lambda\in\mathbb{N}$. Even for the base case
such that $\lambda=1$, evaluating the series in (3) seems to be quite
difficult, and it appears that the coefficient-extraction methods of Wang and
Chu [15, 16] do not apply to (3). We offer explicit solutions to Wang and
Chu’s problem for the  cases whereby $\lambda\in\\{1,2,3,4\\}$, using recent
results from references as in [6, 17], and,  furthermore, our techniques, as
in Section 3, allow us to express (3) as the sum of

$\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{(1+2n-2(\lambda-1))^{2}}$

plus expressions that admit explicit finite sum evaluations for
$\lambda\in\mathbb{N}$, thus providing a complete solution to Wang and Chu’s
open problem on (3). We also take care of two other open problems in [16],
namely, the evaluation of:

$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}}{(1+2n-2\lambda)^{2}},$
(4)
$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{1+2n-2\lambda}.$
(5)

###  1.1 Organization of the article

In Section 2.1, we introduce a full proof for a symbolic evaluation for Wang
and Chu’s series in (3), for the $\lambda=1$ case. Our strategy in Section 2.1
has the advantage of producing many new and nontrivial series/integral
evaluations closely related to (3), and is also of interest in terms of the
many different mathematical techniques involved, including the Wilf–Zeilberger
method [13], and with reference to recent advances in the study of Apéry-type
series and colored multiple zeta values [17]. Furthermore, with a particular
regard to the formula in (18), methods and formulas involved in our proof in
Section 2.1 are required in generalizing our solution for the $\lambda=1$ case
so as to be applicable for $\lambda>1$. In Section 2.2, we consider a method,
based on the recent work of Xu and Zhao [17], for generalizing our evaluation
for (3), according to higher powers of odd harmonic numbers, again in the
$\lambda=1$ case. In Section 3, we succeed in determining a recursion for (3)
of the form described above, and we provide explicit evaluations for (3) for
the $\lambda\in\\{2,3,4\\}$ cases. Many closely related integral/series
identities are introduced in Section 4.

##  2 The $\lambda=1$ case

We adopt notation from [8], writing

$\mathcal{G}:=\Im\left(\text{Li}_{3}\left(\frac{i+1}{2}\right)\right)$

to denote the Catalan-like constant explored in [8], letting
$G:=\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{(2n-1)^{2}}$ denote the “original”
Catalan constant, and writing
$\text{Li}_{s}(z)=z+\frac{z^{2}}{2^{s}}+\frac{z^{3}}{3^{s}}+\cdots$ to denote
the polylogarithm function.

###  2.1 A full proof

We again emphasize the importance of our evaluation in (18) in the following
proof, which we use in a direct way to obtain an explicit solution for the
$\lambda=2$ case that may be generalized for $\lambda>2$.

######  Theorem 1.

The infinite series

$\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}^{2}}{(2k-1)^{2}}$

admits the symbolic form

$-\frac{8G}{\pi}-\frac{8G\ln(2)}{\pi}-\frac{48\mathcal{G}}{\pi}+\frac{9\pi^{2}}{8}+\frac{\pi}{6}+\frac{4\ln^{2}(2)}{\pi}+\frac{3\ln^{2}(2)}{2}.$

###### Proof.

Using the Wilf–Zeilberger method [13], Campbell [6] had recently proved the
following equality, noting the appearance of a copy of the base case of (3)
for $\lambda=1$:

$\displaystyle
4\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}}{2k+1}-\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}^{2}}{2k-1}-2\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}^{2}}{(2k-1)^{2}}$
$\displaystyle\hskip
28.45274pt=\frac{12G}{\pi}+\frac{32\mathcal{G}}{\pi}-\frac{6\ln^{2}(2)}{\pi}-\frac{3\pi^{2}}{4}-\frac{\pi}{4}-\ln^{2}(2).$
(6)

Using a recursive proof approach, Campbell [6] had also evaluated the
“central” series in (6):

$\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}^{2}}{2k-1}=\frac{4G}{\pi}-\frac{\pi}{12}-\frac{2\ln^{2}(2)}{\pi},$
(7)

giving us that the following equality holds:

$\displaystyle
4\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}}{2k+1}-2\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}^{2}}{(2k-1)^{2}}$
(8)
$\displaystyle=\frac{16G}{\pi}+\frac{32\mathcal{G}}{\pi}-\frac{3\pi^{2}}{4}-\frac{\pi}{3}-\frac{8\ln^{2}(2)}{\pi}-\ln^{2}(2).$

So, the problem of evaluating the base case for (3) is equivalent to the
problem of evaluating

$\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}}{2k+1}.$

In this direction, we claim that the following equality holds:

$\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{2O_{k}-H_{k}}{2k+1}=-\frac{16\mathcal{G}}{\pi}+\frac{\pi^{2}}{8}+\frac{8G\ln(2)}{\pi}+\frac{\ln^{2}(2)}{2}.$
(9)

To prove this, we need the following integrals:

$\displaystyle\int_{0}^{\pi/4}(\ln(\sin(t))-\ln(\cos(t)))dt=-G,$ (10)
$\displaystyle\int_{0}^{\pi/4}(\ln^{2}(\sin(t))-\ln^{2}(\cos(t)))dt=\frac{5\pi^{3}}{64}+\frac{1}{16}\pi\ln^{2}(2)-2\mathcal{G}+G\ln(2).$
(11)

The first is a well known expression for Catalan’s constant, whereas the
second one can be found in [8]. We are to also use the following result [12]:

$\int_{0}^{\arcsin(x)}\ln^{n}(\sin
u)\,{du}=\sum_{i=0}^{n}\left[(-1)^{i}i!\binom{n}{i}\ln^{n-i}(x)\sum_{k=0}^{\infty}\frac{\binom{2k}{k}\,x^{2k+1}}{(2k+1)^{i+1}\,2^{2k}}\right].$
(12)

For $n=1$, (12) reduces to:

$\int_{0}^{\arcsin(x)}\ln(\sin
u)\,{du}=\ln(x)\sum_{k=0}^{\infty}\frac{\binom{2k}{k}\,x^{2k+1}}{(2k+1)\,2^{2k}}-\sum_{k=0}^{\infty}\frac{\binom{2k}{k}\,x^{2k+1}}{(2k+1)^{2}\,2^{2k}}.$

We divide by $x$, replace $x$ by $\sin v$, and integrate the result between
$0$ and $\tfrac{\pi}{2}$:

$\displaystyle\int_{0}^{\pi/2}\left(\frac{1}{\sin v}\int_{0}^{t}\ln(\sin
u)\,{du}\right){dv}=$
$\displaystyle\sum_{k=0}^{\infty}\frac{\binom{2k}{k}}{(2k+1)\,2^{2k}}\int_{0}^{\pi/2}\ln(\sin
v)\sin^{2k}v\,dv-\frac{\pi}{2}\sum_{k=0}^{\infty}\frac{\binom{2k}{k}^{2}}{(2k+1)^{2}\,2^{4k}}.$
(13)

In the double integral on the left-hand side, we change the order of
integration:

$\displaystyle\int_{0}^{\pi/2}\left(\frac{1}{\sin v}\int_{0}^{t}\ln(\sin
u)\,{du}\right){dv}$ $\displaystyle=\int_{0}^{\pi/2}\left(\ln(\sin
u)\int_{u}^{\pi/2}\frac{1}{\sin v}\,{dv}\right){du}$
$\displaystyle=-\int_{0}^{\pi/2}\ln(\sin u)\ln(\tan\tfrac{u}{2}){du}.$

Using the substitution $u=2t$ the integral can be rewritten as:

$\displaystyle-2\int_{0}^{\pi/4}\left[\ln
2\cdot(\ln(\sin(t))-\ln(\cos(t)))+\ln^{2}(\sin(t))-\ln^{2}(\cos(t))\right]{dt}$
(14) $\displaystyle\hskip
85.35826pt=-\frac{5}{32}\pi^{3}-\frac{1}{8}\pi\ln^{2}2+4\mathcal{G}$

using (10) and (11). To calculate the integral in the right-hand side of (13),
we replace $\sin v$ by $x$ and use integration by parts. This leads to the
following recurrence:

$I_{k}:=\int_{0}^{\pi/2}\ln(\sin v)\sin^{2k}v\,dv\ \ \mbox{satisfies:}\hskip
28.45274pt\
2kI_{k}=\frac{\pi}{2}\cdot\frac{1}{2k}\cdot\frac{\binom{2k-2}{k-1}}{2^{2k-2}}+(2k-1)I_{k-1}$

with initial value $I_{0}=-\frac{1}{2}\pi\ln 2$. If we put
$I_{k}=\frac{\pi}{2}\frac{1}{2^{2k}}\binom{2k}{k}Z_{k}$, $Z_{k}$ is a solution
of the recurrence $Z_{k}=\frac{1}{2k-1}-\frac{1}{2k}+Z_{k-1}$ with $Z_{0}=-\ln
2$. It is easy to prove that the solution $Z_{k}$ is given by:

$Z_{k}=O_{k}-\frac{1}{2}H_{k}-\ln 2\ \ \Rightarrow\ \
I_{k}=\frac{\pi}{2}\frac{1}{2^{2k}}\binom{2k}{k}(O_{k}-\frac{1}{2}H_{k}-\ln
2).$

Hence, the right-hand side of (13) reduces to:

$\frac{\pi}{2}\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}(O_{k}-\frac{1}{2}H_{k}-\ln
2)\frac{1}{2k+1}-\frac{\pi}{2}\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{(2k+1)^{2}}.$

Multiplying both sides of (13) by $\frac{4}{\pi}$ and using (10), (11) and
(14) we obtain and equivalent form of (9).

So, it remains to evaluate the following series:

$\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{H_{k}}{2k+1}.$
(15)

Using a reindexing argument, we obtain the following:

$\displaystyle\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{H_{k}}{2k+1}=\sum_{k=1}^{\infty}\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k-1}\binom{2(k-1)}{k-1}^{2}\frac{H_{k-1}}{2k-1}$
$\displaystyle=-\frac{4(2G-1)}{\pi}+\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}H_{k}\left(\frac{1}{(2k-1)^{3}}+\frac{2}{(2k-1)^{2}}+\frac{1}{2k-1}\right).$

Using previously known Ramanujan-like series introduced in [4, 10], we obtain
that the series in (15) equals:

$\frac{24}{\pi}-\frac{8G}{\pi}-\frac{24\ln(2)}{\pi}+\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{H_{k}}{(2k-1)^{3}}.$

We claim that

$\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{H_{k}}{(2k-1)^{3}}$
(16)

equals

$\frac{8G}{\pi}-\frac{16G\ln(2)}{\pi}-\frac{16\mathcal{G}}{\pi}+\frac{5\pi^{2}}{8}-\frac{24}{\pi}+\frac{\ln^{2}(2)}{2}+\frac{24\ln(2)}{\pi}.$

To show this, we use an interated integral-based approach, using results from
[17]. We recall that for any real numbers $a,b$ and $1$-forms $f_{j}(t)\,dt$
($j=1,\dots,d$) with $d>1$ we may recursively define the iterated integral

$\int_{a}^{b}f_{1}(t)\,dt\,f_{2}(t)\,dt\cdots
f_{d}(t)\,dt:=\int_{a}^{b}f_{1}(u)\Big{(}\int_{a}^{u}f_{2}(t)\,dt\cdots
f_{d}(t)\,dt\Big{)}\,du.$

By the proof of [17, Theorem 6.1] we may use iterated integrals to compute
(16) as follows:

$\displaystyle I_{s}:=$
$\displaystyle\,\sum_{k>m>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{(2k-1)^{s+2}(2m)}$
$\displaystyle=$ $\displaystyle\,\frac{2}{\pi}\int_{0}^{\pi/2}\sin t\,dt(\cot
t\,dt)^{s}\,\cot^{2}t\,dt\,d(\sec t)\Big{(}\csc t-\cot t\Big{)}\,dt$
$\displaystyle=$ $\displaystyle\,\frac{2}{\pi}\int_{0}^{\pi/2}\sin t\,dt(\cot
t\,dt)^{s}\,(\csc^{2}t-1)\,dt\,d(\sec t)\Big{(}\csc t-\cot
t\Big{)}\,dt=X_{s}-Y_{s}-I_{s-1},$

where

$\displaystyle X_{s}:=$ $\displaystyle\,\frac{2}{\pi}\int_{0}^{\pi/2}\sin
t\,dt(\cot t\,dt)^{s}\,\sec t\,dt\,\Big{(}\csc t-\cot t\Big{)}\,dt,$
$\displaystyle Y_{s}:=$ $\displaystyle\,\frac{2}{\pi}\int_{0}^{\pi/2}\sin
t\,dt(\cot t\,dt)^{s}\,dt\,d(\sec t)\Big{(}\csc t-\cot t\Big{)}\,dt.$

Set 1-forms ${\texttt{a}}=dt/t$, ${\texttt{x}}_{\mu}=dt/(\mu-t)$,
${\texttt{d}}_{\mu,\nu}={\texttt{x}}_{\mu}-{\texttt{x}}_{\nu}$ for any roots
of unity $\mu$ and $\nu$, and
${\texttt{y}}={\texttt{x}}_{i}+{\texttt{x}}_{-i}-{\texttt{x}}_{1}-{\texttt{x}}_{-1}$.
Then

$\displaystyle X_{0}=$
$\displaystyle\,\frac{2}{\pi}\int_{0}^{\pi/2}dt\,\Big{(}\csc t-\cot
t\Big{)}\,dt=\frac{2}{\pi}\int_{0}^{1}i\Big{(}2{\texttt{x}}_{-1}-{\texttt{x}}_{i}-{\texttt{x}}_{-i}\Big{)}{\texttt{d}}_{-i,i},$
$\displaystyle X_{1}=$ $\displaystyle\,\frac{2}{\pi}\int_{0}^{\pi/2}(\cos
t\cot t)\,dt\,\sec t\,dt\,\Big{(}\csc t-\cot t\Big{)}\,dt$ $\displaystyle=$
$\displaystyle\,\frac{2}{\pi}\int_{0}^{\pi/2}(\csc t-\sin t)\,dt\,\sec
t\,dt\,\Big{(}\csc t-\cot t\Big{)}\,dt$ $\displaystyle=$
$\displaystyle\,-X_{0}+\frac{2}{\pi}\int_{0}^{\pi/2}\csc t\,dt\,\sec
t\,dt\,\Big{(}\csc t-\cot t\Big{)}\,dt$ $\displaystyle=$
$\displaystyle\,-X_{0}-(-1)^{s}\frac{2}{\pi}\int_{0}^{1}\Big{(}2{\texttt{x}}_{-1}-{\texttt{x}}_{i}-{\texttt{x}}_{-i}\Big{)}{\texttt{a}}{\texttt{d}}_{-1,1}$

by the change of variables $t\to\sin^{-1}(1-t^{2})/(1+t^{2})$ at the last step
for both $X_{0}$ and $X_{1}$. Similarly,

$\displaystyle Y_{0}:=$ $\displaystyle\,\frac{2}{\pi}\int_{0}^{\pi/2}\sin
t\,dt\,dt\,d(\sec t)\Big{(}\csc t-\cot t\Big{)}\,dt$ $\displaystyle=$
$\displaystyle\,\frac{2}{\pi}\int_{0}^{\pi/2}\cos t\,dt\,d(\sec t)\Big{(}\csc
t-\cot t\Big{)}\,dt$ $\displaystyle=$
$\displaystyle\,\frac{2}{\pi}\int_{0}^{\pi/2}dt\Big{(}\csc t-\cot
t\Big{)}\,dt-\frac{2}{\pi}\int_{0}^{\pi/2}\Big{(}1-\sin t\Big{)}\Big{(}\csc
t\sec t-\csc t\Big{)}\,dt$ $\displaystyle=$
$\displaystyle\,\frac{2}{\pi}\int_{0}^{1}i\Big{(}2{\texttt{x}}_{-1}-{\texttt{x}}_{i}-{\texttt{x}}_{-i}\Big{)}{\texttt{d}}_{-i,i}+i{\texttt{d}}_{-i,i}-2{\texttt{x}}_{-1},$
$\displaystyle Y_{1}:=$ $\displaystyle\,\frac{2}{\pi}\int_{0}^{\pi/2}\sin
t\,dt\cot t\,dt\,dt\,d(\sec t)\Big{(}\csc t-\cot t\Big{)}\,dt$
$\displaystyle=$ $\displaystyle\,\frac{2}{\pi}\int_{0}^{\pi/2}(\cos t\cot
t)\,dt\,dt\,d(\sec t)\Big{(}\csc t-\cot t\Big{)}\,dt$ $\displaystyle=$
$\displaystyle\,-Y_{0}+\frac{2}{\pi}\int_{0}^{\pi/2}\csc t\,dt\,dt\,d(\sec
t)\Big{(}\csc t-\cot t\Big{)}\,dt$ $\displaystyle=$
$\displaystyle\,-Y_{0}+\frac{2}{\pi}\int_{0}^{\pi/2}\csc t\,dt\Big{[}\sec
t\,dt\,\Big{(}\csc t-\cot t\Big{)}\,dt-dt\,\Big{(}\csc t\sec t-\csc
t\Big{)}\,dt\Big{]}$ $\displaystyle=$
$\displaystyle\,-Y_{0}+\frac{2}{\pi}\int_{0}^{1}\Big{(}2{\texttt{x}}_{-1}-{\texttt{x}}_{i}-{\texttt{x}}_{-i}\Big{)}{\texttt{a}}{\texttt{d}}_{-1,1}-\frac{2}{\pi}\int_{0}^{\pi/2}i\Big{(}{\texttt{a}}+2{\texttt{x}}_{-1}\Big{)}{\texttt{d}}_{-i,i}{\texttt{d}}_{-1,1}.$

Thus

$\displaystyle I_{1}=X_{1}-Y_{1}-I_{0}=$
$\displaystyle\,Y_{0}-X_{0}-I_{0}+\frac{2}{\pi}\int_{0}^{1}\Big{(}2{\texttt{x}}_{-1}-{\texttt{x}}_{i}-{\texttt{x}}_{-i}\Big{)}{\texttt{a}}{\texttt{d}}_{-1,1}$
$\displaystyle\,-\frac{2}{\pi}\int_{0}^{1}\Big{(}2{\texttt{x}}_{-1}-{\texttt{x}}_{i}-{\texttt{x}}_{-i}\Big{)}{\texttt{a}}{\texttt{d}}_{-1,1}+\frac{2}{\pi}\int_{0}^{\pi/2}i\Big{(}{\texttt{a}}+2{\texttt{x}}_{-1}\Big{)}{\texttt{d}}_{-i,i}{\texttt{d}}_{-1,1}$
$\displaystyle=$
$\displaystyle\,-I_{0}+\frac{2}{\pi}\int_{0}^{1}i\Big{(}2{\texttt{x}}_{-1}-{\texttt{x}}_{i}-{\texttt{x}}_{-i}\Big{)}{\texttt{d}}_{-i,i}+i{\texttt{d}}_{-i,i}-2{\texttt{x}}_{-1}$
$\displaystyle\,-\frac{2}{\pi}\int_{0}^{1}i\Big{(}2{\texttt{x}}_{-1}-{\texttt{x}}_{i}-{\texttt{x}}_{-i}\Big{)}{\texttt{d}}_{-i,i}+\frac{2}{\pi}\int_{0}^{\pi/2}i\Big{(}{\texttt{a}}+2{\texttt{x}}_{-1}\Big{)}{\texttt{d}}_{-i,i}{\texttt{d}}_{-1,1}.$

Note that $I_{0}$ is given by [17, Example B.8]. Hence we obtain the following
formula by using Au’s Mathematica package of colored multiple zeta values [2]:

$\displaystyle\,\sum_{k>m>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{(2k-1)^{3}(2m)}$
$\displaystyle=$
$\displaystyle\,\frac{2}{\pi}\bigg{(}\frac{5}{32}\pi^{3}-4\mathcal{G}-2G(1+2\ln
2)+6\ln 2+\frac{\pi}{8}\Big{(}8\ln 2+\ln^{2}2-12\Big{)}\bigg{)}.$

Now, using partial fraction decomposition, we see that

$\displaystyle\,\sum_{k>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{(2k-1)^{3}(2k)}$
$\displaystyle=$
$\displaystyle\,\sum_{k>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\left(\frac{1}{(2k-1)^{3}}-\frac{1}{(2k-1)^{2}}+\frac{1}{2k-1}-\frac{1}{2k}\right)$
$\displaystyle=$
$\displaystyle\,\frac{2}{\pi}\left(\frac{\pi}{2}+2G-3\right)-\frac{2}{\pi}\Big{(}2-\frac{\pi}{2}\Big{)}+\frac{2}{\pi}\Big{(}\frac{\pi}{2}-1\Big{)}-\frac{2}{\pi}(\pi\ln
2-2G)$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\Big{(}\frac{3\pi}{2}+4G-\pi\ln 2-6\Big{)}$

by Examples B.1 and B.3 of [17]. Therefore,

$\displaystyle\,\sum_{k>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{H_{k}}{(2k-1)^{3}}$
$\displaystyle=$
$\displaystyle\,2\sum_{k>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{(2k-1)^{3}(2k)}+2\sum_{k>m>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{(2k-1)^{3}(2m)}$
$\displaystyle=$
$\displaystyle\,\frac{4}{\pi}\bigg{(}\frac{5}{32}\pi^{3}-4\mathcal{G}+2G(1-2\ln
2)+6\ln 2+\frac{\pi}{8}\ln^{2}2-6\bigg{)}.$

So, we obtain that:

$\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{H_{k}}{2k+1}=-\frac{16G\ln(2)}{\pi}-\frac{16\mathcal{G}}{\pi}+\frac{5\pi^{2}}{8}+\frac{\ln^{2}(2)}{2}.$
(17)

Thus, according to (9), we find that:

$\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}}{2k+1}=-\frac{4G\ln(2)}{\pi}-\frac{16\mathcal{G}}{\pi}+\frac{3\pi^{2}}{8}+\frac{\ln^{2}(2)}{2}.$
(18)

So, according to our relation for (8), as derived from [6], we obtain the
desired result. ∎

###  2.2 Alternate approaches

The series in (3) may be reduced to a $\mathbb{Q}$-combination of the series
treated in [17]. We offer a sketch of a proof, as below, based on this
alternate approach.

First, for all $k\in\mathbb{N}$ and
${\boldsymbol{\sl{s}}}=(s_{1},\dots,s_{d})\in\mathbb{N}^{d}$ we define the
multiple $t$-sums as

$\displaystyle t_{k}({\boldsymbol{\sl{s}}}):=\sum_{k\geq
k_{1}>\dots>k_{d}>0}\frac{1}{(2k_{1}-1)^{s_{1}}\cdots(2k_{d}-1)^{s_{d}}}.$

They satisfy the stuffle relations. For example,

$\displaystyle t_{k}(1)^{2}=2t_{k}(1,1)+t_{k}(2).$

Since $O_{k}=t_{k}(1)$ the ${\lambda}=1$ case of (2) is

$\displaystyle\sum_{k>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}^{2}}{(2k-1)^{2}}=2\sum_{k>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{t_{k}(1,1)}{(2k-1)^{2}}+\sum_{k>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{t_{k}(2)}{(2k-1)^{2}}.$

Now, for all ${\boldsymbol{\sl{s}}}=(s_{1},\dots,s_{d})\in\mathbb{N}^{d}$ and
$m\in\mathbb{N}$,

$\displaystyle\,\sum_{k>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{t_{k}({\boldsymbol{\sl{s}}})}{(2k-1)^{m}}$
$\displaystyle=$
$\displaystyle\,\sum_{k>k_{2}>\cdots>k_{d}>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{(2k-1)^{s_{1}+m}(2k_{2}-1)^{s_{2}}\cdots(2k_{d}-1)^{s_{d}}}$
$\displaystyle\,+\sum_{k>k_{1}>\cdots>k_{d}>0}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{(2k-1)^{m}(2k_{1}-1)^{s_{1}}\cdots(2k_{d}-1)^{s_{d}}},$

and this can be handled by the approach in [17]. It is not hard to see that
this above procedure can be generalized to treat arbitrary powers of $O_{k}$,
not only its square.

In an unpublished online posting222https://mathoverflow.net/questions/323287,
Martin Nicholson described a technique, using results from [9], that may be
used to evaluate the series shown in (15) and (18). So, in view of the WZ-
derived identity involved in our proof of Theorem 1, along with the rest of
this proof, we may formulate an alternative proof of Theorem 1, using
Nicholson’s arguments to evaluate the series in (15) and (18).

##  3 The general case

Starting from the base case of (3) with $\lambda=1$, we are now able, using a
reindexing argument, to derive an evaluation for

$\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}^{2}}{(2k-3)^{2}}$
(19)

from Theorem 1. In the process we obtain an evaluation of

$\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}}{(2k-1)^{3}}.$
(20)

###  3.1 The $\lambda=2$ case

######  Theorem 2.

The series

$\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}^{2}}{(2k-3)^{2}}$

admits the following symbolic form:

$-\frac{32G}{27\pi}-\frac{32G\ln(2)}{9\pi}-\frac{64\mathcal{G}}{3\pi}+\frac{\pi^{2}}{2}+\frac{11\pi}{162}+\frac{16}{27\pi}+\frac{44\ln^{2}(2)}{27\pi}+\frac{2\ln^{2}(2)}{3}-\frac{52\ln(2)}{27\pi}.$

###### Proof.

Applying a reindexing argument to the $\lambda=1$ case of (3), we find that

$\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}O_{k-1}^{2}\left(\frac{9}{4(2k-3)^{2}}-\frac{3}{4(2k-3)}+\frac{1}{4(2k-1)^{2}}+\frac{3}{4(2k-1)}\right)$

admits the same symbolic form as in the $\lambda=1$ case. Using the
${}_{4}F_{3}(1)$-evaluation

${}_{4}F_{3}\\!\\!\left[\begin{matrix}-\frac{1}{2},-\frac{1}{2},-\frac{1}{2},-\frac{1}{2}\vspace{1mm}\\\
\frac{1}{2},\frac{1}{2},1\end{matrix}\ \Bigg{|}\
1\right]=-\frac{8G}{\pi}-\frac{16\mathcal{G}}{\pi}+\frac{3\pi^{2}}{8}+\frac{8}{\pi}+\frac{\ln^{2}(2)}{2},$
(21)

which was given in an equivalent form in [11], and which may be proved using a
reindexing argument from [11] together with a Fourier–Legendre-based argument
from [9], we obtain the following, using results from [4, 5, 6]: The
expression

$-\frac{1}{2}\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}}{(2k-1)^{3}}+\frac{9}{4}\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}^{2}}{(2k-3)^{2}}$
(22)

admits the symbolic form

$\frac{4G}{3\pi}-\frac{6G\ln(2)}{\pi}-\frac{32\mathcal{G}}{\pi}+\frac{3\pi^{2}}{4}+\frac{11\pi}{72}+\frac{4}{3\pi}+\frac{11\ln^{2}(2)}{3\pi}+\ln^{2}(2)-\frac{22\ln(2)}{3\pi}.$

We proceed to reindex (18):

$\displaystyle\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}}{2k+1}=\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k+1}-\frac{1}{2k+1}}{2k+1}$
$\displaystyle=\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}O_{k}\left(\frac{1}{2k-1}+\frac{2}{(2k-1)^{2}}+\frac{1}{(2k-1)^{3}}\right)-\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{(2k+1)^{2}}$
$\displaystyle=\frac{2\ln(2)}{\pi}+\frac{8(G-\ln(2))}{\pi}+\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}}{(2k-1)^{3}}-\left(-\frac{16\mathcal{G}}{\pi}+\frac{3\pi^{2}}{8}+\frac{\ln^{2}(2)}{2}\right).$

The first two series evaluations are from [16], and the sum of the last series
is given in (28).

In this manner, we obtain:

$\sum_{k=1}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{O_{k}}{(2k-1)^{3}}=-\frac{4G(\ln(2)+2)}{\pi}-\frac{32\mathcal{G}}{\pi}+\frac{6\ln(2)}{\pi}+\frac{3\pi^{2}}{4}+\ln^{2}(2).$
(23)

So, in view of our symbolic form for (22), we obtain the desired result. ∎

###  3.2 Further reindexings

We will need the following Lemma:

######  Lemma 1.

For $m\in\mathbb{N}_{0}$, we have that

$\sum_{n=0}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{1}{2n+1-2m}=-\frac{16^{m}}{2\pi
m^{2}\binom{2m}{m}^{2}}\sum_{k=0}^{m-1}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}=:A(m)$
(24)
$\displaystyle\sum_{n=0}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{1}{(2n+1-2m)^{2}}=$
$\displaystyle\frac{16^{m}}{m^{2}\binom{2m}{m}^{2}}\sum_{k=1}^{m-1}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{k}{2k+1}A(k)$
(25) $\displaystyle+\frac{16^{m}}{\pi
m^{2}\binom{2m}{m}^{2}}\sum_{k=0}^{m-1}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{2k+1}.$

The proof of these two results is similar to that of (2). See also [11].

Note that reindexing proves to be a very powerful tool. Using an alternative
method of reindexing we are able to deduce recurrences for the series (3),
(4), (5). For instance, for (3) we start from

$\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{(1+2n-2\lambda)^{2}}=\sum_{n=0}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{\left(O_{n+1}-\frac{1}{2n+1}\right)^{2}}{(1+2n-2\lambda)^{2}}.$

Expanding the square in the numerator gives us three terms. The third one,
which contains no factor $O_{n}$, can be split using a partial fraction
expansion:

$\frac{1}{(2n+1)^{2}(1+2n-2\lambda)^{2}}=\frac{1}{4\lambda^{3}}\left(\frac{1}{2n+1}-\frac{1}{2n+1-2\lambda}+\frac{\lambda}{(2n+1)^{2}}+\frac{\lambda}{(2n+1-2\lambda)^{2}}\right).$

This leads to four series for which the values are well known [11].

For the first and second terms we use reindexing:

$\displaystyle\sum_{n=0}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n+1}^{2}}{(1+2n-2\lambda)^{2}}$
$\displaystyle=\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n-1}\binom{2n-2}{n-1}^{2}\frac{O_{n}^{2}}{(1+2n-2(\lambda+1))^{2}}$
$\displaystyle=\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{4n^{2}}{(2n-1)^{2}(1+2n-2(\lambda+1))^{2}}O_{n}^{2},$
$\displaystyle-2\sum_{n=0}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}$
$\displaystyle\frac{O_{n+1}}{(2n+1)(1+2n-2\lambda)^{2}}$
$\displaystyle=-2\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{4n^{2}}{(2n-1)^{3}(1+2n-2(\lambda+1))^{2}}O_{n}.$

Now using partial fraction expansions both these series can be written as
linear combinations of series of the form (3), (4), (5), and some other series
which can be found in [16]. For the last one we also need (23). The same
method may be applied to (4) and (5). Rearranging the terms in the results
leads to the following three recurrences:
For the series (4):

$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}$
$\displaystyle\frac{O_{n}}{(1+2n-2(\lambda+1))^{2}}=\frac{4\lambda^{2}}{(1+2\lambda)^{2}}\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}}{(1+2n-2\lambda)^{2}}$
$\displaystyle+\frac{1}{\lambda(1+2\lambda)}\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}}{1+2n-2(\lambda+1)}$
$\displaystyle-\frac{3+2\lambda}{4\lambda^{2}(1+2\lambda)}\sum_{n=0}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{1}{1+2n-2(\lambda+1)}$
$\displaystyle+\frac{4\lambda^{2}}{(1+2\lambda)^{2}}\sum_{n=0}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{1}{(1+2n-2(\lambda+1))^{2}}$
$\displaystyle-\frac{4\ln(2)\lambda+3}{2\lambda^{2}(1+2\lambda)^{2}\pi}.$

The value of the second series at the RHS can be found in [16],  the sums of
the third and fourth series can be found using Lemma 1.

The first few values are given by:

$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}}{(2n-3)^{2}}=-\frac{44\ln\left(2\right)}{27\pi}+\frac{26}{27\pi}+\frac{16G}{9\pi},$
$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}}{(2n-5)^{2}}=-\frac{3388\ln\\!\left(2\right)}{3375\pi}+\frac{2954}{3375\pi}+\frac{256G}{225\pi},$
$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}}{(2n-7)^{2}}=-\frac{92804\ln\\!\left(2\right)}{128625\pi}+\frac{32642}{42875\pi}+\frac{1024G}{1225\pi},$

confirming the results obtained in [16] using Mathematica.
For the series (5):

$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}$
$\displaystyle\frac{O_{n}^{2}}{1+2n-2(\lambda+1)}=\frac{4\lambda^{2}}{(1+2\lambda)^{2}}\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{1+2n-2\lambda}$
$\displaystyle+\frac{1}{\lambda}\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}}{1+2n-2(\lambda+1)}$
$\displaystyle-\frac{1}{4\lambda^{2}}\sum_{n=0}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{1}{1+2n-2(\lambda+1)}$
$\displaystyle-\frac{(24\ln^{2}(2)+\pi^{2})\lambda^{2}+24\ln(2)\lambda+6}{12\lambda^{2}(1+2\lambda)^{2}\pi}.$

Using the case $\lambda=1$ given in (7), we obtain

$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{2n-3}=-\frac{10\ln^{2}(2)}{9\pi}-\frac{5\pi}{108}+\frac{8\ln\\!\left(2\right)}{9\pi}+\frac{16G}{9\pi}-\frac{2}{9\pi},$
$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{2n-5}=-\frac{178\ln^{2}(2)}{225\pi}-\frac{89\pi}{2700}+\frac{208\ln\\!\left(2\right)}{225\pi}+\frac{256G}{225\pi}-\frac{74}{225\pi},$
$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{2n-7}=-\frac{762\ln^{2}(2)}{1225\pi}-\frac{127\pi}{4900}+\frac{3208\ln\\!\left(2\right)}{3675\pi}+\frac{1024G}{1225\pi}-\frac{818}{2205\pi}.$

For the series (3):

$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}$
$\displaystyle\frac{O_{n}^{2}}{(1+2n-2(\lambda+1))^{2}}=\frac{4\lambda^{2}}{(1+2\lambda)^{2}}\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{(1+2n-2\lambda)^{2}}$
$\displaystyle+\frac{1}{\lambda(1+2\lambda)}\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{1+2n-2(\lambda+1)}$
$\displaystyle-\frac{3+2\lambda}{2\lambda^{2}(1+2\lambda)}\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}}{1+2n-2(\lambda+1)}$
$\displaystyle+\frac{1}{\lambda}\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}}{(1+2n-2(\lambda+1))^{2}}$
$\displaystyle+\frac{1+\lambda}{2\lambda^{3}(1+2\lambda)}\sum_{n=0}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{1}{1+2n-2(\lambda+1)}$
$\displaystyle-\frac{1}{4\lambda^{2}}\sum_{n=0}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{1}{(1+2n-2(\lambda+1))^{2}}$
$\displaystyle+\frac{(24\ln^{2}(2)+\pi^{2})\lambda^{2}+36\ln(2)\lambda+12}{12\lambda^{3}(1+2\lambda)^{2}\pi}.$

Note that for this last one we need the results from the previous two
recurrences.

The case $\lambda=2$ is discussed in Theorem 2. For $\lambda=3$ and
$\lambda=4$ we obtain:

$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{(2n-5)^{2}}=\frac{32\ln^{2}(2)}{75}+\frac{8\pi^{2}}{25}+\frac{3388\ln^{2}(2)}{3375\pi}-\frac{512G\ln\\!\left(2\right)}{225\pi}$
$\displaystyle\hskip
142.26378pt+\frac{847\pi}{20250}-\frac{5908\ln\\!\left(2\right)}{3375\pi}-\frac{256G}{3375\pi}-\frac{1024\mathcal{G}}{75\pi}+\frac{2624}{3375\pi},$
$\displaystyle\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{(2n-7)^{2}}=\frac{384\ln^{2}(2)}{1225}+\frac{288\pi^{2}}{1225}+\frac{92804\ln^{2}(2)}{128625\pi}-\frac{2048G\ln\\!\left(2\right)}{1225\pi}$
$\displaystyle\hskip
85.35826pt+\frac{23201\pi}{771750}-\frac{65284\ln\\!\left(2\right)}{42875\pi}+\frac{11264G}{42875\pi}-\frac{12288\mathcal{G}}{1225\pi}+\frac{188144}{231525\pi}.$

##  4 Further evaluations

Although the main purpose of our work is devoted to Wang and Chu’s problem
concerning (3), we remark that many of the formulas involved in our proof of
Theorem 1 are of interest in their own right and may be applied to evaluate
intractable integrals over elementary expressions. The formula in (17), in
particular, together with the integration technique from [4], gives us that

$\int_{0}^{1}\frac{\sin^{-1}(x)\ln\left(1-x^{2}\right)}{x\sqrt{1-x^{2}}}\,dx=4G\ln(2)+8\mathcal{G}-\frac{5\pi^{3}}{16}-\frac{\pi\ln^{2}(2)}{4},$

but current computer algebra systems, including the 2022 version of
Mathematica, cannot evaluate this definite integral and cannot evaluate the
underlying indefinite integral. Note that

$\int_{0}^{1}\frac{\sin^{-1}(x)\ln\left(1-x^{2}\right)}{x\sqrt{1-x^{2}}}\,dx=2\int_{0}^{\pi/2}\frac{\theta}{\sin\theta}\ln(\cos\theta)d\theta.$
(26)

Similarly, using (18) together with the integration method from [5], we may
obtain the following remarkable evaluation:

$\displaystyle\int_{0}^{1}\frac{\sqrt{1-x^{2}}\sin^{-1}(x)\ln(x)}{x}\,dx$
$\displaystyle{\color[rgb]{0,0,0}=\int_{0}^{\pi/2}\frac{\theta\cos^{2}\theta}{\sin\theta}\ln(\sin\theta)d\theta}$
$\displaystyle=-2G-4\mathcal{G}+\frac{\pi^{3}}{32}+2+\frac{\pi\ln^{2}(2)}{8},$

which together with (26) leads to

$G=1+\frac{1}{2}\int_{0}^{\pi/2}\theta\,{\sin\theta}\ln(\sin\theta)d\theta.$

Much of our proof in Section 2.1 relied on the evaluation of series involving
summand factors of the form

$\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{2k+1}$

for $k\in\mathbb{N}_{0}$, so it seems worthwhile to consider further
evaluations for such sums and relation summations. In this regard, setting
$\beta(4)=\sum_{n=0}^{\infty}\frac{(-1)^{n}}{(2n+1)^{4}}$, let us reproduce
the following identities from [8, 9]:

$\displaystyle\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{2k+1}=\frac{4G}{\pi},$
(27)
$\displaystyle\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{(2k+1)^{2}}=-\frac{16\mathcal{G}}{\pi}+\frac{3\pi^{2}}{8}+\frac{\ln^{2}(2)}{2},$
(28)
$\displaystyle\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}\frac{1}{(2k+1)^{3}}=-64\frac{\Im\left(\text{Li}_{4}(1-i)\right)}{\pi}-\pi^{2}\ln(2)-\frac{2}{3}\ln^{3}(2)-48\frac{\beta(4)}{\pi}.$
(29)

Such evaluations are closely related to the evaluation of series of the form
$\sum_{k=1}^{\infty}\frac{\binom{2k}{k}^{2}}{k^{n}\,16^{k}}$. For example,
consider the following analogue of (28).

######  Theorem 3.

The evaluation

${\sum_{k=1}^{\infty}\frac{\binom{2k}{k}^{2}}{k^{2}\,16^{k}}=\frac{3\pi^{2}}{2}-\frac{64\mathcal{G}}{\pi}-6\ln^{2}(2)}$

holds true.

###### Proof.

For our proof, we have to calculate the following:

$\sum_{k=1}^{\infty}\frac{\binom{2k}{k}^{2}}{k^{2}\,16^{k}}\cdot\frac{\pi}{2}=\int_{0}^{\pi/2}2\,{\rm
Li}_{2}(\tfrac{1-\cos(t)}{2}){dt}-\int_{0}^{\pi/2}\ln^{2}(\cos(\tfrac{t}{2})){dt}.$
(30)

The second integral in (30) is equivalent with

$\int_{0}^{\pi/4}\ln^{2}(\cos(\tfrac{t}{2})){dt}=-\frac{7\pi^{3}}{64}+\frac{7}{32}\pi\ln^{2}(2)+\mathcal{G}-\frac{G\ln(2)}{2}$
(31)

which can be found in [8]. For the first one, we use partial integration:

$\displaystyle\int_{0}^{\pi/2}2\,{\rm Li}_{2}(\tfrac{1-\cos(t)}{2}){dt}$
$\displaystyle=\frac{1}{12}\pi^{3}-\frac{1}{2}\pi\ln^{2}(2)+2\int_{0}^{\pi/2}t\cot(\tfrac{t}{2})\ln(\cos^{2}(\tfrac{t}{2})){dt}$
$\displaystyle=\frac{1}{12}\pi^{3}-\frac{1}{2}\pi\ln^{2}(2)+8\int_{0}^{\pi/4}t\cot(t)\ln(\cos^{2}(t)){dt}.$

For this integral we use the substitution $u=\tan(t)$:

$\displaystyle\int_{0}^{\pi/4}t\cot(t)\ln(\cos^{2}(t)){dt}$
$\displaystyle=-\int_{0}^{1}\arctan(u)\ln(1+u^{2})\frac{1}{u(1+u^{2})}{du}$
$\displaystyle=-\int_{0}^{1}\arctan(u)\ln(1+u^{2})\left(\frac{1}{u}-\frac{u}{1+u^{2}}\right){du}.$

The first integral can be found in [8].The second one requires one partial
integration in combination with [8]:

$\displaystyle\int_{0}^{\pi/4}\ln^{2}(\sin(t))dt=-\mathcal{G}+\frac{1}{2}G\ln(2)+\frac{23}{384}\pi^{3}+\frac{9}{32}\pi\ln^{2}(2)$

and (11). Bringing everything together, we obtain the desired result. ∎

######  Remark 1.

Another method to evaluate the series in Theorem 3 is given in [17, Example
B.1] using different basis elements. Moreover, Xu and Zhao proved in loc. cit.
that

$\sum_{k=1}^{\infty}\frac{\binom{2k}{k}^{2}}{k^{3}\,16^{k}}=-\frac{512}{\pi}\Im\left(\text{Li}_{4}\Big{(}\frac{1+i}{2}\Big{)}\right)-6\pi^{2}\ln(2)+8\ln^{3}(2)+384\frac{\beta(4)}{\pi}+4\zeta(3).$
(32)

Noting that

$\sum_{k=1}^{\infty}\frac{\binom{2k}{k}^{2}}{k^{n}\,16^{k}}=\sum_{k=0}^{\infty}\frac{4(2k+1)^{2}\binom{2k}{k}^{2}}{(k+1)^{n+2}\,16^{k+1}}=\sum_{k=0}^{\infty}\frac{4(2(k+1)-1)^{2}\binom{2k}{k}^{2}}{(k+1)^{n+2}\,16^{k+1}},$

we obtain, using previous results, that

$\displaystyle\sum_{k=0}^{\infty}\frac{\binom{2k}{k}^{2}}{(k+1)^{3}\,16^{k}}$
$\displaystyle=\frac{48}{\pi}-16-\frac{32G}{\pi}-16\ln(2),$
$\displaystyle\sum_{k=0}^{\infty}\frac{\binom{2k}{k}^{2}}{(k+1)^{4}\,16^{k}}$
$\displaystyle=\frac{128}{\pi}-48-\frac{128G}{\pi}+64\ln(2)+6\pi^{2}-\frac{256\mathcal{G}}{\pi}-24\ln^{2}(2),$
$\displaystyle\sum_{k=0}^{\infty}\frac{\binom{2k}{k}^{2}}{(k+1)^{5}\,16^{k}}$
$\displaystyle=\frac{320}{\pi}-128-\frac{384G}{\pi}+192\ln(2)+24\pi^{2}+\frac{1024\mathcal{G}}{\pi}-96\ln^{2}(2)$
$\displaystyle-\frac{2048}{\pi}\Im\left(\text{Li}_{4}\Big{(}\frac{1+i}{2}\Big{)}\right)-24\pi^{2}\ln(2)+32\ln^{3}(2)+\frac{1536\beta(4)}{\pi}+16\zeta(3).$

Our applications, as above, of series as in (27)–(29) motivate the study and
application of identities as below.

######  Theorem 4.

For any integer $p\geq 0$, we have

$\displaystyle\sum_{n=0}^{\infty}\binom{2n}{n}\frac{z^{2n+1}}{(2n+1)^{p+1}}=\frac{\theta}{2}\frac{\ln^{p}(2\sin\theta)}{p!}+\frac{1}{4p!}\sum_{j=1}^{p}(-1)^{j-1}\binom{p}{j}\ln^{p-j}(2\sin\theta){\rm
Ls}_{j+1}(2\theta),$ (33)

where $\theta:=\arcsin(2z)$ and

$\displaystyle{\rm
Ls}_{j}^{(k)}(\theta):=-\int_{0}^{\theta}t^{k}\ln^{j-k-1}\left(2\sin\frac{t}{2}\right)dt\quad\text{and}\quad{\rm
Ls}_{j}(\theta):={\rm Ls}_{j}^{(0)}(\theta).$ (34)

###### Proof.

From [3], we have

$\displaystyle\sum_{n=0}^{\infty}\binom{2n}{n}\frac{z^{2n+1}}{2n+1}=\frac{1}{2}\arcsin(2z).$
(35)

Hence, we get

$\displaystyle\sum_{n=0}^{\infty}\binom{2n}{n}\frac{z^{2n+1}}{(2n+1)^{p+1}}=\frac{1}{2(p-1)!}\int_{0}^{z}\frac{\ln^{p-1}\left(\frac{z}{t}\right)\arcsin(2t)}{t}dt$
$\displaystyle=\frac{1}{2(p-1)!}\int_{0}^{z}\frac{\ln^{p-1}\left(\frac{4z}{4t}\right)\arcsin(2t)}{4t}d(4t)$
$\displaystyle=\frac{1}{2(p-1)!}\sum_{j=0}^{p-1}\binom{p-1}{j}(-1)^{j}\ln^{p-1-j}(4z)\int_{0}^{z}\frac{\ln^{j}(4t)\arcsin(2t)}{4t}d(4t)$
$\displaystyle=\frac{1}{2(p-1)!}\sum_{j=0}^{p-1}\binom{p-1}{j}\frac{(-1)^{j}}{j+1}\ln^{p-1-j}(4z)\int_{0}^{z}\arcsin(2t)d\ln^{j+1}(4t)$
$\displaystyle=\frac{1}{2(p-1)!}\sum_{j=0}^{p-1}\binom{p-1}{j}\frac{(-1)^{j}}{j+1}\ln^{p-1-j}(4z)\int_{0}^{\theta}xd\ln^{j+1}(2\sin
x)$
$\displaystyle=\frac{1}{2(p-1)!}\sum_{j=0}^{p-1}\binom{p-1}{j}\frac{(-1)^{j}}{j+1}\ln^{p-1-j}(4z)\left(\theta\ln^{j+1}(2\sin\theta)-\int_{0}^{\theta}\ln^{j+1}(2\sin
x)dx\right)$
$\displaystyle=\frac{1}{2(p-1)!}\theta\ln^{p}(2\sin\theta)\sum_{j=0}^{p-1}\binom{p-1}{j}\frac{(-1)^{j}}{j+1}+\frac{1}{4p!}\sum_{j=1}^{p}(-1)^{j-1}\binom{p}{j}\ln^{p-j}(2\sin\theta){\rm
Ls}_{j+1}(2\theta).$

Then, noting the fact that

$\sum_{j=0}^{p-1}\binom{p-1}{j}\frac{(-1)^{j}}{j+1}=\frac{1}{p},$

we obtain the desired evaluation. ∎

Setting $p\leq 3$ in (33), we deduce (with $\theta:=\arcsin(2z)$)

$\displaystyle\sum_{n=0}^{\infty}\binom{2n}{n}\frac{z^{2n+1}}{2n+1}=\frac{1}{2}\theta,$
(36)
$\displaystyle\sum_{n=0}^{\infty}\binom{2n}{n}\frac{z^{2n+1}}{(2n+1)^{2}}=\frac{1}{2}\theta\ln(2\sin\theta)+\frac{1}{4}{\rm
Cl}_{2}(2\theta),$ (37)
$\displaystyle\sum_{n=0}^{\infty}\binom{2n}{n}\frac{z^{2n+1}}{(2n+1)^{3}}=\frac{1}{4}\theta\ln^{2}(2\sin\theta)+\frac{1}{4}\ln(2\sin\theta){\rm
Cl}_{2}(2\theta)-\frac{1}{8}{\rm Ls}_{3}(2\theta),$ (38)

where we used the relation ${\rm Ls}_{2}(\theta)={\rm Cl}_{2}(\theta)$ and for
any positive integer $m$ the Clausen function ${\rm Cl}_{m}(\theta)$ is
defined by

$\displaystyle{{\rm
Cl}}_{2m-1}(\theta):=\sum_{n=1}^{\infty}\frac{\cos(n\theta)}{n^{2m-1}}\quad\text{and}\quad{{\rm
Cl}}_{2m}(\theta):=\sum_{n=1}^{\infty}\frac{\sin(n\theta)}{n^{2m}}\qquad\forall\theta\in[0,\pi].$

If we replace $z$ by $\frac{\sin\theta}{2}$ in (36), divide by $\sin\theta$,
multiply by $\ln(\sin\theta)$ and integrate between 0 and $\frac{\pi}{2}$, we
get the following result:

$\displaystyle\int_{0}^{\pi/2}\frac{\theta}{\sin\theta}\ln\sin(\theta)\,{d\theta}$
$\displaystyle=\frac{\pi}{2}\sum_{k=0}^{\infty}\left(\frac{1}{16}\right)^{k}\binom{2k}{k}^{2}(O_{k}-\frac{1}{2}H_{k}-\ln
2)\frac{1}{2k+1}$
$\displaystyle=-4\mathcal{G}+\frac{\pi^{3}}{32}+\frac{\pi\ln^{2}(2)}{8}$

using (9). This integral is a companion to (26).

##  5 Some concluding remarks

In this paper, we have studied a few variations of Ramanujan-like series as in
(1), and we have provided a solution to a problem raised by Wang and Chu in
[16]. The main result gives closed evaluations of the series of the form

$\sum_{n=1}^{\infty}\left(\frac{1}{16}\right)^{n}\binom{2n}{n}^{2}\frac{O_{n}^{2}}{(1+2n-2\lambda)^{2}}$

for $\lambda=1,2,3,4$. Our method works for general $\lambda$ although the
computation will be more complicated as $\lambda$ increases. In the process,
we have utilized a combination of several different approaches to calculate
many other series and a few highly non-trivial definite integrals in the
intermediate steps, which are of interest in their own right.

Our last result, Theorem 4, differs from the others in flavor in two aspects:
First, the binomial coefficients appear as a first power instead of a square;
second, the series involves a variable $z$ so that its specializations may
offer many series identities involving log-sine integrals. We wonder if
Theorem 4 can be generalized to evaluate the following two series:

$\displaystyle\sum_{n=0}^{\infty}\binom{2n}{n}\frac{z^{2n+1}}{(2n+1)^{p}}H_{n}=?\quad\text{and}\quad\sum_{n=0}^{\infty}\binom{2n}{n}\frac{z^{2n+1}}{(2n+1)^{p}}H_{2n+1}=?\quad(p\in\mathbb{N}).$
(39)

These would lead to possible solutions of some conjectures of Z.-W. Sun [14].
For example, Prof. Sun conjectured the following two identities of Apéry-like
sums (see [14, Conjectures 10.59 and 10.60]):

$\displaystyle\sum_{n=0}^{\infty}\frac{\binom{2n}{n}}{(2n+1)^{2}(-16)^{n}}\left(5H_{2n+1}+\frac{12}{2n+1}\right)\overset{\text{?}}{=}14\zeta(3),$
(40)
$\displaystyle\sum_{n=0}^{\infty}\frac{\binom{2n}{n}}{(2n+1)^{3}16^{n}}\left(9H_{2n+1}+\frac{32}{2n+1}\right)\overset{\text{?}}{=}40\beta(4)+\frac{5}{12}\pi\zeta(3).$
(41)

### Competing interests statement

There are no competing interests to declare.

### Acknowledgement

Ce Xu is supported by the National Natural Science Foundation of China (Grant
No. 12101008), the Natural Science Foundation of Anhui Province (Grant No.
2108085QA01), and the University Natural Science Research Project of Anhui
Province. The corresponding Grant Number for the last case is KJ2020A0057. J.
Zhao is supported by the Jacobs Prize from The Bishop’s School.

## References

  * [1] V. S. Adamchik, A certain series associated with Catalan’s constant, _Z. Anal. Anwendungen_ 21 (2002), 817–826.
  * [2] K. C. Au, Evaluation of one-dimensional polylogarithmic integral, with applications to infinite series, arXiv:2007.03957. A companion Mathematica package available at researchgate.net/publication/357601353.
  * [3] B. Berndt, _Ramanujan’s Notebooks, Vol.2_. New York: Springer-Verlag, 1989.
  * [4] J. M. Campbell, Ramanujan-like series for $\frac{1}{\pi}$ involving harmonic numbers, _Ramanujan J._ 46 (2018), 373–387.
  * [5] J. M. Campbell, Series containing squared central binomial coefficients and alternating harmonic numbers, _Mediterr. J. Math._ 16 (2019), Paper No. 37, 7.
  * [6] J. M. Campbell, A Wilf-Zeilberger-based solution to the Basel problem with applications, _Discrete Math. Lett._ 10 (2022), 21–27.
  * [7] J. M. Campbell and K.-W. Chen, Explicit identities for infinite families of series involving squared binomial coefficients, _J. Math. Anal. Appl._ 513 (2022), Paper No. 126219, 23.
  * [8] J. M. Campbell, P. Levrie, and A. S. Nimbran, A natural companion to Catalan’s constant, _J. Class. Anal._ 18 (2021), 117–135.
  * [9] M. Cantarini and J. D’Aurizio, On the interplay between hypergeometric series, Fourier-Legendre expansions and Euler sums, _Boll. Unione Mat. Ital._ 12 (2019), 623–656.
  * [10] H. Chen, Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers, _J. Integer Seq._ 19 (2016), Article 16.1.5, 11.
  * [11] A. S. Nimbran, Deriving Forsyth-Glaisher type series for $\frac{1}{\pi}$ and Catalan’s constant by an elementary method, _Math. Student_ 84 (2015), 69–86.
  * [12] A. S. Nimbran and P. Levrie, Series of the form $\sum a_{n}\binom{2n}{n}$, submitted to Math. Student.
  * [13] M. Petkovšek, H. S. Wilf, and D. Zeilberger, $A=B$. A K Peters, Ltd., Wellesley, MA, 1996.
  * [14] Z.-W. Sun, _New Conjectures in Number Theory and Combinatorics_ , Harbin Institute of Technology Press, (in Chinese) 2021.
  * [15] X. Wang and W. Chu, Further Ramanujan-like series containing harmonic numbers and squared binomial coefficients, _Ramanujan J._ 52 (2020), 641–668.
  * [16] X. Wang and W. Chu, Series with harmonic-like numbers and squared binomial coefficients, _Rocky Mountain J. Math._ (2022).
  * [17] C. Xu and J. Zhao, Apéry-type series with summation indices of mixed parities and colored multiple zeta values, II, (2022). arXiv:2203.00777.