# ReactAIvate: A Deep Learning Approach to Predicting Reaction Mechanisms and
Unmasking Reactivity Hotspots

Ajnabiul Hoque0000-0001-9807-3061111Equal contribution. Manajit
Das0000-0001-7709-8809222Equal contribution. Mayank
Baranwal0000-0001-9354-2826 Corresponding Author. Email:
<EMAIL_ADDRESS>Raghavan B. Sunoj0000-0002-6484-2878 Corresponding
Author. Email<EMAIL_ADDRESS>Department of Chemistry, Indian Institute
of Technology Bombay, India  Department of Systems & Control Engineering,
Indian Institute of Technology, India  Tata Consultancy Services Research,
Mumbai, India  Centre for Machine Intelligence and Data Science, Indian
Institute of Technology Bombay, India

###### Abstract

A chemical reaction mechanism (CRM) is a sequence of molecular-level events
involving bond-breaking/forming processes, generating transient intermediates
along the reaction pathway as reactants transform into products. Understanding
such mechanisms is crucial for designing and discovering new reactions. One of
the currently available methods to probe CRMs is quantum mechanical (QM)
computations. The resource-intensive nature of QM methods and the scarcity of
mechanism-based datasets motivated us to develop reliable ML models for
predicting mechanisms. In this study, we created a comprehensive dataset with
seven distinct classes, each representing uniquely characterized elementary
steps. Subsequently, we developed an interpretable attention-based GNN that
achieved near-unity and 96% accuracy, respectively for reaction step
classification and the prediction of reactive atoms in each such step,
capturing interactions between the broader reaction context and local active
regions. The near-perfect classification enables accurate prediction of both
individual events and the entire CRM, mitigating potential drawbacks of
Seq2Seq approaches, where a wrongly predicted character leads to incoherent
CRM identification. In addition to interpretability, our model adeptly
identifies key atom(s) even from out-of-distribution classes. This
generalizabilty allows for the inclusion of new reaction types in a modular
fashion, thus will be of value to experts for understanding the reactivity of
new molecules.

1850

##  1 Introduction

The reliable prediction of chemical reactions holds paramount significance in
pharmaceutical and materials manufacturing, and in understanding many
processes in molecular biology [9, 5, 39]. A chemical reaction entails the
reorganization of atoms and bonds within the initial reactants, resulting in
the creation of novel molecules or compounds as the final products. To
comprehend a chemical reaction, it is essential to know the underlying
chemical transformations. A sequence of these transformation steps, also
called elementary mechanistic steps, is generally expressed in the form of a
chemical reaction mechanism (CRM). Many of these steps may involve a bond-
forming/bond-breaking process characterized by the corresponding transition
state. These elementary steps serve as building blocks for developing novel
reactions and discerning side products. Knowledge of CRM can provide atomic-
level insights into why the products are formed [19, 12].

---
Figure 1: (a) A representative example of oxidative addition template, (b)
complete workflow of our proposed ReactAIvate method, (c) process of reaction
step classification and reactive atom identification using ReactAIvate

One of the ways to identify a CRM is to perform quantum mechanical (QM)
calculations [6, 3]. Such methods are often computationally demanding and
often require substantial human attention rendering it a time-consuming task.
In recent years, numerous ML studies have focused on retrosynthesis and
forward synthesis prediction, with models mostly trained on the USPTO-50k
dataset [38]. While the dataset proves valuable for direct product prediction,
it lacks information about elementary reaction steps, thus lacking
opportunities to understand CRM [44, 10]. In the realm of CRM prediction,
navigating the complexities of multi-step reactions and ensuring atom and
charge balance presents formidable challenges. The conventional transformer-
based sequence-to-sequence (Seq2Seq) models, commonly employed for sequence
generation, prove inadequate in handling the intricate long-term dependencies
inherent in CRM [15]. Their limitations become glaring when even a single
incorrect character introduced during inference can render the entire CRM
meaningless. Beyond the conventional focus on atom and charge balance,
ensuring both semantic and syntactic validity of product SMILES strings –
representing molecules in a sequence-based format – highlights the necessity
for a comprehensive reevaluation at the modeling level in CRM prediction.
Recognizing these challenges, we have endeavored to craft an interpretable,
swift, and dependable alternative for CRM prediction.

Herein, we propose an interpretable attention-based graph neural network (GNN)
model for the elementary reaction step predictions, which are then used to
generate the full CRM (Figure 1). We introduce a dataset containing 3
different families of catalytic reactions comprising 7 distinct elementary
steps. In addition to the key task of identifying the elementary steps, our
model is simultaneously trained to detect the reactive atoms in such steps. In
the case of out-of-distribution (OOD) samples, our model points to the
reactive atoms, besides classifying them into an unseen group. This can act as
a guide to experts for understanding the reactivity. Below we summarize our
key contributions:

1\. CRM dataset: While there are several benchmark datasets, such as
USPTO-50k, comprising single-step chemical reactions, there is currently no
existing dataset on CRMs within the available literature as per our knowledge.
We curate a first-of-its-kind CRM dataset containing elementary mechanistic
steps for transition metal-catalyzed reactions.

2\. CRM identification via reaction step classification: We introduce
ReactAIvate, a graph-attention-based classification model that precisely
identifies the necessary elementary steps for a given combination of
reactants, reagents, and catalysts. Our model’s accurate intermediate-step
classification is pivotal for CRM identification. Notably, our approach
focuses on identifying underlying reaction rules rather than generating exact
product SMILES, simplifying the problem, as our experiments reveal that SMILES
generation can lead to vacuous CRM predictions due to even a single character
mismatch.

3\. Identification of reactive atoms/groups: Our framework is distinctly
trained to minimize a composite of two distinct loss types: (a) Graph-level
loss for predicting reaction classes, and (b) Node-level loss for
distinguishing reactive and non-reactive atoms. As a result, our curated
database includes information on reactive atoms for each specific reaction
step within a CRM. This feature is primarily introduced to offer valuable
insights to domain experts.

4\. Visualizing reactive centers via attention Mechanism: Expanding on the
previous point, we demonstrate that the inclusion of node-level loss
inherently compels the attention mechanism to align with reactive centers
within the molecules involved. This alignment significantly enhances the
visualization of reactive centers.

5\. Generalization to OOD samples: ReactAIvate avoids overconfidence by
introducing an OOD eighth class, enabling accurate identification of scenarios
not covered in training data. This enhances trust in predictions and allows
for the incorporation of new reaction classes when relevant data becomes
available. Attention visualization on OOD samples also reveals potential
reaction centers, showcasing the model’s adaptability.

##  2 Related Work

Recent years have witnessed several interesting applications of machine
learning in predicting molecular properties as well as their reactions [49, 2,
30]. Intriguing ML algorithms have been developed to make complex chemical
problems, such as organic synthesis, increasingly more amenable. In tasks such
as forward or retrosynthesis predictions, predefined reaction templates are
employed to make an intuitive connection from a set of reactants to
product(s). The templates can be obtained by using data-driven approaches [42]
or through encoding by domain experts. Authors in [13] combined these
templates and ML for product prediction and subsequent ranking of the
products. Transformation rules in retrosynthesis tasks were leveraged in [43].
Current trends suggest the use of templates for both forward and
retrosynthesis analysis owing to their efficiency and interpretability. [50,
41, 17, 20].

Interesting alternatives relying on template-free methods have emerged in very
recent times. Transformer-based Seq2Seq generative model has been adopted for
forward and retrosynthesis predictions using SMILES as the molecular
representation [22, 40, 27]. Various models have been proposed that differ in
molecular representation and/or model architecture [45, 53, 7, 23, 48]. For
instance, authors in [48] included both molecular graph and SMILES
representations for retrosynthesis. Although the template-free methods have
gained recent attention, they are known to suffer from the generation of
invalid SMILES. During translation, a single character addition or a missing
one can render the generated reaction invalid.

We note that the current literature although focuses on direct prediction of
products and/or retrosynthesis, there are very limited efforts for CRM
predictions [36]. In [16], authors utilized deep learning on a rather limited,
private dataset of elementary reactions to identify the probable electron
sources and sinks and subsequent ranking of these combinations. More recently,
Authors in [8] proposed ELECTRO, a graph-based generative model for electron
paths movement mimicking ‘arrow-pushing’ diagrams, but did not consider
elementary reaction steps due to the absence of such details in the USPTO-50k
dataset. Consequently, ELECTRO is not suitable to multi-step reactions.

##  3 Preliminaries and Proposed Method

###  3.1 Background

Chemical reaction. A chemical reaction is a process in which one set of
substances (reactants) transforms into another set of substances (products).

Elementary step. The elementary steps in a chemical reaction typically
consists of bond breaking or bond forming involving the
reactants/intermediates, where an electron rich source gets attached to an
electron poor sink. The movement of electrons is usually denoted using ’arrow-
pushing diagrams’ where the arrow direction is from the electron source to the
sink.

Chemical reaction mechanism (CRM). A sequence of elementary steps that
describe the transformation of reactant $\mathbf{R}$ to product $\mathbf{P}$,
through several transient intermediates $\mathbf{I}$, constitutes a full CRM
as shown in eqn(1).

$\displaystyle\mathbf{R}\xrightarrow{stepS_{1}}I_{1}\xrightarrow{stepS_{2}}I_{2},{\dotsm},{I_{n}}\xrightarrow{stepS_{n}}\mathbf{P}$
(1)

where $S_{1:n}$ represents individual elementary steps.

Reactive atom. Reactive atoms are the active atoms involved in an elementary
step that undergoes a large change in their immediate bonding/valency
environment as a result of the reaction.

Reaction templates. Reaction templates are defined as predetermined sets of
chemical transformation rules with specific constraints, such as the presence
of a particular substructure. [26] In Figure 1a, an illustration of a template
is provided using oxidative addition as a representative elementary step in a
reaction. The specified constraints for this template are that any eligible
reactant should possess a substructure featuring an aryl C-X bond (where X =
Cl, Br, I) and that a catalytically active palladium (Pd) metal center be
present.

###  3.2 CRM prediction via ReactAIvate

We develop an interpretable GNN model, ReactAIvate, to predict the elementary
steps, which is further used to devise CRM. An extensive dataset comprising of
seven such mechanistic steps for transition metal-catalyzed reactions is
curated (see Section 4.1 for further details). ReactAIvate is build upon two
elementary tasks:

Reaction step classification (RSC). The primary task of ReactAIvate is to
predict the correct elementary step among the seven identified elementary
steps. These mechanistic steps include ‘oxidative addition’, ‘boron
transmetallation’, ‘acid-base deprotonation’, ‘boronate formation’, ‘substrate
coordination’, ‘transmetallation’, and ‘reductive elimination’.[25] Once the
correct mechanistic step is classified, an off-the-shelf template based
reaction rules are used to predict the product information accurately. The
predicted products form the reactants for the next step, and the subsequent
mechanistic step is identified again. The process is repeated until the
catalyst is regenerated. Our advantage in CRM prediction stems from our focus
on mechanistic step classification, providing a distinct edge over traditional
sequence-based modeling approaches. More specifically, given the set of
reactants $\mathbf{R},$ or intermediates $I_{1}$, $I_{2}$, ..,$I_{n}$,
ReactAIvate predicts labels denoted as $\\{S_{i}\\}_{i=1}^{7}$, where each
$S_{i}$ corresponds to one of the seven mechanistic steps (Figure 1c).

However, while the dataset consists of the seven elementary steps, it doesn’t
encompass the entire range of mechanistic rules that a set of reactants may
undergo in a chemical reaction. To ensure ReactAIvate doesn’t erroneously
predict reactions following different chemical transformation rules or force
chemically non-reactive combinations into the seven elementary steps, we
introduce an eighth class, denoted as $S_{8}$. Our framework is trained to
classify any out-of-distribution (OOD) samples into this eighth class,
bolstering confidence in ReactAIvate’s predictions.

Reactive atom identification (RAI). Reactive step classification alone falls
short in revealing the fundamental mechanisms and rationales underlying a CRM.
Accurate identification of reaction centers within the reactant molecules is
important for forward synthesis. This approach is particularly crucial when
dealing with OOD samples, where the identification of reactive atoms within a
given reaction class holds significant value. These insights serve as a guide
for experts in recognizing feasible elementary reactions. In this work, the
labels for reactive atom classification are derived from reaction templates,
with ‘1’ indicating a reactive atom and ‘0’ indicating a non-reactive atom.

####  3.2.1 ReactAIvate Workflow

We now proceed to outline the operational methodology of our framework,
ReactAIvate.

1\. Molecular representaion via graphs A molecule can be portrayed as a graph,
where atoms and bonds constitute nodes and edges, respectively [31]. Each
unique atom is represented by the following set of features, encompassing nine
types such as atom symbol, formal charge, hybridization, aromaticity, etc.
These collectively result in a total of 39 atom features (see supporting
information for further details). Similarly, edge representation includes
feature vectors corresponding to different bond types (single, double, triple,
etc.).

In our graph representation, $G=(V,E)$, where $V$ is the set of atoms and $E$
is the set of edges, each atom $\nu$ is associated with a feature vector
$X_{\nu}$ in $\mathbb{R}^{D}$, with $D=39$ representing the number of features
for each atom. The tasks at hand involve: (1) classifying elementary steps
into predefined classes. Given an input graph $G$ consisting of reactant
molecules, the goal is to learn the representation vector $h^{*}_{G}$ and a
linear function $g_{1}$ such that the predicted step $\hat{y}_{G}\coloneqq
g_{1}(h^{*}_{G})$ aligns with the true step $y_{G}\in\\{S_{i}\\}{i=1}^{8}$;
(2) identifying reactive atoms, where each atom $\nu$ in $V$ has a label
$y_{\nu}$ ($y_{\nu}\in[0,1]$). The objective is to learn the representation
vector $h^{*}_{\nu}$ and a linear function $g_{2}$ for all $\nu$ such that the
predicted reactivity $\hat{y}_{\nu}\coloneqq g_{2}(h^{*}_{\nu})$ aligns with
the true binary label $y_{\nu}$.

2\. Graph attention network for RSC & RAI After encoding reactant molecules as
graphs, we introduce an attention mechanism to these graphs by generating a
context vector for a target atom $(\nu)$ through attention on its neighboring
atoms ${u}$, where $u\in\mathcal{N}(\nu)$. The initial step in computing the
context vector involves determining attention weights $\alpha_{u\nu}$ between
the state vectors $h_{\nu}$ and $h_{u}$ of the two atoms [47, 52].
Subsequently, a context operation follows, where a linear transformation is
applied to $h_{u}$, the state vectors of neighboring atoms. This is succeeded
by a weighted sum and a non-linear activation function, resulting in
$C_{\nu}$, the context vector for the target atom $\nu$. The calculation of
normalized attention coefficients can be formulated as follows:

$\alpha_{u\nu}=\frac{\exp\left(\texttt{LeakyReLU}\left(\textbf{W}[{h}_{\nu},{h}_{u}]\right)\right)}{\sum_{u\in\mathcal{N}(\nu)}\exp\left(\texttt{LeakyReLU}\left(\textbf{W}[{h}_{\nu},{h}_{u}]\right)\right)},$

where $\alpha_{u\nu}$ signifies the importance (weight) of neighbor atom $u$
to target atom $\nu$, and W is a trainable weight matrix, with

$C_{\nu}=\texttt{ELU}\left({\sum_{u\in\mathcal{N}(\nu)}\alpha_{u\nu}.\textbf{W}.h_{u}}\right).$

Modern GNNs adopt neighborhood aggregation strategies, updating the features
of the target atom iteratively by incorporating the features of its
neighboring atoms. The atom’s representation encapsulates the structural
information within the $k$-hop network around it after $k$ iterations of this
aggregation process. This strategy can be formulated as:

Aggregation phase

$C^{k-1}_{\nu}={\sum_{u\in\mathcal{N}(\nu)}A^{k-1}(h^{k-1}_{u},h^{k-1}_{\nu})},$

Update phase

$h^{k}_{\nu}={\texttt{GRU}^{k-1}(C^{k-1}_{\nu},h^{k-1}_{\nu})},$

where, $h^{k}_{\nu}$ represents the feature vector of atom $\nu$ after the
$k$th layer, with $h^{0}{\nu}=X_{\nu}$ (see Figure 1c). In the aggregation
phase, the graph attention mechanism, $A^{k-1}$, provides the most relevant
information to the target atom from its neighborhoods in the form of the
context vector $C^{k-1}_{\nu}$. In the subsequent step, the update function
$\texttt{GRU}^{k-1}$ (gated recurrent unit) takes the attention context and
the previous state vector of the target atom as input to update the feature
vector from the previous state $h^{k-1}_{\nu}$ to the current state
$h^{k}_{\nu}$.

To generate a graph-level embedding, a virtual _supernode_ is introduced,
connecting with all the atoms in the molecular graph. The feature vector for
this supernode is obtained using sum pooling, expressed as
$h^{0}_{s}=\sum_{\nu}h^{k}_{\nu}$. Through a similar neighborhood aggregation
method, a high-level graph embedding $h^{t}_{s}$ is learned iteratively over
$t$ iterations. At this point, the optimally learned vector $h^{t}_{s}$ can be
considered equivalent to $h^{*}_{G}$ (see Figure 1c), which is subsequently
utilized for the elementary step prediction task. By feeding the graph-level
embedding $h^{*}_{G}$ into a feed-forward neural network (FNN), predicted
values are obtained for the RSC task as:

$\hat{y}_{G}=\texttt{FNN}(h^{*}_{G}).$

ReactAIvate employs cross-entropy loss between the predicted and true labels
for the RSC task:

$\mathcal{L}^{G}_{class}\coloneqq\texttt{CrossEntropyLoss}(\hat{y}_{G},{y}_{G}).$

In our RAI task, the objective is to classify each atom in a molecule as
reactive or non-reactive, enabling the use of standard classification loss
functions. ReactAIvate employs binary cross-entropy (BCE) loss for two-way
classification, with a slight modification. Given that the majority of atoms
in a sample are non-reactive, we modify the BCE loss by introducing a weight
that strongly penalizes incorrect predictions for atoms that are genuinely
reactive. This adjustment helps prevent the model from becoming biased towards
predicting all atoms as non-reactive (see supporting information).

Recall that for an atom $\nu$, the associated embedding vector before the sum
pooling is denoted as $h^{k}_{\nu}$. This embedding vector is considered the
optimal atom representation $h^{*}_{\nu}$. In the node-level classification
task, the updated node feature vector is passed through the FNN to predict
whether the atom is reactive or non-reactive, expressed as:

$\hat{y}_{\nu}=\texttt{FNN}(h^{*}_{\nu}).$

For RAI, we consider the weighted BCE loss as discussed above:

$\mathcal{L}^{G}_{RC}\coloneqq\sum_{\nu}\texttt{WeightedBCELoss}(\hat{y}_{\nu},{y}_{\nu})$

The overall loss for ReactAIvate is made up of the individual losses for RSC
and RAI:

$\mathcal{L}=\mathcal{L}^{G}_{class}+\mathcal{L}^{G}_{RC}$

##  4 Experiments

###  4.1 Dataset details

In this study, we consider a diverse and representative set of transition
metal-catalyzed reactions. These include Suzuki-Miyaura coupling (SMC) [29,
4], Buchwald-Hartwig amination (BHA) [34, 14], and Kumada coupling (KC) [1].
The inclusion of these reactions is motivated by their significance in drugs,
agrochemicals, and pharmaceutical synthesis. To the best of our knowledge,
there is currently no existing database that includes credible mechanisms for
these reactions, compelling the creation of CRM datasets. For this, seven
distinct elementary mechanistic steps are recognized that can account for all
three reaction mechanisms: ‘oxidative addition’, ‘boron transmetallation’,
‘acid-base deprotonation’, ‘boronate formation’, ‘substrate coordination’,
‘transmetallation’, and ‘reductive elimination’ (see supporting information).

We have created essential reaction templates for each elementary mechanistic
step within the three considered reaction datasets in this study. The template
is structured to yield the product based on the given reactants (refer to
Figure 1a). A reaction comprises a set of substrates, catalyst and reagents.
For example, in the case of SMC, key substrates include an aryl halide and a
boronic acid, with a metal-phosphine complex serving as the catalyst and a
base facilitating the reaction. These reaction partners are curated from
primary literature [51, 21]and the PubChem database [23]. By inputting this
set of molecules into the reaction template aligned with the specific reaction
mechanism type, all elementary steps and the complete CRM can be derived. Our
dataset encompasses a total of 100,000 elementary mechanistic steps (see
supporting information for more details).

###  4.2 Training details

The dataset is partitioned into training, validation, and test samples with a
distribution ratio of 70:10:20. ReactAIvate is constructed using PyTorch [32]
with the Adam optimizer [24] and a batch size of 256. Throughout this study,
we maintain consistent hyperparameter values: $k$ (number of attentive message
passing layers for atom embedding) is set to 2, $t$ (number of attentive
message passing layers for molecule embedding) is set to 1, with an $L2$
weight decay of 0.000001, a learning rate of 0.001, and a dropout rate of 0.1.
The number of atom features and graph feature size are specified as 39 and
200, respectively.

###  4.3 Baseline

First, we employed Seq2Seq approaches, such as T5Chem [28] and Transformer
[35] as proxies for baselines, which are considered state-of-the-art (SOTA)
for individual single step reaction prediction. T5Chem is a multi-tasking
model designed for various reaction prediction tasks, leveraging the Text-To-
Text Transfer Transformer (T5), an encoder-decoder model from the transformer
family. The Transformer baseline is built on the original encoder-decoder
framework [46].

In these frameworks, the problem is formulated as a generation task, where the
model is trained to predict the full CRM given the reactants and other
entities. To ensure a fair comparison, we introduce two test datasets, each
comprising 1000 samples. The first is an in-distribution (ID) test dataset,
where the individual reacting partners of each sample belong to the same set
as the training datasets, although the samples themselves are not part of the
training set. The second is an out-of-distribution (OOD) dataset, where the
reaction components of each sample differ structurally from those used in
training the model.

Model | Accuracy (%) | Train
---|---|---
ID test set | OOD test set | time (m)
Seq2Seq |  |  |
T5Chem |  95.60 $\pm$ 0.01 |  0.07 $\pm$ 0.01 | 180
T5Chem(FI) |  98.40 $\pm$ 0.01 |  0.11 $\pm$ 0.02 | 45
Transformer |  11.10 $\pm$ 0.03 |  0.01 $\pm$ 0.00 | 150
Transformer(FI) |  91.80 $\pm$ 0.05 |  0.07 $\pm$ 0.00 | 60
Featurization |  |  |
Morgan(r=2) |  100.00 $\pm$ 0.00 |  31.50 $\pm$ 3.25 | 20
MFF |  100.00 $\pm$ 0.00 |  48.76 $\pm$ 3.60 | 30
PCD |  100.00 $\pm$ 0.00 |  48.57 $\pm$ 2.67 | 10
Graph |  |  |
MPNN |  100.00 $\pm$ 0.00 |  52.80 $\pm$ 9.50 | 10
ReactAIvate | $\mathbf{100.00\pm 0.00}$ | $\mathbf{95.70\pm 0.34}$ | $\mathbf{15}$
Table 1: Performance of different models for CRM prediction on ID and OOD test
molecules. Training time is reported in minutes (m). Bold indicates best
performance.

Next, we developed several DNN models as baselines that are not based on
sequence data such as in seq2seq models. These models use different feature
representations, including fingerprint (Morgan and multiple fingerprint
feature (MFF)) [37], and rdkit-based physicochemical descriptors (PCD)) [33].
Their design targets the same objective as ReactAIvate. Additionally, we
considered a Message-Passing Neural Network (MPNN) model as another benchmark
[18]. It operates by passing messages between nodes (atoms) in a graph
(molecule) to learn the features and interactions. More details about baseline
models are provided in the supporting information.

---
Figure 2: Effect of the inclusion of node-level loss in ReactAIvate
demonstrated through attention visualization. The rightmost bar represents
min-max rescaled attention values.
---
Figure 3: An illustration of the sequential generation of the full CRM for the
Kumada coupling reaction
---
Figure 4: Attention visualization for a sample in (a) non-reactive, (b)
reactive out-of-distribution (OOD) set.

###  4.4 Results

Table 1 illustrates the performance comparison, specifically focusing on
accuracy, i.e., the percentage of samples where the generated CRM precisely
matches the true CRM. For the ID dataset, ReactAIvate exhibits better accuracy
compared to T5Chem and Transformer. Intriguingly, in the OOD dataset,
ReactAIvate offers superior accuracy in CRM generation, reaching upto 96%. In
contrast, both the baselines fail to predict the correct CRM, potentially due
to incorrect predictions of one of the many tokens or atom imbalance,
rendering the entire sequence invalid (see supporting information for more
details). The observation is not specific to CRM prediction, but has also been
reported in recent findings, where several SOTA template-free models
unexpectedly falter when faced with OOD in retrosynthesis prediction [11]. In
addition, unlike the above-mentioned baselines, where T5Chem and Transformer-
based architectures are trained to predict the entire CRM, we have conducted
additional experiments where the baselines are trained to predict just the
forward intermediates (FIs) at each step (exactly what these models were
originally designed for), and then the predicted intermediates serve as the
reactants for the next step. The process is repeated until the catalyst is
regenerated (i.e., end of CRM). The results are reported above in Table 1.
Even for the FI variants of the baselines, their performances remain highly
inadequate for OOD samples. Although T5Chem benefits from pre-training on the
Pubchem database containing 97M molecules, its sequence-centric framework
limits its effectiveness in CRM prediction, especially with OOD instances.

The poor OOD performance of seq2seq/transformer models is not due to poor
hyperparameter tuning but rather the nature of character generation-based
product prediction. For example, if a model is trained on molecules with
Iodine and Chlorine but tested on a molecule where Iodine is replaced with
Fluorine, it struggles because it was never trained to generate sequences
involving Fluorine. Additionally, any incorrect symbol prediction in an
intermediate step can corrupt the entire downstream task for seq2seq models.
On the contrary, ReactAIvate predicts the most appropriate reaction rule
(template) at each step, using these rules to accurately generate the right
products at each step.

The fingerprints and PCD based DNN models exhibited enhanced performance as
compared to sequence-based models on both ID and OOD samples. However, it is
worth noting that these models yield lower accuracy for the OOD set as
compared to ReactAIvate, this discrepancy may stem from the fact that
fingerprints rely on predefined sets of molecular substructures, potentially
limiting their ability to capture nuanced structural details. Whereas MPNN
performed better than fingerprints and PCD featurization, they still fell
short of ReactAIvate. This could be attributed to the advantages conferred by
attention mechanisms and the incorporation of RAI task within ReactAIvate.
This, in turn, implicitly captures the relevance between CRM and RAI, and also
helps with the alignment of highly attentive atoms with likely reaction
centers. Nonetheless, these models provide robust and adequate validation.
These results strongly imply the necessity for a meticulous reconsideration of
the modeling aspect for CRM generation.

Task |  Accuracy (%) | F1-Score
---|---|---
ID-RSC |  100.00 $\pm$ 0.00 |  1.00 $\pm$ 0.00
ID-RAI |  96.40 $\pm$ 0.30 |  0.87 $\pm$ 0.01
OOD-RSC |  98.60 $\pm$ 0.23 |  0.98 $\pm$ 0.01
OOD-RAI |  94.86 $\pm$ 0.40 |  0.85 $\pm$ 0.02
Table 2: Performance of ReactAIvate for RSC and RAI tasks on ID and OOD test
molecules.

We delve deeper into ReactAIvate’s capability for accurately predicting RSC
and RAI concurrently. The classification accuracies and F1-scores for ID and
OOD reactants are presented in Table 2. Beyond accurately classifying the
correct elementary step, ReactAIvate demonstrates proficiency in
distinguishing between reactive and non-reactive atoms. Notably, even in the
OOD dataset, ReactAIvate maintains high efficacy for both RSC and RAI tasks.
This indicates that the model has a broad applicability domain and holds
potential utility for domain experts.

##  5 Discussion and Analysis

Having demonstrated the exceptional performance of our framework ReactAIvate
for both the RSC and RAI tasks, we proceed to delve into the underlying
reasons for its effectiveness through a comprehensive ablation study.
Additionally, we aim to elucidate domain-specific advantages that contribute
to its superior performance.

Significance of incorporating node-level RAI: We begin by emphasizing the
implication of integrating the node-level loss in the reactive atom prediction
task along with the graph-level loss for elementary step prediction. In Figure
2, we compare the attention visualization between ReactAIvate without the
node-level loss and ReactAIvate with both losses, using an oxidative addition
step from the Kumada coupling dataset. In the former model, attention weights
are dispersed across the reactants and catalyst, providing limited utility.
Notably, the crucial reactive metal atom Pd:7, responsible for catalyzing this
step, fails to draw any attention. Thus, attention visualization without the
node-level loss contributes little to understanding chemical reactivity. In
the latter model with the two-level loss, attention weights are more
concentrated around the reactive region of the molecules. For instance, the Pd
metal center and the aryl chlorine bond garner higher attention, (Pd:7, C:20,
Cl:21), representing the true reactive atoms in this elementary step. This
underscores the importance of including the node-level loss, aligning the
model’s attention mechanism more closely with how a chemist would assign
attention to reactive atoms.

Illustration of a full CRM: Moving forward, we assess ReactAIvate’s capability
in predicting a complete CRM for Kumada coupling, serving as a representative
example (Figure 3). In the initial step, the model accurately predicts
oxidative addition ($S_{1}$) as the elementary step, along with the correct
identification of reactive atoms. Notably, the model also predicts two
additional atoms as active, which intriguingly turn out to be reactive in the
subsequent step. In the second step, ReactAIvate once again correctly predicts
the elementary class, transmetallation ($S_{6}$), and accurately identifies
all true reactive atoms. Finally, the model predicts the reductive elimination
($S_{7}$) step as the concluding phase. Both the predicted active atoms and
attention distributions align consistently with the expected mechanism. In
summary, ReactAIvate demonstrates the ability to generate a complete CRM
starting from only reactants and catalysts.

Attention visualization for OOD-class: To gain deeper insights into the model,
we visualize attentions in a sample belonging to the OOD class, as shown in
Figure 4a. In this instance, the pair of molecules is predicted to fall into
the eighth (OOD) class, accurately reflecting that the amine and the aryl
halide (N:27, C:18, Cl:19) would not react in the absence of a catalyst.
Intriguingly, the model predicts chlorine in one molecule and the adjacent
carbon in the other as reactive atoms. Additionally, attentions are primarily
distributed around functional groups. The broad dispersion of attention
underscores ReactAIvate’s challenge in pinpointing the exact reaction
mechanism. This outcome aligns with expectations, as the combination is
chemically non-reactive, and it would be counterintuitive for the model to
highlight specific reaction centers leading to its prediction in the OOD
class.

In our final evaluation, we test ReactAIvate for the identification of
potential reactive atoms in an OOD sample involving an entirely different
reaction mechanism, as depicted in Figure 4b. Please note that the molecule is
not expected to undergo any transformation (both predicted and true class are
$S_{8}$, i.e., no reaction). However, if it is presented with a suitable
reagent, this molecule is anticipated to undergo a $\beta$ hydride elimination
step and the reactive centers get activated. Among the predicted reactive
atoms, two (Pd:13 and $\alpha$-C: 6) correspond to possible true reactive
atoms. Since the molecule in this example doesn’t undergo reaction, the true
reactive centers are presented as an empty list. Moreover, attentions are
dispersed around the reaction center and the reactive atoms. These highlights
the model’s versatility in predicting OOD samples and serves as a guide for
comprehending the reactivity of metal-catalyzed reactions. Consequently, one
can easily incorporate new elementary classes of interest to further broaden
the model’s applicability.

##  6 Conclusion and Future Work

In conclusion, we introduce ReactAIvate, a graph-attention-based Graph Neural
Network (GNN) model designed for interpretable Chemical Reaction Mechanism
(CRM) generation. Our model is trained on a novel dataset comprising seven
distinct elementary mechanistic steps, covering the complete CRM for three
different transition-metal-catalyzed processes. ReactAIvate excels in
accurately classifying elementary steps and recognizing reactive atoms,
demonstrating its capability to construct full CRMs. Notably, the model
exhibits a prudent handling of non-reactive cases, showcasing its reliability
in predictions. ReactAIvate outperforms Seq2Seq baseline models, emphasizing
the limitations of the latter in CRM identification due to minor errors. The
robust OOD classification performance underscores ReactAIvate’s potential for
exploring additional mechanisms with the availability of more data. As part of
future work, we plan create a user-friendly interface for predicting entire
CRMs based on user-inputted SMILES of reactants.

##  7 Data Availability

Data and codes related to this work are publicly available through our Github
repository at https://github.com/alhqlearn/ReactAIvate.

##  8 Acknowledgements

We gratefully acknowledge the generous computing time provided by the
SpaceTime supercomputing facility at IIT Bombay. M.D. expresses gratitude for
the Prime Minister’s Research Fellowship.

## References

  * Ackermann and Althammer [2006] L. Ackermann and A. Althammer.  Air-stable pinp (o) h as preligand for palladium-catalyzed kumada couplings of unactivated tosylates.  _Organic Letters_ , 8(16):3457–3460, 2006.
  * Ahneman et al. [2018] D. T. Ahneman, J. G. Estrada, S. Lin, S. D. Dreher, and A. G. Doyle.  Predicting reaction performance in c–n cross-coupling using machine learning.  _Science_ , 360(6385):186–190, 2018.
  * Bahmanyar et al. [2003] S. Bahmanyar, K. Houk, H. J. Martin, and B. List.  Quantum mechanical predictions of the stereoselectivities of proline-catalyzed asymmetric intermolecular aldol reactions.  _Journal of the American Chemical Society_ , 125(9):2475–2479, 2003.
  * Beletskaya et al. [2019] I. P. Beletskaya, F. Alonso, and V. Tyurin.  The suzuki-miyaura reaction after the nobel prize.  _Coordination Chemistry Reviews_ , 385:137–173, 2019.
  * Benzinger [1971] T. Benzinger.  Thermodynamics, chemical reactions and molecular biology.  _Nature_ , 229(5280):100–102, 1971.
  * Bhattacharya et al. [2023] T. Bhattacharya, S. Ghosh, S. Dutta, S. Guin, A. Ghosh, H. Ge, R. B. Sunoj, and D. Maiti.  Combinatorial ligand assisted simultaneous control of axial and central chirality in highly stereoselective c- h allylation.  _Angewandte Chemie_ , page e202310112, 2023.
  * Born and Manica [2021] J. Born and M. Manica.  Trends in deep learning for property-driven drug design.  _Current medicinal chemistry_ , 28(38):7862–7886, 2021.
  * Bradshaw et al. [2018] J. Bradshaw, M. J. Kusner, B. Paige, M. H. Segler, and J. M. Hernández-Lobato.  A generative model for electron paths.  _arXiv preprint arXiv:1805.10970_ , 2018.
  * Brown and Bostrom [2016] D. G. Brown and J. Bostrom.  Analysis of past and present synthetic methodologies on medicinal chemistry: where have all the new reactions gone? miniperspective.  _Journal of medicinal chemistry_ , 59(10):4443–4458, 2016.
  * Chen and Jung [2022] S. Chen and Y. Jung.  A generalized-template-based graph neural network for accurate organic reactivity prediction.  _Nature Machine Intelligence_ , 4(9):772–780, 2022.
  * Chen and Jung [2024] S. Chen and Y. Jung.  Assessing the extrapolation capability of template-free retrosynthesis models.  _arXiv preprint arXiv:2403.03960_ , 2024.
  * Chin et al. [2022] Y. P. Chin, N. W. See, I. D. Jenkins, and E. H. Krenske.  Computational discoveries of reaction mechanisms: recent highlights and emerging challenges.  _Organic & Biomolecular Chemistry_, 20(10):2028–2042, 2022.
  * Coley et al. [2017] C. W. Coley, R. Barzilay, T. S. Jaakkola, W. H. Green, and K. F. Jensen.  Prediction of organic reaction outcomes using machine learning.  _ACS central science_ , 3(5):434–443, 2017.
  * Dorel et al. [2019] R. Dorel, C. P. Grugel, and A. M. Haydl.  The buchwald–hartwig amination after 25 years.  _Angewandte Chemie International Edition_ , 58(48):17118–17129, 2019.
  * Flam-Shepherd et al. [2022] D. Flam-Shepherd, K. Zhu, and A. Aspuru-Guzik.  Language models can learn complex molecular distributions.  _Nature Communications_ , 13(1):3293, 2022.
  * Fooshee et al. [2018] D. Fooshee, A. Mood, E. Gutman, M. Tavakoli, G. Urban, F. Liu, N. Huynh, D. Van Vranken, and P. Baldi.  Deep learning for chemical reaction prediction.  _Molecular Systems Design & Engineering_, 3(3):442–452, 2018.
  * Genheden et al. [2020] S. Genheden, A. Thakkar, V. Chadimová, J.-L. Reymond, O. Engkvist, and E. Bjerrum.  Aizynthfinder: a fast, robust and flexible open-source software for retrosynthetic planning.  _Journal of cheminformatics_ , 12(1):70, 2020\.
  * Gilmer et al. [2017] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl.  Neural message passing for quantum chemistry.  In _International conference on machine learning_ , pages 1263–1272. PMLR, 2017.
  * Goldstein et al. [1996] E. Goldstein, B. Beno, and K. Houk.  Density functional theory prediction of the relative energies and isotope effects for the concerted and stepwise mechanisms of the diels- alder reaction of butadiene and ethylene.  _Journal of the American Chemical Society_ , 118(25):6036–6043, 1996.
  * Heid et al. [2021] E. Heid, J. Liu, A. Aude, and W. H. Green.  Influence of template size, canonicalization, and exclusivity for retrosynthesis and reaction prediction applications.  _Journal of Chemical Information and Modeling_ , 62(1):16–26, 2021.
  * Huang et al. [2003] X. Huang, K. W. Anderson, D. Zim, L. Jiang, A. Klapars, and S. L. Buchwald.  Expanding pd-catalyzed c- n bond-forming processes: the first amidation of aryl sulfonates, aqueous amination, and complementarity with cu-catalyzed reactions.  _Journal of the American Chemical Society_ , 125(22):6653–6655, 2003.
  * Jin et al. [2017] W. Jin, C. Coley, R. Barzilay, and T. Jaakkola.  Predicting organic reaction outcomes with weisfeiler-lehman network.  _Advances in neural information processing systems_ , 30, 2017.
  * Kim et al. [2021] E. Kim, D. Lee, Y. Kwon, M. S. Park, and Y.-S. Choi.  Valid, plausible, and diverse retrosynthesis using tied two-way transformers with latent variables.  _Journal of Chemical Information and Modeling_ , 61(1):123–133, 2021.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kurosawa and Yamamoto [2003] H. Kurosawa and A. Yamamoto.  _Fundamentals of molecular catalysis_.  Elsevier, 2003.
  * Landrum [2013] G. Landrum.  Rdkit documentation.  _Release_ , 1(1-79):4, 2013.
  * Lin et al. [2020] K. Lin, Y. Xu, J. Pei, and L. Lai.  Automatic retrosynthetic route planning using template-free models.  _Chemical science_ , 11(12):3355–3364, 2020.
  * Lu and Zhang [2022] J. Lu and Y. Zhang.  Unified deep learning model for multitask reaction predictions with explanation.  _Journal of Chemical Information and Modeling_ , 62(6):1376–1387, 2022.
  * Martin and Buchwald [2008] R. Martin and S. L. Buchwald.  Palladium-catalyzed suzuki- miyaura cross-coupling reactions employing dialkylbiaryl phosphine ligands.  _Accounts of chemical research_ , 41(11):1461–1473, 2008.
  * Meng et al. [2023a] Z. Meng, P. Zhao, Y. Yu, and I. King.  Doubly stochastic graph-based non-autoregressive reaction prediction.  _In IJCAI_ , 2023a.
  * Meng et al. [2023b] Z. Meng, P. Zhao, Y. Yu, and I. King.  A unified view of deep learning for reaction and retrosynthesis prediction: Current status and future challenges.  _In IJCAI_ , 2023b.
  * Paszke et al. [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.
  * Riniker and Landrum [2013] S. Riniker and G. A. Landrum.  Open-source platform to benchmark fingerprints for ligand-based virtual screening.  _Journal of cheminformatics_ , 5(1):26, 2013.
  * Ruiz-Castillo and Buchwald [2016] P. Ruiz-Castillo and S. L. Buchwald.  Applications of palladium-catalyzed c–n cross-coupling reactions.  _Chemical reviews_ , 116(19):12564–12649, 2016\.
  * Rush [2018] A. M. Rush.  The annotated transformer.  In _Proceedings of workshop for NLP open source software (NLP-OSS)_ , pages 52–60, 2018.
  * Sakai et al. [2023] M. Sakai, M. Kaneshige, and K. Yasuda.  Learning organo-transition metal catalyzed reactions by graph neural networks.  _Journal of Computational Chemistry_ , 2023.
  * Sandfort et al. [2020] F. Sandfort, F. Strieth-Kalthoff, M. Kühnemund, C. Beecks, and F. Glorius.  A structure-based platform for predicting chemical reactivity.  _Chem_ , 6(6):1379–1390, 2020.
  * Schneider et al. [2016] N. Schneider, N. Stiefl, and G. A. Landrum.  What’s what: The (nearly) definitive guide to reaction role assignment.  _Journal of chemical information and modeling_ , 56(12):2336–2346, 2016.
  * Schupp et al. [2024] M. Schupp, I. Saridakis, D. Kaiser, and N. Maulide.  Chemical synthesis as a discovery platform in immunosuppression and determination of mode of action.  _Nature Synthesis_ , pages 1–13, 2024.
  * Schwaller et al. [2019] P. Schwaller, T. Laino, T. Gaudin, P. Bolgar, C. A. Hunter, C. Bekas, and A. A. Lee.  Molecular transformer: a model for uncertainty-calibrated chemical reaction prediction.  _ACS central science_ , 5(9):1572–1583, 2019\.
  * Segler and Waller [2017a] M. H. Segler and M. P. Waller.  Modelling chemical reasoning to predict and invent reactions.  _Chemistry–A European Journal_ , 23(25):6118–6128, 2017a.
  * Segler and Waller [2017b] M. H. Segler and M. P. Waller.  Neural-symbolic machine learning for retrosynthesis and reaction prediction.  _Chemistry–A European Journal_ , 23(25):5966–5971, 2017b.
  * Segler et al. [2018] M. H. Segler, M. Preuss, and M. P. Waller.  Planning chemical syntheses with deep neural networks and symbolic ai.  _Nature_ , 555(7698):604–610, 2018.
  * Tavakoli et al. [2023] M. Tavakoli, Y. T. T. Chiu, A. Shmakov, A. M. Carlton, D. Van Vranken, and P. Baldi.  Ai for interpretable chemistry: Predicting radical mechanistic pathways via contrastive learning.  _arXiv preprint arXiv:2311.01118_ , 2023.
  * Ucak et al. [2022] U. V. Ucak, I. Ashyrmamatov, J. Ko, and J. Lee.  Retrosynthetic reaction pathway prediction through neural machine translation of atomic environments.  _Nature communications_ , 13(1):1186, 2022.
  * Vaswani et al. [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Veličković et al. [2017] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Lio, and Y. Bengio.  Graph attention networks.  _arXiv preprint arXiv:1710.10903_ , 2017.
  * Wan et al. [2022] Y. Wan, B. Liao, C.-Y. Hsieh, and S. Zhang.  Retroformer: Pushing the limits of interpretable end-to-end retrosynthesis transformer.  _arXiv preprint arXiv:2201.12475_ , 2022.
  * Wang et al. [2022] Y. Wang, Y. Gu, C. Lou, Y. Gong, Z. Wu, W. Li, Y. Tang, and G. Liu.  A multitask gnn-based interpretable model for discovery of selective jak inhibitors.  _Journal of cheminformatics_ , 14(1):16, 2022\.
  * Wei et al. [2016] J. N. Wei, D. Duvenaud, and A. Aspuru-Guzik.  Neural networks for the prediction of organic chemistry reactions.  _ACS central science_ , 2(10):725–732, 2016.
  * Wolfe and Buchwald [1997] J. P. Wolfe and S. L. Buchwald.  Improved functional group compatibility in the palladium-catalyzed amination of aryl bromides.  _Tetrahedron Letters_ , 38(36):6359–6362, 1997\.
  * Xiong et al. [2019] Z. Xiong, D. Wang, X. Liu, F. Zhong, X. Wan, X. Li, Z. Li, X. Luo, K. Chen, H. Jiang, et al.  Pushing the boundaries of molecular representation for drug discovery with the graph attention mechanism.  _Journal of medicinal chemistry_ , 63(16):8749–8760, 2019.
  * Zheng et al. [2019] S. Zheng, J. Rao, Z. Zhang, J. Xu, and Y. Yang.  Predicting retrosynthetic reactions using self-corrected transformer neural networks.  _Journal of chemical information and modeling_ , 60(1):47–55, 2019.


# Asymptotics of a Mathieu-Gaussian series

R. B. Paris
Division of Computing and Mathematics,
Abertay University, Dundee DD1 1HG, UK

###### Abstract

We consider the asymptotic expansion of the functional series

$S_{\mu,\gamma}(a;\lambda)=\sum_{n=1}^{\infty}\frac{n^{\gamma}e^{-\lambda
n^{2}/a^{2}}}{(n^{2}+a^{2})^{\mu}}$

for real values of the parameters $\gamma$, $\lambda>0$ and $\mu\geq 0$ as
$|a|\to\infty$ in the sector $|\arg\,a|<\pi/4$. For general values of $\gamma$
the expansion is of algebraic type with terms involving the Riemann zeta
function and a terminating confluent hypergeometric function. Of principal
interest in this study is the case corresponding to even integer values of
$\gamma$, where the algebraic-type expansion consists of a finite number of
terms together with a contribution comprising an infinite sequence of
increasingly subdominant exponentially small expansions. This situation is
analogous to the well-known Poisson-Jacobi formula corresponding to the case
$\mu=\gamma=0$. Numerical examples are provided to illustrate the accuracy of
these expansions.

Mathematics subject classification (2010): 30E15, 33E20, 33C15, 34E05, 41A60

Keywords: Mathieu series, asymptotic expansions, exponential asymptotics,
Mellin transform method

1\. Introduction

The functional series

$\sum_{n=1}^{\infty}\frac{n}{(n^{2}+a^{2})^{\mu}}$ (1.1)

in the case $\mu=2$ was introduced by Mathieu in his 1890 book [2] dealing
with the elasticity of solid bodies. The asymptotic expansion for large $a$ of
more general functional series has been discussed in [4] and [11]. More
recently, Gerhold and Tomovski [1] extended the asymptotic study of such
Mathieu series by introducing in (1.1) the factor $z^{n}$, where $|z|\leq 1$.
From this result they were able to deduce, in particular, the large-$a$
expansions of the trigonometric Mathieu series

$\sum_{n=1}^{\infty}\frac{n\sin
nx}{(n^{2}+a^{2})^{\mu}},\qquad\sum_{n=1}^{\infty}\frac{n\cos
nx}{(n^{2}+a^{2})^{\mu}}.$

Subsequently, the above trigonometric series were generalised to include the
oscillatory Bessel functions $J_{\nu}(x)$ and $Y_{\nu}(x)$ with argument
proportional to $n/a$, and their large-$a$ asymptotics determined in [6]. In
addition, this last study also considered the inclusion of the modified Bessel
function $K_{\nu}(x)$ of similar argument, which contains the decaying
exponential as a special case.

The asymptotic expansion we consider in this paper is the Mathieu series
coupled with a Gaussian exponential of the form

$S_{\mu,\gamma}(a;\lambda):=\sum_{n=1}^{\infty}\frac{n^{\gamma}e^{-\lambda
n^{2}/a^{2}}}{(n^{2}+a^{2})^{\mu}}\qquad(\mu\geq 0,\ \lambda>0)$ (1.2)

for $|a|\rightarrow\infty$ in the sector $|\arg\,a|<\pi/4$. It will be
supposed throughout that $\gamma$ is real, although the analysis is easily
modified to incorporate complex $\gamma$. We shall employ the Mellin transform
approach used in [4, 6, 11], where our interest will be primarily concerned
with even integer values of $\gamma$ (positive or negative). We shall find
that the asymptotic expansion of $S_{\mu,\gamma}(a;\lambda)$ with these
parameter values for large complex $a$ in the sector $|\arg\,a|<\pi/4$
consists of a finite algebraic expansion together with an infinite sequence of
increasingly subdominant exponentially small contributions.

It is interesting that the apparently simple series (1.2) should possess such
an intricate asymptotic structure in the case of even integer values of
$\gamma$. This is also found to be the case when $\lambda=0$ in (1.2); see [5]
for details. A well-known related series corresponding to $\mu=\gamma=0$ is
the Poisson-Jacobi formula [10, p. 124]

$S_{0,0}(a;\lambda)=\sum_{n=1}^{\infty}e^{-\lambda
n^{2}/a^{2}}=\frac{a}{2}\sqrt{\frac{\pi}{\lambda}}-\frac{1}{2}+a\sqrt{\frac{\pi}{\lambda}}\sum_{n=1}^{\infty}e^{-\pi^{2}n^{2}a^{2}/\lambda}.$
(1.3)

This sum is also seen to consist of a finite algebraic contribution together
with an infinite sum of exponentially small terms when $|a|\to\infty$ in the
sector $|\arg\,a|<\pi/4$.

In the application of the Mellin transform method to the series in (1.2) and
its alternating variant we shall require the following estimates for the gamma
function and the Riemann zeta function. For real $\sigma$ and $t$, we have the
estimates

$\Gamma(\sigma\pm it)=O(t^{\sigma-\frac{1}{2}}e^{-\frac{1}{2}\pi
t}),\qquad|\zeta(\sigma\pm
it)|=O(t^{\Omega(\sigma)}\log^{\alpha}t)\quad(t\rightarrow+\infty),$ (1.4)

where $\Omega(\sigma)=0$ ($\sigma>1$),
$\mbox{${\textstyle\frac{1}{2}}$}-\mbox{${\textstyle\frac{1}{2}}$}\sigma$
($0\leq\sigma\leq 1$), $\mbox{${\textstyle\frac{1}{2}}$}-\sigma$ ($\sigma<0$)
and $\alpha=1$ ($0\leq\sigma\leq 1$), $\alpha=0$ otherwise [9, p. 95]. The
zeta function $\zeta(s)$ has a simple pole of unit residue at $s=1$ and the
evaluations for positive integer $k$

$\zeta(0)=-\mbox{${\textstyle\frac{1}{2}}$},\quad\zeta(-2k)=0,\quad\zeta(2k)=\frac{(2\pi)^{2k}}{2(2k)!}|B_{2k}|\
\ (k\geq 1),$ (1.5)

where $B_{k}$ are the Bernoulli numbers. Finally, we have the well-known
functional relation satisfied by $\zeta(s)$ given by [3, p. 603]

$\zeta(s)=2^{s}\pi^{s-1}\zeta(1-s)\Gamma(1-s)\sin\mbox{${\textstyle\frac{1}{2}}$}\pi
s.$ (1.6)

2\. An integral representation

The generalised Mathieu series defined in (1.2) can be written as

$S_{\mu,\gamma}(a;\lambda)=a^{-\delta}\sum_{n=1}^{\infty}h(n/a),\qquad
h(x):=\frac{x^{\gamma}e^{-\lambda
x^{2}}}{(1+x^{2})^{\mu}},\quad\delta:=2\mu-\gamma,$ (2.1)

where the parameter $\delta$ is real and $\lambda>0$. We employ a Mellin
transform approach as discussed in [7, Section 4.1.1]. The Mellin transform of
$h(x)$ is ${\cal H}(s)=\int_{0}^{\infty}x^{s-1}h(x)\,dx$, where [3, (13.4.4)]

${\cal H}(s)=\int_{0}^{\infty}\frac{x^{\gamma+s-1}e^{-\lambda
x^{2}}}{(1+x^{2})^{\mu}}\,dx=\mbox{${\textstyle\frac{1}{2}}$}\Gamma(\mbox{${\textstyle\frac{\gamma+s}{2}}$})\,U(\mbox{${\textstyle\frac{\gamma+s}{2}}$},1+\mbox{${\textstyle\frac{\gamma+s}{2}}$}-\mu,\lambda)$
(2.2)

in the half-plane $\Re(s)>-\gamma$, with $U(a,b,z)$ being the confluent
hypergeometric function of the second kind. The transform ${\cal H}(s)$ can be
represented alternatively in the form

${\cal H}(s)=\mbox{${\textstyle\frac{1}{2}}$}\\{{\cal H}_{1}(s)+{\cal
H}_{2}(s)\\},$ (2.3)

where

${\cal
H}_{1}(s)=\frac{\Gamma(\frac{\gamma+s}{2})\Gamma(\mu-\frac{\gamma+s}{2})}{\Gamma(\mu)}\,{}_{1}F_{1}(\mbox{${\textstyle\frac{\gamma+s}{2}}$};1+\mbox{${\textstyle\frac{\gamma+s}{2}}$}-\mu;\lambda),$
(2.4) ${\cal
H}_{2}(s)=\lambda^{\mu-(\gamma+s)/2}\Gamma(\mbox{${\textstyle\frac{\gamma+s}{2}}$}-\mu)\,{}_{1}F_{1}(\mu;1-\mbox{${\textstyle\frac{\gamma+s}{2}}$}+\mu;\lambda).$
(2.5)

Using the Mellin inversion theorem (see, for example, [7, p. 118]), we find

$S_{\mu,\gamma}(a;\lambda)=\frac{a^{-\delta}}{2\pi
i}\sum_{n=1}^{\infty}\int_{c-\infty i}^{c+\infty i}{\cal
H}(s)(n/a)^{-s}ds=\frac{a^{-\delta}}{2\pi i}\int_{c-\infty i}^{c+\infty
i}{\cal H}(s)\zeta(s)a^{s}ds,$ (2.6)

where $\zeta(s)$ is the Riemann zeta function and $c>\max\\{1,-\gamma\\}$. The
inversion of the order of summation and integration is justified by absolute
convergence provided $c$ satisfies this condition. From the estimates in (1.4)
and the fact that the confluent hypergeometric functions111The ${}_{1}F_{1}$
function appearing in ${\cal H}_{1}(s)$ can be written as
$e^{\lambda}\,{}_{1}F_{1}(1-\mu;1+\mbox{${\textstyle\frac{1}{2}}$}(\gamma+s)-\mu;-\lambda)$
by Kummer’s transformation [3, p. 325]. appearing in ${\cal H}_{1}(s)$ and
${\cal H}_{2}(s)$ are both $O(1)$ as $\Im(s)\to\pm\infty$, the integral in
(2.6) then defines $S_{\mu,\gamma}(a;\lambda)$ for complex $a$ in the sector
$|\arg\,a|<\pi/4$. The integration path in (2.6) lies to the right of the
simple pole of $\zeta(s)$ at $s=1$ and the poles of $\Gamma((\gamma+s)/2)$ at
$s=-\gamma-2k$ ($k=0,1,2,\ldots\,$), these being the only poles of the
integrand since the $U$ function in (2.2) has no poles; see Appendix A for a
demonstration of this fact.

We consider the integral in (2.6) taken round the rectangular contour with
vertices at $c\pm iT$ and $-c^{\prime}\pm iT$, where $c^{\prime}>0$ and $T>0$.
The contribution from the upper and lower sides of the rectangle $s=\sigma\pm
iT$, $-c^{\prime}\leq\sigma\leq c$, vanishes as $T\rightarrow\infty$ provided
$|\arg\,a|<\pi/4$, since from (1.4), the modulus of the integrand is
controlled by $O(T^{\Omega(\sigma)+(\sigma-\delta-1)/2}\log\,T\,e^{-\Delta
T})$, where $\Delta=\pi/4-|\arg\,a|$. Displacement of the integration path to
the left over the pole at $s=1$ and those of ${\cal H}_{1}(s)$ at
$s=-2k-\gamma$ (when $\mu>0$) then yields

$S_{\mu,\gamma}(a;\lambda)-a^{1-\delta}{\cal H}(1)\hskip 227.62204pt$ $\hskip
113.81102pt\sim
a^{-2\mu}\sum_{k=0}^{\infty}\frac{(-)^{k}(\mu)_{k}}{a^{2k}}\zeta(-2k-\gamma)\,{}_{1}F_{1}(-k;1-\mu-k;\lambda),$
(2.7)

where

${\cal
H}(1)=\frac{1}{2}\Gamma(\mbox{${\textstyle\frac{1+\gamma}{2}}$})\,U(\mbox{${\textstyle\frac{1+\gamma}{2}}$},1+\mbox{${\textstyle\frac{1+\gamma}{2}}$}-\mu,\lambda).$
(2.8)

When $\mu=0$, we have ${\cal H}_{1}(s)\equiv 0$ and the poles from ${\cal
H}_{2}(s)$ at $s=-2k-\gamma$ yield

$S_{0,\gamma}(a;\lambda)-\frac{1}{2}\biggl{(}\frac{a^{2}}{\lambda}\biggr{)}^{(\gamma+1)/2}\,\Gamma(\frac{\gamma+1}{2})\sim\sum_{k=0}^{\infty}\frac{(-\lambda)^{k}}{k!\,a^{2k}}\,\zeta(-2k-\gamma),$
(2.9)

where we have used the fact that $U(\alpha,1+\alpha,x)=x^{-\alpha}$. Both
expansions (2.7) and (2.9) hold as $|a|\to\infty$ in $|\arg\,a|<\pi/4$. The
hypergeometric functions appearing in the sum in (2.7) are polynomials in
$\lambda$ of degree $k$ when $\mu$ is non-integer; for positive integer $\mu$
they can be expressed by Kummer’s theorem as $e^{\lambda}$ multiplied by a
polynomial in $\lambda$ of degree $\mu-1$.

The expansion (2.7) holds for general values of $\mu>0$ and $\gamma$. If
$\gamma$ equals an odd negative integer, the pole at $s=1$ has both a double
pole contribution (resulting from ${\cal H}_{1}(s)$) and a simple pole
contribution (resulting from ${\cal H}_{2}(s)$). An example of the expansion
when $\gamma=-1$ is discussed in Appendix B. We remark that when $\gamma=0$ in
(2.9) the expansion for $S_{0,0}(a;\lambda)$ correctly reduces to the first
two terms in the Poisson-Jacobi formula (1.3), but does not account for the
exponentially small contribution. Further consideration of this case is
discussed at the end of Section 3.

Finally, we note that when $\gamma=2p$ is an even integer, there will be a
finite number of poles of the integrand of the sequence $s=-2k-2p$ on account
of the trivial zeros of $\zeta(s)$ at $s=-2,-4,\ldots\ $. This results in the
number of terms in the asymptotic series in (2.7) and (2.9) being either
finite or zero. This situation is the main subject of this paper. We shall
show that, in addition to a finite algebraic contribution, there is a sequence
of increasingly subdominant exponentially small terms in the large-$a$ limit.
This is analogous to the exponentially small contribution appearing on the
right-hand side of the Poisson-Jacobi formula (1.3).

3\. The exponentially small contribution to $S_{\mu,\gamma}(a;\lambda)$ when
$\gamma=2p$

Let $\gamma=2p$ be an even integer and $\mu\geq 0$. Then the quantity $\delta$
defined in (2.1) is $\delta=2(\mu-p)$. The number of poles of the sequence
$s=-2k-2p$ is finite (when $p\leq 0$) or zero (when $p\geq 1$) on account of
the trivial zeros of $\zeta(s)$. Then we have upon displacement of the
integration path

$S_{\mu,\gamma}(a;\lambda)=a^{1-\delta}{\cal
H}(1)+H_{\mu,\gamma}(a;\lambda)+J(a;\lambda),$ (3.1)

where

$H_{\mu,\gamma}(a;\lambda)=\left\\{\begin{array}[]{ll}0&(p\geq 1)\\\
-\mbox{${\textstyle\frac{1}{2}}$}a^{-\delta}&(p=0)\\\
a^{-\delta}\sum_{k=0}^{|p|}\zeta(2k)R_{k}(\mu,|p|)a^{2k}&(p\leq-1).\end{array}\right.$
(3.2)

The quantity $R_{k}(\mu,|p|)$ denotes the residue of ${\cal H}(s)$ at $s=2k$,
$0\leq k\leq|p|$ when $p\leq-1$ given by

$R_{k}(\mu,|p|)=\frac{(-)^{q-k}}{(q-k)!}\,U(k-q,1+k-q-\mu,\lambda)\qquad(q=-p,\
0\leq k\leq q).$

Routine calculations show that when $p=-1,-2$, for example, we have

$\displaystyle R_{0}(\mu,1)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!-(\mu+\lambda),\quad R_{1}(\mu,1)=1,$ $\displaystyle
R_{0}(\mu,2)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\mbox{${\textstyle\frac{1}{2}}$}\mu(1+\mu)+\mu\lambda+\mbox{${\textstyle\frac{1}{2}}$}\lambda^{2},\quad
R_{1}(\mu,2)=-(\mu+\lambda),\quad R_{2}(\mu,2)=1.$

The integral $J(a;\lambda)$ appearing in (3.1) is defined by

$J(a;\lambda)=\frac{a^{-\delta}}{2\pi i}\int_{c-\infty i}^{c+\infty i}{\cal
H}(-s)\zeta(-s)a^{-s}ds\qquad(c>0),$

where we have replaced $s$ by $-s$. Use of the functional relation for
$\zeta(s)$ in (1.6) followed by expansion of $\zeta(1+s)$ (permissible since
$c>0$) leads to

$\displaystyle J(a;\lambda)$ $\displaystyle=$
$\displaystyle-\frac{a^{-\delta}}{2\pi i}\int_{c-\infty i}^{c+\infty i}{\cal
H}(-s)\zeta(1+s)\Gamma(1+s)\frac{\sin\mbox{${\textstyle\frac{1}{2}}$}\pi
s}{\pi}\,(2\pi a)^{-s}ds$ $\displaystyle=$
$\displaystyle-\frac{a^{-\delta}}{\pi}\sum_{k\geq
1}\frac{1}{k}\cdot\frac{1}{2\pi i}\int_{c-\infty i}^{c+\infty i}{\cal
H}(-s)\Gamma(1+s)\sin\mbox{${\textstyle\frac{1}{2}}$}\pi s\,(2\pi ka)^{-s}ds.$

From (2.2) and an application of Kummer’s transformation [3, (13.2.40)] we
have

${\cal
H}(-s)=\mbox{${\textstyle\frac{1}{2}}$}\Gamma(p-\mbox{${\textstyle\frac{1}{2}}$}s)\,\lambda^{\mu-p+s/2}U(\mu,1+\mu-p+\mbox{${\textstyle\frac{1}{2}}$}s,\lambda),$

whence

$J(a;\lambda)=\frac{(-)^{p}a^{-\delta}}{2\lambda^{p-\mu}}\sum_{k\geq
1}\frac{1}{k}\cdot\frac{1}{2\pi i}\int_{c-\infty i}^{c+\infty
i}\frac{2^{-s}\Gamma(1+s)}{\Gamma(1-p+\mbox{${\textstyle\frac{1}{2}}$}s)}\,X_{k}^{-s/2}\,U(\mu,1\\!+\\!\mu\\!-\\!p\\!+\\!\mbox{${\textstyle\frac{1}{2}}$}s,\lambda)\,ds,$
(3.3)

where

$X_{k}:=\frac{\pi^{2}k^{2}a^{2}}{\lambda}.$ (3.4)

We first consider the case $p=0,1,2,\ldots\,$. Since there are no poles of the
integrand in (3.3) in $\Re(s)>0$ the integration path can be displaced as far
to the right as we please such that $|s|$ is everywhere large on the new path.
The quotient of gamma functions in the integrand can then be expanded as [7,
p. 53], [3, (5.11.19)]

$\frac{2^{-s}\sqrt{\pi}\Gamma(1+s)}{\Gamma(1\\!-\\!p\\!+\\!\mbox{${\textstyle\frac{1}{2}}$}s)}=\frac{\Gamma(\mbox{${\textstyle\frac{1}{2}}$}+\mbox{${\textstyle\frac{1}{2}}$}s)\Gamma(1+\mbox{${\textstyle\frac{1}{2}}$}s)}{\Gamma(1\\!-\\!p\\!+\\!\mbox{${\textstyle\frac{1}{2}}$}s)}=\sum_{j=0}^{p}(-)^{j}c_{j}\Gamma(\mbox{${\textstyle\frac{1}{2}}$}s+\vartheta-j),$
(3.5)

where $\vartheta=p+\mbox{${\textstyle\frac{1}{2}}$}$. The coefficients are
given explicitly by

$c_{j}=\frac{1}{j!}(-p)_{j}(-p+\mbox{${\textstyle\frac{1}{2}}$})_{j}=\frac{(-2p)_{2j}}{2^{2j}j!}.$
(3.6)

We observe that $c_{j}=0$ for $j>p$ so that the above sum of gamma functions
terminates and so is exact. Substitution of the expansion (3.5) in (3.3),
combined with the integral representation [3, (13.4.4)]

$U(\mu,1+\mu-p+\mbox{${\textstyle\frac{1}{2}}$}s,\lambda)=\frac{1}{\Gamma(\mu)}\int_{0}^{\infty}e^{-\lambda
t}t^{\mu-1}(1+t)^{-p-s/2}dt\qquad(\mu>0),$

then shows that, provided $\mu>0$,

$J(a;\lambda)=\frac{(-)^{p}a^{-\delta}\lambda^{p-\mu}}{2\sqrt{\pi}\Gamma(\mu)}\sum_{k\geq
1}\frac{1}{k}\sum_{j=0}^{p}(-)^{j}c_{j}\hskip 142.26378pt$ $\hskip
85.35826pt\times\int_{0}^{\infty}\frac{e^{-\lambda
t}t^{\mu-1}}{(1+t)^{p}}\biggl{(}\frac{1}{2\pi i}\int_{c-\infty i}^{c+\infty
i}\Gamma(\mbox{${\textstyle\frac{1}{2}}$}s+\vartheta-j)\biggl{(}\frac{X_{k}}{1+t}\biggr{)}^{\\!-s/2}ds\biggr{)}dt.$

The inner integral appearing in $J(a;\lambda)$ can be evaluated by making use
of the well-known result

$\frac{1}{2\pi
i}\int_{L}\Gamma(s+\alpha)z^{-s}ds=z^{\alpha}e^{-z}\qquad(|\arg\,z|<\mbox{${\textstyle\frac{1}{2}}$}\pi),$
(3.7)

where $L$ is a path parallel to the imaginary $s$-axis lying to the right of
all the poles of $\Gamma(s+\alpha)$; see, for example, [7, Section 3.3.1].
Evaluation of the inner integral when $|\arg\,a|<\pi/4$ then produces the
final exact result

$J(a;\lambda)=\frac{(-)^{p}\pi^{\mu}e^{\lambda}}{\Gamma(\mu)}\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!\mu\\!-\\!2p\\!-\\!1/2}\sum_{k\geq 1}k^{2p}e^{-2\pi
ka}\sum_{j=0}^{p}(-)^{j}c_{j}\biggl{(}\frac{\lambda}{\pi^{2}k^{2}a^{2}}\biggr{)}^{\\!j}I_{jk},$

where

$I_{jk}=\int_{0}^{\infty}\frac{t^{\mu-1}e^{-\psi(t)}}{(1+t)^{2p-j+1/2}}\,dt,\qquad\psi(t):=\lambda(1+t)+\frac{X_{k}}{1+t}-2\pi
ka.$ (3.8)

In the special case $\mu=0$, we have from (3.3) (since $U(0,b,z)=1$) that

$J(a;\lambda)=\frac{(-)^{p}a^{2p}}{2\sqrt{\pi}\lambda^{p}}\sum_{k\geq
1}\frac{1}{k}\sum_{j=0}^{p}(-)^{j}c_{j}\cdot\frac{1}{2\pi i}\int_{c-\infty
i}^{c+\infty
i}\Gamma(\mbox{${\textstyle\frac{1}{2}}$}+\vartheta-j)\,X_{k}^{-s/2}\,ds.$

Evaluation of the integral by means of (3.7) then produces

$J(a;\lambda)=(-)^{p}\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!-2p-1/2}\sum_{k\geq
1}k^{2p}e^{-\pi^{2}k^{2}a^{2}/\lambda}\sum_{j=0}^{p}(-)^{j}c_{j}\biggl{(}\frac{\lambda}{\pi^{2}k^{2}a^{2}}\biggr{)}^{\\!j}\qquad(\mu=0).$

Then we have the following theorem:

###### Theorem 1

$\\!\\!\\!.$ For $\mu\geq 0$, $\lambda>0$, $\delta=2\mu-\gamma$ and
$\gamma=2p$, where $p$ is a non-negative integer, we have when
$|\arg\,a|<\pi/4$

$S_{\mu,\gamma}(a;\lambda)=a^{1-\delta}{\cal
H}(1)-\mbox{${\textstyle\frac{1}{2}}$}a^{-\delta}\delta_{0p}+J(a;\lambda),$

where ${\cal H}(1)$ is defined in (2.8) and $\delta_{0p}$ is the Kronecker
symbol. The exponentially small contribution $J(a;\lambda)$ is given exactly
by the double sums

$J(a;\lambda)=\frac{(-)^{p}\pi^{\mu}e^{\lambda}}{\Gamma(\mu)}\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!\mu\\!-\\!2p\\!-\\!1/2}\sum_{k\geq 1}k^{2p}e^{-2\pi
ka}\sum_{j=0}^{p}(-)^{j}c_{j}\biggl{(}\frac{\lambda}{\pi^{2}k^{2}a^{2}}\biggr{)}^{\\!j}I_{jk}\qquad(\mu>0),$
(3.9)

and

$J(a;\lambda)=(-)^{p}\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!-2p-1/2}\sum_{k\geq
1}k^{2p}e^{-\pi^{2}k^{2}a^{2}/\lambda}\sum_{j=0}^{p}(-)^{j}c_{j}\biggl{(}\frac{\lambda}{\pi^{2}k^{2}a^{2}}\biggr{)}^{\\!j}\qquad(\mu=0),$
(3.10)

where the coefficients $c_{j}=(-2p)_{j}/(2^{2j}j!)$ and the integrals $I_{jk}$
are defined in (3.8).

Remark 1. When $\mu=p=0$, we find from (2.8), (3.1) and (3.10) (since
$U(\mbox{${\textstyle\frac{1}{2}}$},\mbox{${\textstyle\frac{3}{2}}$},\lambda)=\lambda^{-1/2}$)
the result

$S_{0,0}(a;\lambda)=\frac{a}{2}\sqrt{\frac{\pi}{\lambda}}-\frac{1}{2}+a\sqrt{\frac{\pi}{\lambda}}\sum_{k\geq
1}e^{-\pi^{2}k^{2}a^{2}/\lambda},$

which is the Poisson-Jacobi formula stated in (1.3).

When $p=1$, we have
$U(\mbox{${\textstyle\frac{3}{2}}$},\mbox{${\textstyle\frac{5}{2}}$},\lambda)=\lambda^{-3/2}$
and

$S_{0,2}(a;\lambda)=\frac{a^{3}\sqrt{\pi}}{4\lambda^{3/2}}-\biggl{(}\frac{\pi
a^{2}}{\lambda}\biggr{)}^{\\!5/2}\sum_{k\geq
1}\biggl{(}k^{2}-\frac{\lambda}{2\pi^{2}a^{2}}\biggr{)}e^{-\pi^{2}k^{2}a^{2}/\lambda}.$

We observe that this last case can also be obtained by differentiation of the
Poisson-Jacobi formula with respect to $\lambda$, since

$S_{0,2p}(a;\lambda)=(-)^{p}a^{2p}\,\frac{\partial^{p}}{\partial\lambda^{p}}S_{0,0}(a;\lambda)\qquad(p\geq
1).$

Remark 2. When $p=-1,-2,\ldots\ $, the expansion (3.5) does not terminate and
becomes an inverse factorial expansion. Then we have the exponentially small
contribution given by

$J(a;\lambda)\sim\frac{(-)^{p}\pi^{\mu}e^{\lambda}}{\Gamma(\mu)}\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!\mu\\!-\\!2p\\!-\\!1/2}\sum_{k\geq 1}k^{2p}e^{-2\pi
ka}\sum_{j=0}^{\infty}(-)^{j}c_{j}\biggl{(}\frac{\lambda}{\pi^{2}k^{2}a^{2}}\biggr{)}^{\\!j}I_{jk}\qquad(\mu>0)$
(3.11)

as $|a|\to\infty$ in $|\arg\,a|<\pi/4$.

4\. Alternative form of expansion for positive integer values of $\mu$

Let $\mu=m$ be a positive integer and $\gamma=2p$. We split ${\cal H}(s)$ into
its two constituent parts given by (2.3) and write
$J(a;\lambda)=J_{1}(a;\lambda)+J_{2}(a;\lambda)$ in an obvious manner.

4.1 Evaluation of $J_{1}(a;\lambda).$ Then we have

$J_{1}(a;\lambda)=-\frac{a^{-\delta}}{2\pi}\sum_{k\geq
1}\frac{1}{k}\cdot\frac{1}{2\pi i}\int_{c-\infty i}^{c+\infty i}{\cal
H}_{1}(-s)\Gamma(1+s)\sin\mbox{${\textstyle\frac{1}{2}}$}\pi s\,(2\pi
ka)^{-s}ds,$

where from (2.4) with an application of Kummer’s transformation

$\displaystyle{\cal H}_{1}(-s)$ $\displaystyle=$
$\displaystyle\frac{(-)^{m-p}\pi
e^{\lambda}}{\Gamma(m)\,\sin\mbox{${\textstyle\frac{1}{2}}$}\pi
s}\,\frac{\Gamma(p-\mbox{${\textstyle\frac{1}{2}}$}s)}{\Gamma(1-m+p-\mbox{${\textstyle\frac{1}{2}}$}s)}\,{}_{1}F_{1}(1-m;1-m+p-\mbox{${\textstyle\frac{1}{2}}$}s;-\lambda)$
$\displaystyle=$ $\displaystyle\frac{(-)^{m-p}\pi
e^{\lambda}}{\Gamma(m)\,\sin\mbox{${\textstyle\frac{1}{2}}$}\pi
s}\,\sum_{r=0}^{m-1}\frac{(1-m)_{r}(-\lambda)^{r}}{r!}\,\frac{\Gamma(p-\mbox{${\textstyle\frac{1}{2}}$}s)}{\Gamma(1-m+p-\mbox{${\textstyle\frac{1}{2}}$}s+r)}$
$\displaystyle=$ $\displaystyle\frac{(-)^{p-1}\pi
e^{\lambda}}{\Gamma(m)\,\sin\mbox{${\textstyle\frac{1}{2}}$}\pi
s}\,\sum_{r=0}^{m-1}\frac{(1-m)_{r}\lambda^{r}}{r!}\,\frac{\Gamma(m-p-r+\mbox{${\textstyle\frac{1}{2}}$}s)}{\Gamma(1-p+\mbox{${\textstyle\frac{1}{2}}$}s)}.$

If we make the change of summation index $r\to m-1-\ell$ and use the fact that
$(1-m)_{m-1-\ell}=(-)^{m-1+\ell}\Gamma(m)/\ell!$, we find

${\cal H}_{1}(-s)=\frac{(-)^{m-p}\pi
e^{\lambda}}{\sin\mbox{${\textstyle\frac{1}{2}}$}\pi
s}\,\sum_{\ell=0}^{m-1}A_{\ell}(1-p+\mbox{${\textstyle\frac{1}{2}}$}s)_{\ell},\qquad
A_{\ell}:=\frac{(-)^{\ell}\lambda^{m-1-\ell}}{(m-1-\ell)!\ell!}.$ (4.1)

The Pochhammer symbol appearing in (4.1) can be written in the form

$(1-p+\mbox{${\textstyle\frac{1}{2}}$}s)_{\ell}=\sum_{r=0}^{\ell}B_{r\ell}(s+1)_{r},$

where

$B_{00}=1,\quad B_{01}=\mbox{${\textstyle\frac{1}{2}}$}-p,\quad
B_{11}=\mbox{${\textstyle\frac{1}{2}}$}$
$B_{02}=(\mbox{${\textstyle\frac{1}{2}}$}-p)_{2},\quad
B_{12}=\mbox{${\textstyle\frac{3}{4}}$}-p,\quad
B_{22}=\mbox{${\textstyle\frac{1}{4}}$},$
$B_{03}=(\mbox{${\textstyle\frac{1}{2}}$}-p)_{3},\quad
B_{13}=\mbox{${\textstyle\frac{3}{8}}$}(5-10p+4p^{2}),\quad
B_{23}=\mbox{${\textstyle\frac{3}{4}}$}(1-p),\quad
B_{33}=\mbox{${\textstyle\frac{1}{8}}$},\ \ldots\ .$ (4.2)

Then we obtain

$J_{1}(a;\lambda)=(-)^{k-p-1}\frac{e^{\lambda}a^{-\delta}}{2}\sum_{\ell=0}^{m-1}\sum_{r=0}^{\ell}A_{\ell}B_{r\ell}\sum_{k\geq
1}\frac{1}{k}\cdot\frac{1}{2\pi i}\int_{c-\infty i}^{c+\infty
i}\Gamma(1+s+r)(2\pi ka)^{-s}ds.$ (4.3)

The integrals appearing in $J_{1}(a;\lambda)$ can be evaluated by (3.7) to
produce the final exact result

$J_{1}(a;\lambda)=(-)^{m-p-1}a^{2p+1}\pi
e^{\lambda}\sum_{r=0}^{m-1}\sum_{\ell=r}^{m-1}A_{\ell}B_{r\ell}\frac{\sigma_{r}(a)(2\pi
a)^{r}}{a^{2m-r}}$ (4.4)

for positive integer $m$. Here we have defined the sums

$\sigma_{r}(a):=e^{2\pi a}\sum_{k\geq 1}k^{r}e^{-2\pi ka},$ (4.5)

which have the evaluations for $0\leq r\leq 3$

$\sigma_{0}(a)=\frac{e^{\pi a}}{2\sinh\pi a},\quad\sigma_{1}(a)=\frac{e^{2\pi
a}}{4\sinh^{2}\pi a},\quad\sigma_{2}(a)=\frac{e^{2\pi a}\cosh\pi
a}{4\sinh^{3}\pi a},$ $\sigma_{3}(a)=\frac{e^{2\pi a}(2+\cosh\pi
a)}{8\sinh^{4}\pi a},\ldots\ .$

Note that $J_{1}(a;\lambda)\equiv 0$ when $m=0$, since ${\cal H}_{1}(s)$
vanishes for these values.

4.2 Evaluation of $J_{2}(a;\lambda).$ We have

$J_{2}(a;\lambda)=-\frac{a^{-\delta}}{2\pi}\sum_{k\geq
1}\frac{1}{k}\cdot\frac{1}{2\pi i}\int_{c-\infty i}^{c+\infty i}{\cal
H}_{2}(-s)\Gamma(1+s)\sin\mbox{${\textstyle\frac{1}{2}}$}\pi s\,(2\pi
ka)^{-s}ds,$

where from (2.5)

$\displaystyle{\cal H}_{2}(s)$ $\displaystyle=$
$\displaystyle(-)^{m-p-1}\frac{\pi\lambda^{m-p-s/2}}{\sin\mbox{${\textstyle\frac{1}{2}}$}\pi
s}\,\frac{{}_{1}F_{1}(m;1+m-p+\mbox{${\textstyle\frac{1}{2}}$}s;\lambda)}{\Gamma(1+m-p+\mbox{${\textstyle\frac{1}{2}}$}s)}$
$\displaystyle=$
$\displaystyle(-)^{m-p-1}\frac{\pi\lambda^{m-p-s/2}}{\sin\mbox{${\textstyle\frac{1}{2}}$}\pi
s}\,\sum_{r=0}^{\infty}\frac{(m)_{r}\lambda^{r}}{r!\Gamma(1\\!+\\!m\\!-\\!p\\!+\\!r\\!+\\!\mbox{${\textstyle\frac{1}{2}}$}s)}.$

Then we obtain

$J_{2}(a;\lambda)=\frac{(-\lambda)^{m-p}a^{-\delta}}{\sqrt{\pi}}\sum_{k\geq
1}\frac{1}{k}\sum_{r=0}^{\infty}\frac{(m)_{r}\lambda^{r}}{r!}\,K_{kr},$ (4.6)

where

$K_{kr}=\frac{\sqrt{\pi}}{4\pi i}\int_{c-\infty i}^{c+\infty
i}\frac{2^{-s}\Gamma(1+s)}{\Gamma(1\\!+\\!m\\!-\\!p\\!+\\!r\\!+\\!\mbox{${\textstyle\frac{1}{2}}$}s)}\,X_{k}^{-s/2}ds$
(4.7)

with $X_{k}$ defined in (3.4). Since there no poles of the integrand in
$\Re(s)>0$ the integration path in (4.7) can be displaced as far to the right
as we please such that $|s|$ is everywhere large on the new path. The quotient
of gamma functions in the integrand can then be expanded in a manner similar
to that in (3.5) to find [7, p. 53]

$\frac{2^{-s}\sqrt{\pi}\Gamma(1+s)}{\Gamma(1\\!+\\!m\\!-\\!p\\!+\\!r\\!+\\!\mbox{${\textstyle\frac{1}{2}}$}s)}=\frac{\Gamma(\mbox{${\textstyle\frac{1}{2}}$}+\mbox{${\textstyle\frac{1}{2}}$}s)\Gamma(1+\mbox{${\textstyle\frac{1}{2}}$}s)}{\Gamma(1\\!+\\!m\\!-\\!p\\!+\\!r\\!+\\!\mbox{${\textstyle\frac{1}{2}}$}s)}\hskip
113.81102pt$ $\hskip
85.35826pt=\sum_{j=0}^{M-1}(-)^{j}{\hat{c}}_{j}(r)\Gamma(\mbox{${\textstyle\frac{1}{2}}$}+\vartheta-j)+\rho_{M}(s)\Gamma(\mbox{${\textstyle\frac{1}{2}}$}s+\vartheta-M),$

where $M$ is a positive integer,
$\vartheta=\mbox{${\textstyle\frac{1}{2}}$}+p-m-r$ and $\rho_{M}(s)=O(1)$ as
$|s|\to\infty$ in $|\arg\,s|<\pi$. The coefficients are given explicitly by

${\hat{c}}_{j}(r)=\frac{1}{j!}(m-p+r)_{j}(m-p+r+\mbox{${\textstyle\frac{1}{2}}$})_{j}=\frac{(2m-2p+2r)_{2j}}{2^{2j}j!}.$
(4.8)

When $m=r=0$, the coefficients ${\hat{c}}_{j}(r)$ reduce to $c_{j}$ in (3.6).
Then

$K_{kr}=\sum_{j=0}^{M-1}(-)^{r}c_{j}(r)\cdot\frac{1}{4\pi i}\int_{c-\infty
i}^{c+\infty
i}\Gamma(\mbox{${\textstyle\frac{1}{2}}$}s+\vartheta-j)X_{k}^{-s/2}ds+R_{M,r}$
$=X_{k}^{\vartheta}e^{-X_{k}}\biggl{\\{}\sum_{j=0}^{M-1}(-)^{j}c_{j}(r)X_{k}^{-j}+O(X_{k}^{-M})\biggr{\\}}$
(4.9)

by (3.7), where the remainder term

$R_{M,r}=\frac{1}{4\pi i}\int_{c-\infty i}^{c+\infty
i}\rho_{M}(s)\Gamma(\mbox{${\textstyle\frac{1}{2}}$}s+\vartheta-j)X_{k}^{-s/2}ds=O(X_{k}^{\vartheta-M}e^{-X_{k}})$

as $|a|\to\infty$ in $|\arg\,a|<\pi/4$ by Lemma 2.7 in [7, p. 71].

From (4.6) and (4.9), it then follows that

$J_{2}(a;\lambda)\sim(-)^{m-p}\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!\delta-1/2}\sum_{k\geq
1}\frac{e^{-\pi^{2}k^{2}a^{2}/\lambda}}{k^{\delta}}\sum_{r=0}^{\infty}\sum_{j=0}^{\infty}(-)^{j}c_{j}(r)\,\frac{(m)_{r}\lambda^{r}}{r!}\biggl{(}\frac{\lambda}{\pi^{2}k^{2}a^{2}}\biggr{)}^{\\!r+j}$
$=(-)^{m-p}\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!\delta-1/2}\sum_{r=0}^{\infty}(-)^{r}C_{r}\biggl{(}\frac{\lambda}{\pi^{2}a^{2}}\biggr{)}^{\\!r}\sum_{k\geq
1}\frac{e^{-\pi^{2}k^{2}a^{2}/\lambda}}{k^{2r+\delta}},$ (4.10)

where the double sum over $r$ and $j$ has been summed ‘diagonally’ (see [8, p.
58]) and the coefficients $C_{r}$ are given by

$C_{r}\equiv
C_{r}(m,p,\lambda)=\sum_{n=0}^{r}\frac{(-\lambda)^{n}(m)_{n}}{n!}\,{\hat{c}}_{r-n}(n).$

Use of (4.8) shows that the $C_{r}$ may be expressed in terms of a terminating
${}_{2}F_{2}$ hypergeometric series:

$C_{r}=\sum_{n=0}^{r}\frac{(-\lambda)^{r}(m)_{n}}{n!(r-n)!}\,\frac{2^{2n-2r}(2m-2p)_{2r}}{(2m-2p)_{2n}}=\frac{(2m-2p)_{2r}}{2^{2r}r!}\sum_{n=0}^{r}\frac{(-r)_{n}(m)_{n}\lambda^{n}}{n!(m-p)_{n}(m-p+\mbox{${\textstyle\frac{1}{2}}$})_{n}}$
$=\frac{(2m-2p)_{2r}}{2^{2r}r!}\,{}_{2}F_{2}(-r,m;m-p,m-p+\mbox{${\textstyle\frac{1}{2}}$};\lambda),$
(4.11)

where we have made use of the results $(a+2n)_{2r-2n}=(a)_{2r}/(a)_{2n}$ and
$(2a)_{2n}=2^{2n}(a)_{n}(a+\mbox{${\textstyle\frac{1}{2}}$})_{n}$.

Then we have the following theorem:

###### Theorem 2

$\\!\\!\\!.$ Let $\mu=m$ and $\gamma=2p$, where $m\geq 1$, $p$ are integers,
and $\delta=2(m-p)$ with $\lambda>0$. Then we have the representation

$S_{\mu,\gamma}(a;\lambda)=a^{1-\delta}{\cal
H}(1)+H_{\mu,\gamma}(a;\lambda)+J(a;\lambda),$ (4.12)

where ${\cal H}(1)$ and $H_{\mu,\gamma}(a;\lambda)$ are defined in (2.8) and
(3.2). The exponentially small contribution $J(a;\lambda)$ has the expansion

$J(a;\lambda)+(-)^{m-p}\pi a^{2p+1}e^{\lambda-2\pi
a}\sum_{r=0}^{m-1}\sum_{\ell=r}^{m-1}A_{\ell}B_{r\ell}\,\frac{\sigma_{r}(a)(2\pi)^{r}}{a^{2m-r}}\hskip
113.81102pt$ $\hskip 85.35826pt\sim(-)^{m-p}\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!\delta-1/2}\sum_{r=0}^{\infty}(-)^{r}C_{r}\biggl{(}\frac{\lambda}{\pi^{2}a^{2}}\biggr{)}^{\\!r}\sum_{k\geq
1}\frac{e^{-\pi^{2}k^{2}a^{2}/\lambda}}{k^{2r+\delta}}$ (4.13)

as $|a|\to\infty$ in $|\arg\,a|<\pi/4$. The coefficients $A_{\ell}$,
$B_{r\ell}$ and $C_{r}$ are defined in (4.1), (4.2) and (4.11) and the
functions $\sigma_{r}(a)$ are given by (4.5). When $m=0$ the double sum on the
left-hand side of (4.13) vanishes.

5\. Examples

We present some examples of the expansion of $S_{\mu,\gamma}(a;\lambda)$
stated in Theorem 2.

Example 1. Let $\mu=0$ and $\gamma=2p$. Then from (4.12) and (4.13) we have,
for $p=1,2,\ldots\ $,

$S_{0,2p}(a;\lambda)=\sum_{n=1}^{\infty}n^{2p}e^{-\lambda
n^{2}/a^{2}}\sim\frac{1}{2}\biggl{(}\frac{\lambda}{a^{2}}\biggr{)}^{-p-1/2}\Gamma(p+\mbox{${\textstyle\frac{1}{2}}$})\hskip
85.35826pt$ $+(-)^{p}\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!-2p-1/2}\sum_{r=0}^{\infty}(-)^{r}C_{r}\biggl{(}\frac{\lambda}{\pi^{2}a^{2}}\biggr{)}^{\\!r}\sum_{k\geq
1}\frac{e^{-\pi^{2}k^{2}a^{2}/\lambda}}{k^{2r-2p}}$ (5.1)

as $|a|\to\infty$ in $|\arg\,a|<\pi/4$, where from (4.11)
$C_{r}=(-2p)_{2r}/(2^{2r}r!)$. The case $p=0$ is covered in Remark 1.

In the case $p\leq-1$ we let $p=-q$ and note that ${\cal
H}(s)=\mbox{${\textstyle\frac{1}{2}}$}\lambda^{q-s/2}\Gamma(\mbox{${\textstyle\frac{1}{2}}$}s-q)$
when $\mu=0$. The residue of ${\cal H}(s)$ at $s=2q-2k$ ($k=0,1,2,\ldots$) is
given by $(-)^{k}/k!$. Then

$S_{0,-2q}(a;\lambda)=\sum_{n=1}^{\infty}\frac{e^{-\lambda
n^{2}/a^{2}}}{n^{2q}}\sim\frac{1}{2}\biggl{(}\frac{\lambda}{a^{2}}\biggr{)}^{\\!q-1/2}\Gamma(\mbox{${\textstyle\frac{1}{2}}$}-q)+\sum_{k=0}^{q}\frac{(-)^{k}}{k!}\biggl{(}\frac{\lambda}{a^{2}}\biggr{)}^{\\!k}\zeta(2q-2k)$
$+(-)^{q}\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!2q-1/2}\sum_{r=0}^{\infty}\frac{(-)^{r}(2q)_{2r}}{2^{2r}r!}\biggl{(}\frac{\lambda}{\pi^{2}a^{2}}\biggr{)}^{\\!r}\sum_{k\geq
1}\frac{e^{-\pi^{2}k^{2}a^{2}/\lambda}}{k^{2q+2r}}$ (5.2)

as $|a|\to\infty$ in $|\arg\,a|<\pi/4$.

Example 2. Let $\mu=1$ and $\gamma=0$ (so that $\delta=2$). Then we have

$S_{1,0}(a;\lambda)=\sum_{n=1}^{\infty}\frac{e^{-\lambda
n^{2}/a^{2}}}{n^{2}+a^{2}}=\frac{\pi
e^{\lambda}}{2a}\mbox{erfc}\sqrt{\lambda}-\frac{1}{2a^{2}}+J(a;\lambda),$
(5.3)

since [3, (13.6.8)]

${\cal
H}(1)=\mbox{${\textstyle\frac{1}{2}}$}\sqrt{\pi}\,U(\mbox{${\textstyle\frac{1}{2}}$},\mbox{${\textstyle\frac{1}{2}}$},\lambda)=\mbox{${\textstyle\frac{1}{2}}$}\pi
e^{\lambda}\,\mbox{erfc}\sqrt{\lambda},$

where erfc is the complementary error function. From (4.13), the exponentially
small contribution is

$J(a;\lambda)-\frac{\pi e^{\lambda}}{2a}\,\frac{e^{-\pi a}}{\sinh\pi
a}\sim\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!3/2}\sum_{r=0}^{\infty}(-)^{r}C_{r}\biggl{(}\frac{\lambda}{\pi^{2}a^{2}}\biggr{)}^{\\!r}\sum_{k\geq
1}\frac{e^{-\pi^{2}k^{2}a^{2}/\lambda}}{k^{2r+2}}$ (5.4)

as $|a|\to\infty$ in $|\arg\,a|<\pi/4$, where from (4.11) the coefficients
$C_{r}$ are given by

$C_{r}=\frac{2\Gamma(r+\mbox{${\textstyle\frac{3}{2}}$})}{\sqrt{\pi}}\,{}_{1}F_{1}(-r;\mbox{${\textstyle\frac{3}{2}}$};\lambda).$

The first few $C_{r}$ are therefore

$C_{0}=1,\quad C_{1}=\mbox{${\textstyle\frac{3}{2}}$}-\lambda,\quad
C_{2}=\mbox{${\textstyle\frac{15}{4}}$}-5\lambda+\lambda^{2},$
$C_{3}=\mbox{${\textstyle\frac{105}{8}}$}-\mbox{${\textstyle\frac{105}{4}}$}\lambda+\mbox{${\textstyle\frac{21}{2}}$}\lambda^{2}-\lambda^{3}.$

We remark that the case $\mu=1$, $\gamma=2$ can be obtained directly from
Theorem 2, but also follows from the Poisson-Jacobi formula (1.3), together
with (5.3) and (5.4), since

$\displaystyle S_{1,2}(a;\lambda)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{n^{2}e^{-\lambda
n^{2}/a^{2}}}{n^{2}+a^{2}}=\sum_{n=1}^{\infty}\biggl{(}1-\frac{a^{2}}{n^{2}+a^{2}}\biggr{)}e^{-\lambda
n^{2}/a^{2}}$ $\displaystyle=$ $\displaystyle
S_{0,0}(a;\lambda)-a^{2}S_{1,0}(a;\lambda).$

Example 3. Let $\mu=2$ and $\gamma=0$ (so that $\delta=4$). Then we have

$S_{2,0}(a;\lambda)=\sum_{n=1}^{\infty}\frac{e^{-\lambda
n^{2}/a^{2}}}{(n^{2}+a^{2})^{2}}=\frac{\sqrt{\pi}}{2a^{3}}\,U(\mbox{${\textstyle\frac{1}{2}}$},-\mbox{${\textstyle\frac{1}{2}}$},\lambda)-\frac{1}{2a^{4}}+J(a;\lambda).$
(5.5)

From (4.1) and (4.2), the coefficients $A_{0}=\lambda$, $A_{1}=-1$,
$B_{00}=1$, $B_{01}=B_{11}=\mbox{${\textstyle\frac{1}{2}}$}$, so that

$J_{1}(a;\lambda)=-\frac{\pi e^{\lambda-\pi a}}{2a^{3}\sinh\pi
a}\biggl{\\{}\lambda-\frac{1}{2}-\frac{\pi ae^{\pi a}}{2\sinh\pi
a}\biggr{\\}}$

and

$J_{2}(a;\lambda)\sim\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!7/2}\sum_{r=0}^{\infty}(-)^{r}C_{r}\biggl{(}\frac{\lambda}{\pi^{2}a^{2}}\biggr{)}^{\\!r}\sum_{k\geq
1}\frac{e^{-\pi^{2}k^{2}a^{2}/\lambda}}{k^{2r+4}}.$

From (4.11) the coefficients $C_{r}$ are given by

$C_{r}=\frac{4(r+1)}{3\sqrt{\pi}}\,\Gamma(r+\mbox{${\textstyle\frac{5}{2}}$})\,{}_{1}F_{1}(-r;\mbox{${\textstyle\frac{5}{2}}$};\lambda)$

so that the first few coefficients are therefore

$C_{0}=1,\quad C_{1}=5-2\lambda,\quad
C_{2}=\mbox{${\textstyle\frac{105}{4}}$}-21\lambda+3\lambda^{2},$
$C_{3}=\mbox{${\textstyle\frac{315}{2}}$}-189\lambda+54\lambda^{2}-4\lambda^{3}.$

Then the exponentially small contribution is

$J(a;\lambda)+\frac{\pi e^{\lambda-\pi a}}{2a^{3}\sinh\pi
a}\biggl{\\{}\lambda-\frac{1}{2}-\frac{\pi ae^{\pi a}}{2\sinh\pi
a}\biggr{\\}}\hskip 142.26378pt$ $\hskip
113.81102pt\sim\biggl{(}\frac{\lambda}{\pi
a^{2}}\biggr{)}^{\\!7/2}\sum_{r=0}^{\infty}(-)^{r}C_{r}\biggl{(}\frac{\lambda}{\pi^{2}a^{2}}\biggr{)}^{\\!r}\sum_{k\geq
1}\frac{e^{-\pi^{2}k^{2}a^{2}/\lambda}}{k^{2r+4}}$ (5.6)

as $|a|\to\infty$ in $|\arg\,a|<\pi/4$.

As mentioned in the previous example, the sums with $\mu=2$ and $p=2,4$ can be
obtained directly from Theorem 2, but also from the identities

$\displaystyle S_{2,1}(a;\lambda)$ $\displaystyle=$ $\displaystyle
S_{1,0}(a;\lambda)-a^{2}S_{2,0}(a;\lambda)$ $\displaystyle S_{2,2}(a;\lambda)$
$\displaystyle=$ $\displaystyle
S_{0,0}(a;\lambda)-2a^{2}S_{2,1}(a;\lambda)-3a^{4}S_{2,0}(a;\lambda).$

6\. Numerical results and concluding remarks

The expansion of the exponentially small contribution (when $\gamma=2p$) given
in Theorem 1 is exact for $\mu\geq 0$. It is possible to employ an asymptotic
expansion for the integrals $I_{jk}$, but this would necessarliy introduce an
error. However, we have evaluated these integrals to high numerical precision
and have thereby verified the expansion (3.10) of $J(a;\lambda)$ for several
parameter values to 50 decimal precision.

We present some numerical examples of the large-$a$ expansion of
$S_{\mu,\gamma}(a;\lambda)$ given in Theorem 2 to demonstrate the accuracy of
our results. We subtract from the sum $S_{\mu,\gamma}(a;\lambda)$ the finite
terms appearing in (4.12) by defining

${\hat{S}}_{\mu,\gamma}(a;\lambda):=S_{\mu,\gamma}(a;\lambda)-\\{a^{1-\delta}{\cal
H}(1)+H_{\mu,\gamma}(a;\lambda)+J_{1}(a;\lambda)\\}$ (6.1)

and comparing it with the exponentially small asymptotic expansion
$J_{2}(a;\lambda)$ in (4.10). We stress that the contribution
$J_{1}(a;\lambda)$ in (4.4) is an exact result when $\mu$ is an integer. In
Table 1 we show the values of the absolute relative error in the high-
precision computation of ${\hat{S}}_{\mu,\gamma}(a;\lambda)$ from (1.1) using
the asymptotic expansion for $J_{2}(a;\lambda)$ for different truncation index
$r$. The values of $\mu$ and $\gamma$ chosen correspond to the examples given
in Section 5. The final entry in each column gives the value of
${\hat{S}}_{\mu,\gamma}(a;\lambda)$. It is seen that the exponentially small
contribution to $S_{\mu,\gamma}(a;\lambda)$ when $\gamma$ is an even integer
agrees well with the expansion given in Theorem 2.

Table 1: The absolute relative error in the computation of ${\hat{S}}_{\mu,\gamma}(a;\lambda)$ from (6.1) for different $\mu$, $\gamma$ and truncation index $r$ in the asymptotic expansion $J_{2}(a;\lambda)$ when $\lambda=2$ and $a=3$. $r$ | $\mu=0,\ \gamma=-2$ | $\mu=1,\ \gamma=0$ | $\mu=2,\ \gamma=0$
---|---|---|---
0 | $3.307\times 10^{-02}$ | $1.007\times 10^{-02}$ | $2.464\times 10^{-02}$
1 | $1.823\times 10^{-03}$ | $1.070\times 10^{-03}$ | $1.572\times 10^{-03}$
2 | $1.408\times 10^{-04}$ | $5.900\times 10^{-05}$ | $3.764\times 10^{-04}$
5 | $2.438\times 10^{-07}$ | $7.124\times 10^{-08}$ | $3.325\times 10^{-07}$
10 | $9.421\times 10^{-11}$ | $3.101\times 10^{-12}$ | $7.198\times 10^{-10}$
15 | $3.138\times 10^{-13}$ | $6.596\times 10^{-14}$ | $2.317\times 10^{-12}$
20 | $4.678\times 10^{-15}$ | $4.335\times 10^{-16}$ | $5.392\times 10^{-14}$
${\hat{S}}_{\mu,\gamma}$ | $-9.3737097\times 10^{-22}$ | $-9.7822227\times 10^{-22}$ | $+4.7287147\times 10^{-24}$

It is worth mentioning that $J(a;\lambda)$ given in Theorems 1 and 2 appears
to comprise two different types of exponentially small terms, namely
$\exp(-2\pi ka)$ in Theorem 1 and both $\exp(-2\pi ka)$ and
$\exp(-\pi^{2}k^{2}a^{2}/\lambda)$, $k\geq 1$ in Theorem 2. However, a closer
examination of the integrals $I_{jk}$ appearing in Theorem 1 reveals that they
also contain the more subdominant terms $\exp(-\pi^{2}k^{2}a^{2}/\lambda)$. To
see this we consider

$e^{-2\pi ka}I_{jk}=e^{-2\pi
ka}\int_{0}^{\infty}\frac{t^{\mu-1}e^{-\psi(t)}}{(1+t)^{\beta}}dt,\qquad\beta=2p-j+\mbox{${\textstyle\frac{1}{2}}$}.$

The phase function $\psi(t)$ in (3.8) has a saddle point at $t=t_{s}$, where
$1+t_{s}=\sqrt{X_{k}/\lambda}=\pi ka/\lambda$ and $\psi(t_{s})=0$,
$\psi^{\prime\prime}(t_{s})=2X_{k}/(1+t_{s})^{3}=2\lambda^{2}/(\pi ka)$. For
large complex $a$ in the sector $|\arg\,a|<\pi/4$, the integration path is
chosen to emanate from the origin in the direction $\arg\,t=\phi=\pi-2\arg\,a$
to the singularity at $t=-1$ and thence along the path of steepest descent
through $t_{s}$ to infinity in $\Re(t)>0$. The contribution from the saddle is
controlled by

$2e^{-2\pi
ka}\sqrt{\frac{\pi}{2\psi^{\prime\prime}(t_{s})}}\,\frac{t_{s}^{\mu-1}}{(1+t_{s})^{\beta}}=\sqrt{\frac{\pi}{\lambda}}\biggl{(}\frac{\pi
ka}{\lambda}\biggr{)}^{\\!\mu-\beta-1/2}e^{-2\pi ka}$

while that from the neighbourhood of the origin is approximately

$e^{-X_{k}+i\mu\phi}\int_{0}^{\infty}e^{-|X_{k}|\tau}\tau^{\mu-1}d\tau=O(X_{k}^{-\mu}e^{-X_{k}}),$

which produces the more subdominant exponential terms.

Finally, we note that the alternating version of (1.1) can be expressed in
terms of $S_{\mu,\gamma}(a;\lambda)$ since

$\sum_{n=1}^{\infty}\frac{(-)^{n-1}n^{\gamma}}{(n^{2}+a^{2})^{\mu}}\,e^{-\lambda
n^{2}/a^{2}}=S_{\mu,\gamma}(a;\lambda)-2^{1-\delta}S_{\mu,\gamma}(\mbox{${\textstyle\frac{1}{2}}$}a;\lambda).$

Application of Theorems 1 and 2 then enables the large-$a$ expansion of the
alternating series to be determined.

Appendix A: The pole structure of ${\cal H}(s)$

The function ${\cal H}(s)$ defined in (2.3), (2.4) and (2.5) has pokes at
$s=-2k-\gamma$, $k=0,1,2,\ldots\ $ and apparent poles at $s=\pm 2k+\delta$,
$\delta=2\mu-\gamma$. We shall show in this appendix that ${\cal H}(s)$ is
regular at these last points. We have

${\cal H}(s)=\frac{\pi
G(s)}{2\sin\pi(\mu-\mbox{${\textstyle\frac{\gamma+s}{2}}$})},$

where

$G(s):=\frac{\Gamma(\mbox{${\textstyle\frac{\gamma+s}{2}}$})}{\Gamma(\mu)}{\bf
F}(\mbox{${\textstyle\frac{\gamma+s}{2}}$};1\\!+\\!\mbox{${\textstyle\frac{\gamma+s}{2}}$}\\!-\\!\mu;\lambda)-\lambda^{\mu-(\gamma+s)/2}{\bf
F}(\mu;1\\!-\\!\mbox{${\textstyle\frac{\gamma+s}{2}}$}\\!+\\!\mu;\lambda).$

Here ${\bf F}$ denotes the normalised confluent hypergeometric function
defined by

${\bf F}(a;b,z)=\frac{1}{\Gamma(b)}\,{}_{1}F_{1}(a;b;z),$

which is defined for all values of the parameter $b$.

Let $s_{k}=2k+\delta$ so that
$\mbox{${\textstyle\frac{1}{2}}$}(\gamma+s_{k})=k+\mu$. Then

$\displaystyle G(s_{k})$ $\displaystyle=$ $\displaystyle(\mu)_{k}{\bf
F}(\mu+k;k+1;\lambda)-\lambda^{-k}\sum_{r=k}^{\infty}\frac{(\mu)_{r}\lambda^{r}}{r!\Gamma(1+r-k)}$
$\displaystyle=$ $\displaystyle(\mu)_{k}{\bf
F}(\mu+k;k+1;\lambda)-\sum_{r=0}^{\infty}\frac{(\mu)_{r+k}\lambda^{r}}{r!\Gamma(1+k+r)}$
$\displaystyle=$ $\displaystyle(\mu)_{k}{\bf
F}(\mu+k;k+1;\lambda)-(\mu)_{k}\sum_{r=0}^{\infty}\frac{(\mu+k)_{r}\lambda^{r}}{r!\Gamma(1+k+r)}\equiv
0.$

Hence ${\cal H}(s)$ is regular at $s_{k}=2k+\delta$.

A similar argument when $s_{k}=-2k+\delta$ shows that

$\displaystyle G(s_{k})$ $\displaystyle=$
$\displaystyle\frac{\Gamma(\mu-k)}{\Gamma(\mu)}\sum_{r=k}^{\infty}\frac{(\mu-k)_{r}\lambda^{r}}{r!\Gamma(1-k+r)}-\lambda^{k}{\bf
F}(\mu;1+k;\lambda)$ $\displaystyle=$
$\displaystyle\frac{\Gamma(\mu-k)}{\Gamma(\mu)}\sum_{r=0}^{\infty}\frac{(\mu-k)_{r+k}\lambda^{r+k}}{r!\Gamma(1+k+r)}-\lambda^{k}{\bf
F}(\mu;1+k;\lambda)$ $\displaystyle=$
$\displaystyle\lambda^{k}\sum_{r=0}^{\infty}\frac{(\mu)_{r}\lambda^{r}}{r!\Gamma(1+k+r)}-\lambda^{k}{\bf
F}(\mu;1+k;\lambda)\equiv 0,$

so that ${\cal H}(s)$ is also regular at the points $s_{k}=-2k+\delta$.

Appendix B: The expansion in the case $\gamma=-1$

We consider the large-$a$ expansion of $S_{\mu,\gamma}(a;\lambda)$ given in
(2.7) in the special case $\gamma=-1$ when the singularity of the integrand in
(2.6) at $s=1$ is a double pole. We set $s=1+\epsilon$, with $\epsilon\to 0$.
Then

${\cal
H}_{1}(s)=\frac{\Gamma(\mbox{${\textstyle\frac{1}{2}}$}\epsilon)\Gamma(\mu-\mbox{${\textstyle\frac{1}{2}}$}\epsilon)}{2\Gamma(\mu)}\,{}_{1}F_{1}(\mbox{${\textstyle\frac{1}{2}}$}\epsilon;1-\mu+\mbox{${\textstyle\frac{1}{2}}$}\epsilon;\lambda)$
$=\frac{1}{\epsilon}\biggl{\\{}1+\mbox{${\textstyle\frac{1}{2}}$}\epsilon(\psi(1)-\psi(\mu))+O(\epsilon^{2})\biggr{\\}}\,{}_{1}F_{1}(\mbox{${\textstyle\frac{1}{2}}$}\epsilon;1-\mu+\mbox{${\textstyle\frac{1}{2}}$}\epsilon;\lambda),$

where

${}_{1}F_{1}(\mbox{${\textstyle\frac{1}{2}}$}\epsilon;1-\mu+\mbox{${\textstyle\frac{1}{2}}$}\epsilon;\lambda)=1+\frac{\epsilon\lambda}{2(1-\mu)}\biggl{\\{}1+\frac{1!\lambda}{(2-\mu)2!}+\frac{2!\lambda^{2}}{(2-\mu)_{2}3!}+\cdots\biggr{\\}}+O(\epsilon^{2})$
$=1+\frac{\epsilon\lambda}{2(1-\mu)}\,{}_{2}F_{2}(1,1;2,2-\mu;\lambda)+O(\epsilon^{2}).$

Using the fact that
$\zeta(1+\epsilon)=\epsilon^{-1}(1+\gamma_{{}_{E}}\epsilon+O(\epsilon^{2}))$
and $\psi(1)=-\gamma_{{}_{E}}$, where $\gamma_{{}_{E}}$ is the Euler-
Mascheroni constant, we obtain the residue resulting from ${\cal H}_{1}(s)$ at
$s=1$ given by

$a\biggl{\\{}\log\,a+\frac{1}{2}\gamma_{{}_{E}}-\frac{1}{2}\psi(\mu)+\frac{\lambda}{2(1-\mu)}\,{}_{2}F_{2}(1,1;2,2-\mu;\lambda)\biggr{\\}}\qquad(\mu\neq
1,2,\ldots).$

The residue resulting from ${\cal H}_{2}(s)$ is

$a{\cal
H}_{2}(1)=\frac{1}{2}\lambda^{\mu}\Gamma(1-\mu)\,{}_{1}F_{1}(\mu;1+\mu;\lambda)\qquad(\mu\neq
1,2,\ldots).$

Hence, provided $\mu\neq 1,2,\ldots\ $,

$S_{\mu,-1}(a;\lambda)=\sum_{n=1}^{\infty}\frac{e^{-\lambda
n^{2}/a^{2}}}{n(n^{2}+a^{2})^{\mu}}\sim
a^{1-\delta}\biggl{\\{}\log\,a+\frac{1}{2}\gamma_{{}_{E}}-\frac{1}{2}\psi(\mu)+\frac{\lambda}{2(1-\mu)}\,{}_{2}F_{2}(1,1;2,2-\mu;\lambda)\biggr{\\}}$
$+a^{1-\delta}\sum_{k=1}^{\infty}\frac{(-)^{k}(\mu)_{k}}{k!}\zeta(1-2k)\,{}_{1}F_{1}(-k;1-\mu-k;\lambda)a^{-2k}$
(B.1)

as $a\to\infty$ in $|\arg\,a|<\pi/4$.

When $\mu$ is a positive integer a limiting process is required. To
illustrate, we consider only the case $\mu=1$. We find that ${\cal
H}_{1}(1+\epsilon)=e^{\lambda}/\epsilon$ and

${\cal
H}_{2}(1+\epsilon)=-\frac{\lambda}{\epsilon(1-\mbox{${\textstyle\frac{1}{2}}$}\epsilon)}\Gamma(1+\mbox{${\textstyle\frac{1}{2}}$}\epsilon)\,{}_{1}F_{1}(1;2-\mbox{${\textstyle\frac{1}{2}}$}\epsilon;\lambda)$
$=-\frac{\lambda}{\epsilon}\biggl{\\{}1+\frac{1}{2}\epsilon(1-\gamma_{{}_{E}}-\log\,\lambda)+O(\epsilon^{2})\biggr{\\}}\,{}_{1}F_{1}(1;2-\mbox{${\textstyle\frac{1}{2}}$}\epsilon;\lambda),$

where

$\displaystyle{}_{1}F_{1}(1;2-\mbox{${\textstyle\frac{1}{2}}$}\epsilon;\lambda)\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!1\\!+\\!\frac{\lambda}{2}\biggl{(}1\\!+\\!\frac{\epsilon}{4}\biggr{)}\\!+\\!\frac{\lambda^{2}}{2\cdot
3}\biggl{(}1\\!+\\!\frac{\epsilon}{4}\\!+\\!\frac{\epsilon}{6}\biggr{)}+\cdots+O(\epsilon^{2})$
$\displaystyle=$
$\displaystyle\\!\\!\\!{}_{1}F_{1}(1;2;\lambda)+\frac{\epsilon}{2}\biggl{\\{}\frac{\lambda}{2}\frac{1}{2}\\!+\\!\frac{\lambda^{2}}{2\cdot
3}\biggl{(}\frac{1}{2}\\!+\\!\frac{1}{3}\biggr{)}\\!+\\!\frac{\lambda^{3}}{2\cdot
3\cdot
4}\biggl{(}\frac{1}{2}\\!+\\!\frac{1}{3}\\!+\\!\frac{1}{4}\biggr{)}\\!+\\!\cdots\biggr{\\}}+O(\epsilon^{2})$
$\displaystyle=$
$\displaystyle\\!\\!\\!\frac{e^{\lambda}-1}{\lambda}\\!+\\!\frac{\epsilon}{2}\sum_{n=1}^{\infty}\frac{\lambda^{n}\tau(n)}{(2)_{n}}+O(\epsilon^{2}),\qquad\tau(n):=\sum_{r=1}^{n}\frac{1}{r+1}.$

Then we obtain the expansion

$S_{1,-1}(a;\lambda)=\sum_{n=1}^{\infty}\frac{e^{-\lambda
n^{2}/a^{2}}}{n(n^{2}+a^{2})}\sim\frac{a^{1-\delta}}{2}\biggl{\\{}e^{\lambda}(\log\,\lambda+\gamma_{{}_{E}}-1)+\log\,(a^{2}/\lambda)+\gamma_{{}_{E}}+1-\lambda\sum_{n=1}^{\infty}\frac{\lambda^{n}\tau(n)}{(2)_{n}}\biggr{\\}}$
$+a^{1-\delta}\sum_{k=1}^{\infty}\frac{(-)^{k}(\mu)_{k}}{k!}\zeta(1-2k)e_{k}(\lambda)a^{-2k}$
(B.2)

as $a\to\infty$ in $|\arg\,a|<\pi/4$, where

$e_{k}(\lambda):={}_{1}F_{1}(-k;-k;\lambda)=\sum_{n=0}^{k}\frac{\lambda^{n}}{n!}.$

## References

  * [1] S. Gerhold and Z. Tomovski, Asymptotic expansion of Mathieu power series and trigonometric Mathieu series, 2019. [arXiv:1906.02055].
  * [2] E.L. Mathieu, Traité de Physique Mathématique. VI–VII: Théorie de l’Elasticité des Corps Solides (Part 2), Gauthier-Villars, Paris, 1890.
  * [3] F.W.J. Olver, D.W. Lozier, R.F. Boisvert and C.W. Clark (eds.), NIST Handbook of Mathematical Functions, Cambridge University Press, Cambridge, 2010.
  * [4] R.B. Paris, The asymptotic expansion of a generalised Mathieu series, Appl. Math. Sci. 125 (2013) 6209–6216.
  * [5] R.B. Paris, Exponentially small expansions associated with a generalised Mathieu series, Math. Aeterna 6 (2016) 181–196.
  * [6] R.B. Paris, Asymptotic expansions of Mathieu-Bessel series I, II. 2019. [arXiv:1907.01812, 1909.09805].
  * [7] R.B. Paris and D. Kaminski, Asymptotics and Mellin-Barnes Integrals, Cambridge University Press, Cambridge, 2001.
  * [8] L.J. Slater, Generalized Hypergeometric Functions, Cambridge University Press, Cambridge, 1966.
  * [9] E.C. Titchmarsh, The Theory of the Riemann Zeta-Function, revised by D.R. Heath-Brown, Oxford University Press, Oxford, 1986.
  * [10] E. T. Whittaker and G. N. Watson, Modern Analysis, Cambridge University Press, Cambridge, 1952.
  * [11] V.P. Zastavnyi, Asymptotic expansions of several series and their application, Ukrainian Math. Bull. 6 (2009) 549–569.