# Rxn Hypergraph: a Hypergraph Attention Model for Chemical Reaction
Representation

Mohammadamin Tavakoli ∗, Alexander Shmakov ∗, Francesco Ceccarelli, Pierre
Baldi

School of Information and Computer Sciences,
University of California, Irvine
<EMAIL_ADDRESS>

###### Abstract

It is fundamental for science and technology to be able to predict chemical
reactions and their properties. To achieve such skills, it is important to
develop good representations of chemical reactions, or good deep learning
architectures that can learn such representations automatically from the data.
There is currently no universal and widely adopted method for robustly
representing chemical reactions. Most existing methods suffer from one or more
drawbacks, such as: (1) lacking universality; (2) lacking robustness; (3)
lacking interpretability; or (4) requiring excessive manual pre-processing.
Here we exploit graph-based representations of molecular structures to develop
and test a hypergraph attention neural network approach to solve at once the
reaction representation and property-prediction problems, alleviating the
aforementioned drawbacks. We evaluate this hypergraph representation in three
experiments using three independent data sets of chemical reactions. In all
experiments, the hypergraph-based approach matches or outperforms other
representations and their corresponding models of chemical reactions while
yielding interpretable multi-level representations.

Over the past few years, artificial intelligence has refashioned organic
chemistry. Numerous problems such as chemical reaction prediction, synthesis
route planning, drug design, etc., have been benefited from the advancement of
deep learning methods (Coley et al., 2017; Fooshee et al., 2018; Schwaller et
al., 2019, 2020; Yang et al., 2019; Tavakoli et al., 2021). For instance,
accurate reaction yield predictions would help chemists to choose synthesis
routes across high-yielding chemical reactions. As another example, estimating
reaction rates via deep learning methods could circumvent the time and expense
required to experimentally measure the reaction rate of reactions. Not to
mention that because the true solution-phase reaction rates are bounded by the
rate of molecular diffusion, it is impossible to measure the accurate rate
experimentally. Lastly, predicting the outcome of chemical reactions using
deep learning methods would automate and accelerate the demanding processes of
drug design and discovery. All these problems require accurate predictions of
chemical reactions’ properties using deep learning models, consequently,
necessitating optimizing every aspect of such deep learning models. One
important facet that can significantly affect training dynamics and model
performance is the representation of the input data.

In an attempt to find a suitable representation of chemical reactions for
input into deep learning models, several methods have been proposed but each
of them suffers from certain shortcomings. These shortcomings may be
summarized with the following properties:

  1. 1.

Lack of universality Several representations are derived based on predefined
pattern matching algorithms (Schneider et al., 2015; Rogers and Hahn, 2010).
Since there is no learning process involved in deriving these representations,
they cannot be automatically adjusted for different predictive tasks.

  2. 2.

Lack of robustness Atoms and molecules have no inherent ordering within a
chemical reaction. Therefore, a robust representation must be invariant to
permutations on atoms and molecules within a given reaction. Methods built
upon the text representation of chemical reactions are an example of models
with the lack of robustness (Schwaller et al., 2021a, b). Using such methods,
one can obtain different outcomes by only permuting atoms and molecules in the
text representation of one single reaction.

  3. 3.

Non-interpretability It is vital for chemists to understand the reasoning
behind a reaction-level prediction. For example, capturing the correlation
between the presence and absence of certain functional groups or the
interaction between specific electrophiles and nucleophiles would provide
useful insight for chemists (Kadish et al., 2021; Mood et al., 2020). Thus it
is important that a representation can provide means to interpret the final
predictions.

  4. 4.

Need for expensive computations Lastly, some other representations require
hand-crafted implementations and processes in order to be used as the input of
a predictive model. These hand-crafted processes usually include running
pattern and subgraph matching algorithms (Kayala et al., 2011) or performing
massive data augmentations (Bjerrum, 2017) which are extremely time-consuming.

In what follows, we review these methods and evaluate them based on the four
mentioned shortcomings. Then we propose an augmented graph representation of a
chemical reaction called rxn-hypergraph which is designed to fix these
shortcomings and improve the proposed methods. Finally, through a set of
experiments on the classification and plausibility-based ranking of chemical
reactions, we empirically show the viability of our rxn-hypergraph
representation.

## Related Work

One commonly used method for numerically representing chemical reactions was
introduced in (Schneider et al., 2015). This representation is called
reactionFP and it is derived from the fingerprints of the molecules involved
in the chemical reaction. The reactionFP can be described as follows:

$w_{1}(\sum_{P_{i}}FP(P_{i})-\sum_{R_{i}}FP(R_{i}))+w_{2}(\sum_{A_{i}}FP(A_{i}))$
(1)

Where $A$, $R$, and $P$ represent the molecular agents, reactants, and
products respectively. $w_{1}$ and $w_{2}$ are two, potentially learnable,
parameters that adjust the contribution between the agent molecules and the
reactive molecules within the final representation. $FP(.)$ is a function that
outputs a traditional fingerprint of a given molecule such as ECFP4 (Rogers
and Hahn, 2010). Since extracting the traditional molecular fingerprints is
based on the presence or absence of predefined patterns and involves no
learning process, this representation lacks universality and cannot be
adjusted for a variety of reaction-level predictive tasks.

However, reactionFP was an early robust representation since Equation 1 is
trivially invariant to the permutation of the molecules involved in the
reaction. Additionally, many traditional molecular fingerprints such as ECFP4
(Rogers and Hahn, 2010) or AtomPair (Carhart et al., 1985) are also invariant
to atomic order within molecules, presenting a robust representation.
Nevertheless, since the molecular fingerprints are obtained using non-
invertible hashing mechanisms (Glen et al., 2006), this method cannot be
easily interpreted to discover high-level patterns after learning.
Additionally, according to Equation 1, obtaining reactionFP requires extra
hand-crafted computations including: (1) an accurate atom mapping between
reactants and products to identify the role of molecules in a chemical
reaction; and (2) extracting the vectorized form of traditional fingerprints.

Graph-based methods have also been used to extract more informative reaction
representations. To predict the outcome of chemical reactions, Fooshee et al.
(2018) used a neural network to rank a set of potential mechanistic reactions
based on their thermodynamic plausibility. They form two separate count bit
vector representations, one for the reactant and one for the product molecules
by recursively counting the predefined paths and trees of different sizes
rooted at each atom. Then, through a mutual information feature selection
stage, they extract the most informative set of these count bits for a given
downstream task. Finally, the chemical reaction may be represented as the
difference between the count vector of reactants and the count vector of
products. Since this method is based on a set of predefined patterns, it lacks
universality in the sense that it cannot capture necessary information for
different tasks. Although it is not discussed in Fooshee et al. (2018), the
mutual information feature selection step can potentially provide an
additional interpretable view of the final prediction. However, this requires
further expensive computations, not to mention that already massive
computation required for obtaining the count bit vectors.

Finally, the most successful reaction representation makes use of SMIRKS
(Weininger et al., 1989) of chemical reactions (Schwaller et al., 2019, 2021c,
2021a). SMIRKS is a well-defined domain-specific language with special
characters and grammatical rules for describing chemical transformations in
the SMILES strings (Weininger, 1988). Authors in Schwaller et al. (2021a)
deployed commonly used methods from natural language processing (NLP) to
encode the SMIRKS of chemical reactions into a continuous vector. They train
bi-directional encoder representation from transformers (BERT) (Devlin et al.,
2018) models to obtain task-specific reaction representations. They also
trained large sequence to sequence transformer models for masked language
model prediction (MLM) on the reaction SMIRKS. These models may also be used
to extract pre-trained representations of any given chemical reactions
(Pesciullesi et al., 2020). These transformer-based methods are highly
accurate across multiple reaction level predictive tasks (Schwaller et al.,
2021b, 2020). Such models provide a universal reaction representation that can
be used for numerous reaction-level predictive tasks. Additionally,
transformer architectures provide a character-level interpretable framework
over the SMILES strings. However, one major downside of these NLP models is
that their input, SMILES and SMIRKS strings, are not permutation invariant
with respect to the order of atoms and molecules. Depending on the canonical
SMILES parsing algorithm, atom labeling, and a few other details, a single
reaction may be correctly represented with many different SMIRKS. Some
permutation invariance may be recovered through massive data augmentation,
where all possible representations of the input reactions are generated and
randomly selected during training (Bjerrum, 2017). It has been shown the
performance of transformer models is highly dependent on the data
augmentation, where they show surprisingly poor performance without the data
augmentation (Bjerrum, 2017; Schwaller et al., 2019).

## Methods

In this section, we describe the reasons for constructing the chemical
reaction hypergraph (rxn-hypergraph) and why this hypergraph would yield an
abstract and powerful representation of a chemical reaction. Then we explain
the process of constructing the (rxn-hypergraph) for a given reactions, and
finally, we discuss how to train neural networks using this hypergraph
representation of chemical reactions.

### Why Rxn-hypergraph?

Graph neural networks and their variants have become the preeminent tool for
learning patterns and relations from graph-structured data (Duvenaud et al.,
2015; Schlichtkrull et al., 2018a; Kipf and Welling, 2016; Tavakoli and Baldi,
2020; Wieder et al., 2020; De Cao and Kipf, 2018). The main operation of graph
neural networks is to recursively update the representation of nodes using a
message-passing scheme only between the nodes and their neighbors. Then a
read-out function can be applied to the set of nodes’ representations to
provide an abstract representation of the entire graph.

There are several essential reaction-level properties where predicting them
would be highly beneficial to the entire field of chemoinformatics. However,
applying a graph neural network to the graph structure of chemical reactions
in its most raw form would not lead to a rich representation that captures
different properties of the reaction. Particularly, reactions are a more
general form of a graph that consists of multiple disconnected graph
components (molecules). The absence of a message-passing route between these
components would result in node representations that are independent of the
nodes and connections within the other graph components that are involved in
the reaction. Consequently, applying any form of a read-out function to these
independent node representations would result in a non-informative
representation of the entire reaction.

On the contrary, self attentional models (e.g. transformers) have been
impressively successful in reaction-level property predictions (Schwaller et
al., 2021b). The input to these models is the SMIRKS representation of a
chemical reaction in which atoms are represented by alphabetical characters
(tokens). Although this form of input representation is not invariant to the
permutation, it provides a suitable structure for the transformer
architectures. The key reason behind the success of such models can be found
in applying multiple layers of self-attention mechanisms to every pair of
input tokens. By this means, the representation of each token will be updated
by attending to all other tokens including the atoms within other molecules.

Inspired by this crucial factor in representing a reaction, we form a
hypergraph structure of a reaction by constructing efficient message passing
routes between every pair of atoms. These routes would improve the
representation of the reaction by: (1) enabling message-passing schemes
between every possible pair of atoms so the atom representations are updated
with respect to all other atoms within the reaction, and (2) providing a
learnable read-out function (i.e. pooling mechanism) which can attend to
different parts in different levels of the reaction which are informative for
a specific predictive task.

### Constructing Rxn-hypergraph

A chemical reaction with $N$ reactants and $M$ product molecules is described
by two distinct sets of disconnected graph components $R$ and $P$.
$R=\\{G_{i}^{r}\\}_{i=1}^{N}$ represents the set of reactant molecules, and
$P=\\{G_{i}^{p}\\}_{i=1}^{M}$ represents the set of product molecules.
Molecule $G_{i}$ with $n$ atoms (regardless of begin a reactant or product
molecule) is a graph $G_{i}=(V_{i},E_{i},A,S)$, where
$V_{i}=\\{a_{j}^{i}\\}_{j=1}^{n}$ is the set of nodes (atoms) and
$E_{i}=\\{(a_{u}^{i},a_{v}^{i})\\}$ is the set of edges (bonds). The set of
possible labels for the vertices in $A$ correspond to atom types (e.g. C, O),
and the set of possible labels for the edges in $S$ correspond to edge types
(single, double, triple, and aromatic). The idea behind forming the rxn-
hypergraph of a chemical reaction is to efficiently construct new message
passing routes between these disconnected graphs components (molecules) and
form one connected hypergraph to represent the entire reaction.

To form this hypergraph, we begin by unifying all the $G_{i}$s into one graph
$G=(V,E,A,S)$ where $V=\bigcup V_{i}^{r}+\bigcup V_{i}^{p}$ and $E=\bigcup
E_{i}^{r}+\bigcup E_{i}^{p}$, while $A$ and $S$ remains the same. For each of
the disconnected graph components $G_{i}$ (molecules), we add a hypernode to
the graph as a mol-hypernode $m_{i}$. Then we add two types of new edges to
the $G$: (1) a set of bidirectional edges connecting every atom to the mol-
hypernode of their parent molecule, mol-
atom$=\\{(m_{i},a_{j}^{i}),(a_{j}^{i},m_{i})\\}$, and (2) a set of
bidirectional edges connecting every pair of mol-hypernodes on either side of
the reaction, mol-mol$=\\{(m_{i},m_{j}),(m_{j},m_{i})\\}$. The edges of type
mol-mol would form two fully connected subgraphs between the mol-hypernodes on
each side of the reaction. We further augment $G$ by adding two more
hypernodes as rxn-hypernodes $x^{r}$ and $x^{p}$, one for the reactant and one
for the product side of the reaction. Then we add a new type of edge to the
graph: a set of unidirectional edges from each mol-hypernode to the rxn-
hypernode of the same side of the reaction, mol-
rxn$=\\{(m_{i}^{r},x^{r}),(m_{j}^{p},x^{p})\\}$. This augmented version of
graph $G$ is what we refer to as the rxn-hypergraph. Figure 1 shows a
schematic drawing of the rxn-hypergraph. It can also be represented it as
follows:

$\begin{split}&\textit{rxn-hypergraph}=(V^{*},E^{*},A^{*},S^{*})\\\
&\text{where:}\\\
&V^{*}=V\cup\\{m_{i}^{r}\\}_{i=1}^{N}\cup\\{m_{i}^{p}\\}_{i=1}^{M}\cup\\{x^{r},x^{p}\\},\\\
&E^{*}=E\cup\\{(m_{i},a_{j}^{i}),(a_{j}^{i},m_{i})\\}\cup\\{(m_{i}^{r},m_{j}^{r}),(m_{j}^{r},m_{i}^{r})\\}\\\
&\cup\\{(m_{i}^{p},m_{j}^{p}),(m_{j}^{p},m_{i}^{p})\\}\cup\\{(m_{i}^{r},x^{r}),(m_{j}^{p},x^{p})\\},\\\
&A^{*}=A\cup\\{\textit{mol-hypernode},\textit{rxn-hypernode}\\},\\\
&S^{*}=S\cup\\{\textit{atom-mol},\textit{mol-atom},\textit{mol-
mol},\textit{mol-rxn}\\}.\end{split}$ (2)

Figure 1: A depiction of the rxn-hypergraph corresponding to the reaction at
the bottom.

### Relational Graph Attention

Now that the rxn-hypergraph is formed, there is an intermolecular path of
length three ($a_{u}^{i}\rightarrow m_{i}\rightarrow a_{v}^{j}\rightarrow
m_{j}$) between every pair of atoms. Thus, we may construct robust contextual
node representations by applying more than three layers of graph convolution
neural networks to the rxn-hypergraph since this would ensure the receptive
field of every atom would include all other atoms. Also, the presence of the
unidirectional paths between mol-hypernodes and a rxn-hypernode, provides a
form pooling mechanism where the rxn-hypernode can fully represent one side of
the reaction. Assuming that both rxn-hypernodes carry an abstract
representation of the reactants and products molecules, merging them into one
vector would represent the entire chemical reaction.

To obtain this representation, we model the entire rxn-hypergraph as a
relational graph where the relations are represented by $S^{*}$ in the
previous section. We apply and compare two forms of graph neural networks:
relational graph convolution (RGCN) (Schlichtkrull et al., 2018b) and
relational graph attention (RGAT) (Veličković et al., 2018; Brody et al.,
2021). We summarize the layer-wise operations for each of these graph neural
networks.

First, a layer-wise operation of the relational graph convolution can be
described as follows:

$h_{i}^{l+1}=W_{o}h_{i}^{l}+\frac{1}{|\mathcal{N}_{s}(i)|}(\sum_{s\in
S^{*}}\sum_{j\in\mathcal{N}_{s}(i)}W_{n}h_{j}^{l})$ (3)

Where $h_{i}^{l}$ is the vector representation of atom $a_{i}$ and layer $l$.
$W_{o}$ and $W_{n}$ are two learnable weight matrices, and
$\mathcal{N}_{s}(i)$ represent the set of nodes adjacent to $a_{i}$ through
edge type $s$.

The layer-wise operation of our version of relational graph attention and how
the attention scores are computed are as follows. The term $\alpha_{ij}^{s}$
is the amount of attention that node $a_{i}$ would pay to its neighbor node
$a_{j}$ (through edge type $s$).

$h_{i}^{l+1}=\alpha_{ii}Wh_{i}^{l}+(\sum_{s\in
S^{*}}\sum_{j\in\mathcal{N}(i)}\alpha_{ij}Wh_{j}^{l})$ (4)

$\alpha_{(ij)}^{s}=\frac{\exp(A(Wh_{i}||Wh_{j}))}{\sum_{k\in\mathcal{N}_{s}(i)}\exp(A(Wh_{i}||Wh_{k}))}$
(5)

To train a reaction-level predictive model, we apply $L\geq 3$ layers of
relational graph attention/convolution to the rxn-hypergraph. Then final
representation of the reaction would be $X=f(x^{r^{L}},x^{p^{L}})$ where $f$
is a selected summarizing function to combine the two sides of the reaction.
We use subtraction ($f(x,y)=x-y$) and concatenation ($f(x,y)=x||y$) as the
summarizing function. Taking $X$ as the final latent representation of our
reaction, we apply a classification or regression head (typically a feed
forward MLP) to perform reaction-level classification or regression. In the
Experiments section, we evaluate the viability of our proposed method by
performing several reaction-level prediction tasks.

### Interpretability

Applying any form of graph neural networks on the rxn-hypergraph would provide
an interpretable framework using the standard attribution methods such as
Integrated gradient (Sundararajan et al., 2017) and Class Activation Map (CAM)
(Zhou et al., 2016). However, the structure of the rxn-hypergraph is capable
of providing a more in-depth interpretation. The attention weights obtained
from applying GAT can be used as a measure of importance for propagating
information relevant to the predictive task. Given that, we define three types
of interpretations for the predictions of a GAT model trained using rxn-
hypeergraph. (1) The atom-rxn interpretation: the multiplication of the
attention scores ($\alpha$) of the edges in path from an atom to the
corresponding rxn-hypernode (e.g. $a_{u}^{i}\rightarrow m_{i}^{r}\rightarrow
x^{r}$). These scores can be considered as the contribution of each atom to
the final representation of the reaction. (2) The node-node interpretation:
the average of the attention scores of the edges over all the GAT layers used
in architecture. These scores can measure the importance of bonds between
adjacent atoms in the final representation of the reaction. Also the relative
importance of the reacting molecules in the final representation of the
reaction can be measured by node-node scores of the edges between the mol-
hyernodes and their corresponding rxn-hypernode (e.g. $m_{i}^{r}\rightarrow
x^{r}$). This score can indicate the difference between reactants and
reagents. (3) The intermolecular atom-atom interpretation: the multiplication
of the attention scores of the edges in the path between pair of atoms of
different molecules (e.g. $a_{u}^{i}\rightarrow m_{i}^{r}\rightarrow
m_{j}^{r}\rightarrow a_{v}^{j}$). These scores can be considered as a measure
of the pair-wise correlation between atoms.

## Experiments

To evaluate the viability of the hypergraph representation in various
different circumstances, we perform three experiments on different reaction-
level predictive tasks using multiple datasets of chemical reactions.

### Data

For the first experiment, we use the dataset of chemical reactions from the US
patents office (USPTO)(Lowe, 2012) to train a model to classify these
reactions into the top 50 highly populated classes of chemical reactions. All
of these reactions were classified according to reaction classes presented in
Carey et al. (2006); Roughley and Jordan (2011) and the RSC’s RXNO ontology
(Kraut et al., 2013) using the NameRxn tool. Similar to the classification
experiment presented in Schneider et al. (2015), we randomly sample 1000
reactions for each of the 50 reaction classes. Then we split each class into a
subclass of 200 random reactions for training and 800 for testing. This
results in 10,000 training reactions and 40,000 test reactions which are
uniformly distributed over 50 classes of chemical reactions.

For the second experiments, we use an in-house curated dataset of reaction
mechanisms (reaction with a single transition state). This dataset consists of
three classes of mechanistic reactions: (1) over 11,000 polar reactions. The
reactions wherein a pair of electrons would transfer from an electron donor
orbital to an electron acceptor orbital; (2) 2800 of radical reactions. The
reactions that involve a radical species; and (3) 2600 pericyclic reactions.
The reactions wherein the transition state has a cyclic geometry. We split
this dataset into a 80 percent training dataset and 10 percent validation and
10 percent holdout for final testing.

Lastly, we redo the reaction-level experiment described in Fooshee et al.
(2018) as the reaction ranker stage. In this stage, the authors in Fooshee et
al. (2018) used a set of over 11,000 productive polar mechanistic reactions to
train a ranker system which ranks a reaction mechanism based on their
thermodynamic plausibility.

### Reaction Representations

For each experiment, we compare the performance of different models that are
trained using the corresponding reaction representation. These representations
which are introduced in the Related Work section are: (1) reactionFP
representation based on AP, Morgan2, and TT fingerprints; (2) the
representations from the transformer models on reaction SMIRKS, rxnfp (both
pre-trained and trained from scratch); and (3) representations from the
relational graph convolution/attention using rxn-hypergraph.

### Training and Hyperparameter Optimization

We train the graph attention/convolution layers using either cross-entropy or
mean squared error, depending on the task. We use the ADAM optimizer (Kingma
and Ba, 2017) and anneal the learning rate with an exponential schedule across
the training duration. Training is performed across 4 GPUs for a period of 500
epochs for each of the experiments explained below.

Additionally, we used SHERPA (Hertel et al., 2020) to optimize the
hyperparameters associated with each predictive task, guided by Bayesian
Optimization for each parameter. Specifically, we optimized the number of
graph layers $L$, size of the latent representation of nodes $D$, learning
rate, learning rate decay, and $L2$ weight regularization term. The final
parameters for each experiment are presented in Table 4.

### Classification of USPTO reactions

The results of this classification experiment are reported in Table 1. Both
transformer representation (rxnfp) and RGAT on rxn-hypergraph achieve the
highest accuracies. It is important to mention that this classification scheme
is highly dependent on the presence of certain molecules and compounds on the
reactant side of the reaction. Such textual dependencies are not representing
the underlying chemistry of the reaction which implies that a highly accurate
model that uses the text representations might not learn the actual chemistry
of the reactions and take advantage of the presence of absence of these
textual signatures.

Table 1: Comparing testing accuracy of different representation on a test set of 40,000 chemical reactions from USPTO dataset of chemical reactions. Representation | Network | Accuracy
---|---|---
reactionFP | AP | 0.854
Morgan2 | 0.850
TT | 0.852
Transformers | pre-trained rxnfp | 0.862
rxnfp | 0.925
rxn-hypergraph | RGCN | 0.909
RGAT | 0.928

### Classification of Mechanistic Reactions

Here we classify the reactions into three classes of polar, radical, and
pericyclic reactions. Since the classification scheme is only based on the
pair of reacting orbitals (i.e. single transition states), the underlying
chemistry might not be complicated for models to learn. Nevertheless, the
graph attentional models on rxn-hypergraph are outperforming other models and
representations.

For this particular experiment, we also show the interpretations of the final
predictions using the three metrics described in the Interpretability section.
Figures 2, 3, and 4 are illustrating the interpretation results. In each
figure the reaction with labeled atoms is depicted at the top. The atom-rxn
scores are shown in the middle while the node-node and intermolecular atom-
atom scores are shown at the bottom right and bottom left of the figures.

Figure 2: Interpretability plots for a radical reaction for the task of
classifying the mechanistic reactions. Figure 3: Interpretability plots for a
pericyclic reaction for the task of classifying the mechanistic reactions.
Figure 4: Interpretability plots for a polar reaction for the task of
classifying the mechanistic reactions.

### Plausibility Ranking of Polar Mechanistic Reactions

This experiment was first introduced in Fooshee et al. (2018), where they rank
a set of possible reactions from the interactions between one set of reactant
molecules. We precisely follow the procedure of ranking described in Fooshee
et al. (2018), we use Siamese network (Bromley et al., 1993; Chicco, 2021) to
train ranker models for reactionFP and transformer models. In each branch of
the Siamese network, we used the hyperparameters used in the previous section.

#### Ranking Network Architecture

We learn pairwise rankings between different mechanisms using the DirectRanker
architecture introduced in Köppel et al. (2020). We first compute the learned
latent representations of a pair of reactions. Afterward, we use the
DirectRanker method of subtracting the two latent representations before
feeding them through a bias-free fully-connected layer and a sign-preserving
anti-symmetric nonlinearity which produces a scalar value for each pair of
events. We train this network using mean squared error to predict which of the
two input reactions is more plausible than the other, targets of $1$ or $-1$.

After training this pairwise ranker, we still need a method for determining a
final ranked order on a list of plausible mechanisms. To accomplish this, we
use the ranked-pair voting method (Tideman, 1987). This algorithm allows us to
convert a pair-wise ranking matrix between all viable reaction mechanisms to
an ordered list by preference. Ranked pairs construct a directed acyclic graph
based on the sorted pair-wise score between different elements. The acyclic
property is maintained by ignoring any pairs which would introduce a cycle.
After all of the pairs are exhausted, we produce the final ordering by
following a topological ordering on the nodes of the graph. The resulting
ordering is guaranteed to obey certain criterion which is useful for this task
such as independence to irrelevant alternatives, which means that extra
implausible reactions which the network is unsure of will not spoil the top
rankings.

Table 2: Comparing testing accuracy of different representation on a test set of 5,455 chemical reactions from an in house dataset of mechanistic reactions. Representation | Network | Accuracy
---|---|---
reactionFP | AP | 0.915
Morgan2 | 0.915
TT | 0.913
Transformers | pre-trained rxnfp | 0.955
rxnfp | 0.974
rxn-hypergraph | RGCN | 0.988
RGAT | 0.990
Table 3: The top1, 2, 5, and 10 prediction accuracy of plausibility ranking models with different representations of chemical reaction. All the metrics are computed for a test set of 200 real world mechanistic reactions. Representation |  | top1 | top2 | top5 | top10
---|---|---|---|---|---
reactionFP | AP | 58.01 | 67.05 | 77.60 | 84.33
Morgan2 | 59.03 | 69.14 | 78.24 | 85.01
TT | 58.41 | 68.22 | 77.31 | 84.14
Transformers | pre-trained rxnfp | 81.31 | 84.19 | 89.60 | 92.33
rxnfp | 89.14 | 93.22 | 96.09 | 98.55
rxn-hypergraph | RGCN | 82.67 | 92.57 | 97.03 | 99.01
RGAT | 84.23 | 96.06 | 98.57 | 99.28

Out of three experiment, ranking reaction based in thermodynamic plausibility
requires a deeper understanding of the underlying chemistry. In this task, the
correlation between the textual signatures and plausibility of a reaction is
minimal. This potentially explains the results presented in Table 3 where the
RGAT model on rxn-hypergraph is outperforming other models especially the
models based on the text representations with a considerable margin.

Task | USPTO Classification | Mechanism Classification | Polar Plausibility
---|---|---|---
Number of Graph Attention Layers | 10 | 10 | 5
Latent Representation Dimension | 128 | 128 | 64
Learning Rate | 4.11 $\times 10^{-4}$ | 4.11 $\times 10^{-4}$ | 1.14 $\times 10^{-2}$
Learning Rate Decay | 0.999995 | 0.999995 | 0.99986
$L2$ Weight Regularization | 9.75 $\times 10^{-5}$ | 9.75 $\times 10^{-5}$ | 7.28 $\times 10^{-5}$
Table 4: Table of hyperparameters selected through Bayesian Optimization using
SHERPA (Hertel et al., 2020).

## Discussion and Conclusion

We proposed rxn-hypergraph representation of a chemical reaction which is
suitable for training graph neural networks for reaction-level predictive
task. The key idea behind forming the rxn-hypergraph is to construct efficient
message passing routes which provides a platform for (1) updating atom
representation based on the atom and molecules if the other reacting
molecules, and (2) a global pooling mechanism. Rxn-hypergraph is designed to
be a universal and permutation invariant representation that adapt to any
downstream predictive task. There is no manual and hand-crafted pre-processing
stages involved in computing the rxn-hypergraph and it provides different
levels of interpretability. There are two potential aspects of this work which
are left to the future work: (1) several other demanding and useful reaction-
level predictive tasks such as yield prediction, and reaction rate constant
prediction that can be benefited by rxn-hypergraph, and (2) more complicated
and expressive attention mechanisms such as multiplicative attentions and
transformers can be applied to rxn-hypergraph which might result in more
powerful reaction representations.

## References

  * Coley et al. [2017] Connor W Coley, Regina Barzilay, Tommi S Jaakkola, William H Green, and Klavs F Jensen.  Prediction of organic reaction outcomes using machine learning.  _ACS central science_ , 3(5):434–443, 2017.
  * Fooshee et al. [2018] David Fooshee, Aaron Mood, Eugene Gutman, Mohammadamin Tavakoli, Gregor Urban, Frances Liu, Nancy Huynh, David Van Vranken, and Pierre Baldi.  Deep learning for chemical reaction prediction.  _Molecular Systems Design & Engineering_, 3(3):442–452, 2018.
  * Schwaller et al. [2019] Philippe Schwaller, Teodoro Laino, Théophile Gaudin, Peter Bolgar, Christopher A Hunter, Costas Bekas, and Alpha A Lee.  Molecular transformer: a model for uncertainty-calibrated chemical reaction prediction.  _ACS central science_ , 5(9):1572–1583, 2019\.
  * Schwaller et al. [2020] Philippe Schwaller, Riccardo Petraglia, Valerio Zullo, Vishnu H Nair, Rico Andreas Haeuselmann, Riccardo Pisoni, Costas Bekas, Anna Iuliano, and Teodoro Laino.  Predicting retrosynthetic pathways using transformer-based models and a hyper-graph exploration strategy.  _Chemical science_ , 11(12):3316–3325, 2020.
  * Yang et al. [2019] Kevin Yang, Kyle Swanson, Wengong Jin, Connor Coley, Philipp Eiden, Hua Gao, Angel Guzman-Perez, Timothy Hopper, Brian Kelley, Miriam Mathea, et al.  Analyzing learned molecular representations for property prediction.  _Journal of chemical information and modeling_ , 59(8):3370–3388, 2019.
  * Tavakoli et al. [2021] Mohammadamin Tavakoli, Aaron Mood, David Van Vranken, and Pierre Baldi.  Quantum mechanics and machine learning synergies: Graph attention neural networks to predict chemical reactivity.  _arXiv preprint arXiv:2103.14536_ , 2021.
  * Schneider et al. [2015] Nadine Schneider, Daniel M Lowe, Roger A Sayle, and Gregory A Landrum.  Development of a novel fingerprint for chemical reactions and its application to large-scale reaction classification and similarity.  _Journal of chemical information and modeling_ , 55(1):39–53, 2015.
  * Rogers and Hahn [2010] David Rogers and Mathew Hahn.  Extended-connectivity fingerprints.  _Journal of chemical information and modeling_ , 50(5):742–754, 2010.
  * Schwaller et al. [2021a] Philippe Schwaller, Daniel Probst, Alain C Vaucher, Vishnu H Nair, David Kreutter, Teodoro Laino, and Jean-Louis Reymond.  Mapping the space of chemical reactions using attention-based neural networks.  _Nature Machine Intelligence_ , 3(2):144–152, 2021a.
  * Schwaller et al. [2021b] Philippe Schwaller, Alain C Vaucher, Teodoro Laino, and Jean-Louis Reymond.  Prediction of chemical reaction yields using deep learning.  _Machine Learning: Science and Technology_ , 2(1):015016, 2021b.
  * Kadish et al. [2021] Dora Kadish, Aaron D Mood, Mohammadamin Tavakoli, Eugene S Gutman, Pierre Baldi, and David L Van Vranken.  Methyl cation affinities of canonical organic functional groups.  _The Journal of Organic Chemistry_ , 86(5):3721–3729, 2021.
  * Mood et al. [2020] Aaron Mood, Mohammadamin Tavakoli, Eugene Gutman, Dora Kadish, Pierre Baldi, and David L Van Vranken.  Methyl anion affinities of the canonical organic functional groups.  _The Journal of organic chemistry_ , 85(6):4096–4102, 2020.
  * Kayala et al. [2011] Matthew A Kayala, Chloé-Agathe Azencott, Jonathan H Chen, and Pierre Baldi.  Learning to predict chemical reactions.  _Journal of chemical information and modeling_ , 51(9):2209–2222, 2011.
  * Bjerrum [2017] Esben Jannik Bjerrum.  Smiles enumeration as data augmentation for neural network modeling of molecules.  _arXiv preprint arXiv:1703.07076_ , 2017.
  * Carhart et al. [1985] Raymond E Carhart, Dennis H Smith, and R Venkataraghavan.  Atom pairs as molecular features in structure-activity studies: definition and applications.  _Journal of Chemical Information and Computer Sciences_ , 25(2):64–73, 1985.
  * Glen et al. [2006] Robert C Glen, Andreas Bender, Catrin H Arnby, Lars Carlsson, Scott Boyer, and James Smith.  Circular fingerprints: flexible molecular descriptors with applications from physical chemistry to adme.  _IDrugs_ , 9(3):199, 2006.
  * Weininger et al. [1989] David Weininger, Arthur Weininger, and Joseph L Weininger.  Smiles. 2. algorithm for generation of unique smiles notation.  _Journal of chemical information and computer sciences_ , 29(2):97–101, 1989.
  * Schwaller et al. [2021c] Philippe Schwaller, Benjamin Hoover, Jean-Louis Reymond, Hendrik Strobelt, and Teodoro Laino.  Extraction of organic chemistry grammar from unsupervised learning of chemical reactions.  _Science Advances_ , 7(15):eabe4166, 2021c.
  * Weininger [1988] David Weininger.  Smiles, a chemical language and information system. 1. introduction to methodology and encoding rules.  _Journal of chemical information and computer sciences_ , 28(1):31–36, 1988.
  * Devlin et al. [2018] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_ , 2018.
  * Pesciullesi et al. [2020] Giorgio Pesciullesi, Philippe Schwaller, Teodoro Laino, and Jean-Louis Reymond.  Transfer learning enables the molecular transformer to predict regio-and stereoselective reactions on carbohydrates.  _Nature communications_ , 11(1):1–8, 2020.
  * Duvenaud et al. [2015] David Duvenaud, Dougal Maclaurin, Jorge Aguilera-Iparraguirre, Rafael Gómez-Bombarelli, Timothy Hirzel, Alán Aspuru-Guzik, and Ryan P Adams.  Convolutional networks on graphs for learning molecular fingerprints.  _arXiv preprint arXiv:1509.09292_ , 2015.
  * Schlichtkrull et al. [2018a] Michael Schlichtkrull, Thomas N Kipf, Peter Bloem, Rianne Van Den Berg, Ivan Titov, and Max Welling.  Modeling relational data with graph convolutional networks.  In _European semantic web conference_ , pages 593–607. Springer, 2018a.
  * Kipf and Welling [2016] Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ , 2016.
  * Tavakoli and Baldi [2020] Mohammadamin Tavakoli and Pierre Baldi.  Continuous representation of molecules using graph variational autoencoder.  _arXiv preprint arXiv:2004.08152_ , 2020.
  * Wieder et al. [2020] Oliver Wieder, Stefan Kohlbacher, Mélaine Kuenemann, Arthur Garon, Pierre Ducrot, Thomas Seidel, and Thierry Langer.  A compact review of molecular property prediction with graph neural networks.  _Drug Discovery Today: Technologies_ , 2020.
  * De Cao and Kipf [2018] Nicola De Cao and Thomas Kipf.  Molgan: An implicit generative model for small molecular graphs.  _arXiv preprint arXiv:1805.11973_ , 2018.
  * Schlichtkrull et al. [2018b] Michael Schlichtkrull, Thomas N. Kipf, Peter Bloem, Rianne van den Berg, Ivan Titov, and Max Welling.  Modeling relational data with graph convolutional networks.  In Aldo Gangemi, Roberto Navigli, Maria-Esther Vidal, Pascal Hitzler, Raphaël Troncy, Laura Hollink, Anna Tordai, and Mehwish Alam, editors, _The Semantic Web_ , pages 593–607, Cham, 2018b. Springer International Publishing.  ISBN 978-3-319-93417-4.
  * Veličković et al. [2018] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio.  Graph attention networks, 2018.
  * Brody et al. [2021] Shaked Brody, Uri Alon, and Eran Yahav.  How attentive are graph attention networks?, 2021.
  * Sundararajan et al. [2017] Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  In _International Conference on Machine Learning_ , pages 3319–3328. PMLR, 2017.
  * Zhou et al. [2016] Bolei Zhou, Aditya Khosla, Agata Lapedriza, Aude Oliva, and Antonio Torralba.  Learning deep features for discriminative localization.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 2921–2929, 2016.
  * Lowe [2012] Daniel Mark Lowe.  _Extraction of chemical structures and reactions from the literature_.  PhD thesis, University of Cambridge, 2012.
  * Carey et al. [2006] John S Carey, David Laffan, Colin Thomson, and Mike T Williams.  Analysis of the reactions used for the preparation of drug candidate molecules.  _Organic & biomolecular chemistry_, 4(12):2337–2347, 2006.
  * Roughley and Jordan [2011] Stephen D Roughley and Allan M Jordan.  The medicinal chemist’s toolbox: an analysis of reactions used in the pursuit of drug candidates.  _Journal of medicinal chemistry_ , 54(10):3451–3479, 2011.
  * Kraut et al. [2013] Hans Kraut, Josef Eiblmaier, Guenter Grethe, Peter low, Heinz Matuszczyk, and Heinz Saller.  Algorithm for reaction classification.  _Journal of chemical information and modeling_ , 53(11):2884–2895, 2013.
  * Kingma and Ba [2017] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * Hertel et al. [2020] Lars Hertel, Julian Collado, Peter Sadowski, Jordan Ott, and Pierre Baldi.  Sherpa: Robust hyperparameter optimization for machine learning.  _SoftwareX_ , 12:100591, 2020.  ISSN 2352-7110.  doi: https://doi.org/10.1016/j.softx.2020.100591.  URL https://www.sciencedirect.com/science/article/pii/S2352711020303046.
  * Bromley et al. [1993] Jane Bromley, James W Bentz, Léon Bottou, Isabelle Guyon, Yann LeCun, Cliff Moore, Eduard Säckinger, and Roopak Shah.  Signature verification using a “siamese” time delay neural network.  _International Journal of Pattern Recognition and Artificial Intelligence_ , 7(04):669–688, 1993.
  * Chicco [2021] Davide Chicco.  Siamese neural networks: An overview.  _Artificial Neural Networks_ , pages 73–94, 2021.
  * Köppel et al. [2020] Marius Köppel, Alexander Segner, Martin Wagener, Lukas Pensel, Andreas Karwath, and Stefan Kramer.  Pairwise learning to rank by neural networks revisited: Reconstruction, theoretical analysis and practical performance.  In Ulf Brefeld, Elisa Fromont, Andreas Hotho, Arno Knobbe, Marloes Maathuis, and Céline Robardet, editors, _Machine Learning and Knowledge Discovery in Databases_ , pages 237–252, Cham, 2020. Springer International Publishing.  ISBN 978-3-030-46133-1.
  * Tideman [1987] T. N. Tideman.  Independence of clones as a criterion for voting rules.  _Social Choice and Welfare_ , 4(3):185–206, 1987\.  ISSN 01761714, 1432217X.  URL http://www.jstor.org/stable/41105866.


# Late-time Cosmology of scalar field assisted $f(\mathcal{G})$ gravity

S.A. Venikoudis, K.V. Fasoulakos, F.P. Fronimos
venikoudis@gmail.comkvfasoulakos<EMAIL_ADDRESS>Department
of Physics, Aristotle University of Thessaloniki, Thessaloniki 54124, Greece

###### Abstract

In this work we present the late-time behaviour of the Universe in the context
of Einstein-Gauss-Bonnet gravitational theory. The theory involves a scalar
field, which represents low-effective quantum corrections, assisted by a
function $f(\mathcal{G})$ solely depending from the Gauss-Bonnet topological
invariant $\mathcal{G}$. It is considered that the dark energy serves as the
impact of all geometric terms, which are included in the gravitational action
and the density of dark energy acts as a time dependent cosmological constant
evolving with an infinitesimal rate and driving the Universe into an
accelerating expansion. We examine two cosmological models of interest. The
first involves a canonical scalar field in the presence of a scalar potential
while the second, involves a scalar field which belongs to a generalized class
of theories $f(\phi,X)$ namely the k-essence scalar field in the absence of
scalar potential. As it is proved, the aforementioned models are in
consistency with the latest Planck data and in relatively good agreement with
the $\Lambda$CDM standard cosmological model. The absence of dark energy
oscillations at the early stages of matter dominated era, which appear in
alternative scenarios of cosmological dynamics in the context of modified
$f(R)$ gravitational theories, indicates an advantage of the theory for the
interpretation of late-time phenomenology.

###### pacs:

04.50.Kd, 95.36.+x, 98.80.-k, 98.80.Cq,11.25.-w

##  I Introduction

Observations on Supernovae type Ia [1] have shown that the Universe expands
with an accelerating rate. However, the cosmological dynamics during the late-
time era of our Universe remains one of the most intriguing mystery for modern
theoretical cosmologists until today. The most successful interpretation for
the cosmological evolution is given by the $\Lambda$CDM model, which is the
most well-established theory in consistency with the observations from the
data of Cosmic Microwave Background. In the context of the aforementioned
model, $\Lambda$ serves as a cosmological constant, which dominates in late-
time era compared with the relativistic and non-relativistic matter. However,
several problems have yet to be solved with the most important of them, the
understanding of the nature of the cosmological constant, which is related
with the dark energy. In addition, based on the Einstein-Hilbert gravitational
theory the equation of state parameter of $\Lambda$ is considered that remains
constant, but there is lack of certainty from observational perspective.

Beyond the Einstein’s gravity many modified gravitational theories have been
arisen as extensions, in order to provide a solid explanation of the late-time
phenomenon, see Ref. [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
37, 38, 39, 40, 41, 42, 43, 44]. These theories involve higher order
gravitational terms which originate from curvature invariants. The main
advantage of these theories is that, they can provide unification models
between the post quantum primordial era namely, the inflationary era and the
late-time era compatible with the latest Planck data and with the $\Lambda$CDM
model. On the other hand, $f(R)$ theories have a serious drawback under
certain circumstances. As the numerical value of redshift increases, dark
energy oscillations seams to appear in several models, for instance power-law
models with an exponent $n$ satisfying the condition $0<n<1$, specifically
during the last stages of the matter domination [45]. These oscillations
originate from the existence of higher derivatives of the Hubble rate. If one
uses nontrivial models such as singular power-law models $f(R)$ with exponents
$n\leq-1$ oscillations seem to disappear, see Ref [46].

In this paper it is demonstrated the phenomenology of the late-time Universe
of a given scalar-tensor theory in the context of string-inspired gravity, see
Ref. [47, 48, 49, 50, 50, 51, 52, 53]. The theory involves a scalar field
which is minimally coupled with the Ricci scalar and a function
$f(\mathcal{G})$ solely depending from the Gauss-Bonnet topological invariant
$\mathcal{G}$, representing higher curvature corrections and assists the
accelerating expansion. For the sake of simplicity, no coupling between the
scalar field and the Gauss-Bonnet topological invariant is assumed in the
present article. Combining both a scalar field and a tensor Gauss-Bonnet model
has the advantage that a unified description between the early and late time
may be achieved. One plausible scenario could involve a dominant scalar field
during the inflationary era where one for instance can implement the usual
prescription of slow-roll for potential driven inflation while in the late era
the extra curvature corrections, if introduced appropriately, can drive the
accelerated expansion of the universe. Under this assumption, we shall examine
two models of interest, the former involves a scalar field with a canonical
kinetic term and scalar potential, while the latter involves a k-essence
scalar field with a combination of linear and higher order kinetic term in the
absence of scalar potential for simplicity. In both cases the function
$f(\mathcal{G})$ has non linear terms that effectively drive the late-time.
The motivation for the investigation of theories with k-essence scalar field
is based on the recent striking observations from GW170817 [54]. The validity
of these theories has been ascertained, since they can provide viable
phenomenological models in primordial and late-time era [55, 56].

An important feature of the present theory is that it can be compatible with
latest observations and replicate the $\Lambda$CDM model. While the latter is
showcased in the present article, a thorough study of tensor perturbations is
a quite challenging topic and is thus not available, in principle however
tensor perturbations should be quite similar to Ref. [57]. The important
aspect that is worth being highlighted is that the models at hand are free of
ghost instabilities given that the Ricci scalar appears in a linear form. More
involved theories with terms of the form $f(R,\mathcal{G})$ are known to
suffer from ghosts therefore to avoid this possibility, we shall limit our
work to only $R+f(\mathcal{G})$ models and focus mainly on the late-time
evolution of the universe.

Our goal of this work is to examine the phenomenology of late-time
cosmological era, considering that the dark energy serves as the impact of all
the involved geometric terms, in the gravitational action, something which is
highlighted subsequently. Specifically, we shall examine the effect of the
scalar field, which serves as low-energy quantum corrections assisted by the
function $f(\mathcal{G})$. Our analysis is based on Ref. [58, 59] where
similar procedures have been performed. First of all, it is convenient to
apply two variable changes in Friedmann and continuity equation of the scalar
field, in order to proceed with the numerical solution of these equations. The
first change involves the replacement of the cosmic time with redshift.
Specifically, we introduce a differential operator which connects the cosmic
time with redshift given that the Hubble’s parameters is a function of
redshift. As a result, the involved physical quantities are expressed in terms
of redshift which yields a more manageable set of equations that need to be
solved. Furthermore, we introduce an auxiliary statefinder function which is
connected with the density of the dark energy to male the Friedmann equation
dimensionless. By utilizing the solutions for the statefinder function and the
scalar field one can extract the current numerical values for several
statefinder functions in order to properly study the dynamics of the universe
and make comparisons with the $\Lambda$CDM, but also trace back their
evolution near the matter dominated era. In the present article, we shall
focus on several statefinder parameters for the sake of consistency, namely
the deceleration parameter, the jerk, the snap and parameter Om as well as the
equation of state of dark energy and the density parameter of dark energy
$\Omega_{DE}$. The last step in order to ascertain the validity of each model.
As showcased, the models can be made to replicate the effects of $\Lambda$CDM
while an infinitesimal evolution with respect to redshift is showcased.

This paper is organized as follows: In section II the theoretical framework of
late-time Cosmology of scalar-$f(G)$ gravity is presented thoroughly as
mentioned before for the first model of interest. In section III we discuss
the validity of a model defining the function $f(\mathcal{G})$ and the scalar
potential. In sections IV and V it is demonstrated briefly the cosmological
evolution in the presence of a k-essence scalar field and $f(\mathcal{G})$ in
the absence of scalar potential, as it is usual when the theory involves
higher order kinetic terms. Finally, the conclusions of our work follow at the
end of the paper.

##  II late-time cosmology of scalar - $f(\mathcal{G})$ gravity

We commence our study by properly specifying the gravitational action. The
first proposed model involves the case of $f(\mathcal{G})$ gravity accompanied
by a canonical scalar field, corresponding to the following action,

$\centering
S=\int{d^{4}x\sqrt{-g}\left(\frac{R}{2\kappa^{2}}-\frac{1}{2}\nabla_{\mu}\phi\nabla^{\mu}\phi-V(\phi)+\frac{f(\mathcal{G})}{2}+\mathcal{L}_{m}\right)}\,,\@add@centering$
(1)

where R denotes the Ricci scalar, $\mathcal{L}_{m}$ represents the Lagrangian
density of the perfect fluids and $f(\mathcal{G})$ is an arbitrary function
for the time being which depends solely on the Gauss-Bonnet topological
invariant $\mathcal{G}$. Furthermore, g is the determinant of the metric
tensor $g^{\mu\nu}$, $\kappa=\frac{1}{M_{P}}$ is the gravitational constant
while, $M_{P}$ denotes the reduced Planck mass. In addition, the terms
$\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi$ and $V(\phi)$
correspond to the kinetic term of the scalar field and the scalar potential
respectively. The Gauss-Bonnet topological invariant is defined as a
combination of higher order curvature terms namely,
$\mathcal{G}=R_{\mu\nu\sigma\rho}R^{\mu\nu\sigma\rho}-4R_{\mu\nu}R^{\mu\nu}+R^{2}$
with $R_{\mu\nu\sigma\rho}$, being the Riemann curvature tensor, while
$R_{\mu\nu}$ and $R$ represents the Ricci tensor and scalar respectively.
Greek indices run from $0$ to $3$ whereas Latin describes only the spatial
components. As mentioned before the $f(G)$ function can not has a linear form
as $\int{d^{4}x\sqrt{-g}\mathcal{G}}=0$ since it is a total derivative. The
case of keeping a linear term by means of non minimal coupling with the scalar
field, although interesting and quite powerful for both the inflationary and
the late era, it shall not be considered in the present article. Throughout
our analysis, we shall assume that the geometric background corresponds to
that of a flat Friedman-Robertson-Walker (FRW) metric, meaning that the line
element has the following form,

$\centering ds^{2}=-dt^{2}+a(t)^{2}\delta_{ij}dx^{i}dx^{j},\,\@add@centering$
(2)

where $a(t)$ denotes the scale factor of the Universe and the metric tensor
has the form of $g_{\mu\nu}=diag(-1,a(t)^{2},a(t)^{2},a(t)^{2})$. Due to the
fact that the metric is assumed to be flat, the curvature terms, being the
Ricci scalar and Gauss-Bonnet topological invariant, are written explicitly in
terms of Hubble’s parameter $H=\frac{\dot{a}}{a}$. Specifically,
$R=6(2H^{2}+\dot{H})$ and $\mathcal{G}=24H^{2}(H^{2}+\dot{H})$ where the dot
as usual implies differentiation with respect to cosmic time $t$. Furthermore,
we shall assume that the aforementioned scalar field is homogeneous so as to
simplify our study. This is a reasonable assumption which leads to the kinetic
term being written as $X=-\frac{1}{2}\dot{\phi}^{2}$. Therefore, by
implementing the variation principle in the gravitational action with respect
to the metric tensor $g^{\mu\nu}$ and the scalar field $\phi$, the field
equations for gravity and the continuity equation for the scalar field can be
extracted. The aforementioned equations are,

$\centering\frac{3H^{2}}{\kappa^{2}}=\frac{1}{2}\dot{\phi}^{2}+V(\phi)+\frac{\mathcal{G}f_{\mathcal{G}}-f}{2}-12\dot{f}_{\mathcal{G}}H^{3}+\rho_{m}\,,\@add@centering$
(3)
$\centering-\frac{2\dot{H}}{\kappa^{2}}=\dot{\phi}^{2}+4H^{2}(\ddot{f}_{\mathcal{G}}-H\dot{f}_{\mathcal{G}})+8\dot{f}_{\mathcal{G}}H\dot{H}+\rho_{m}+P_{m}\,,\@add@centering$
(4) $\centering\ddot{\phi}+3H\dot{\phi}+\frac{dV}{d\phi}=0,\ \@add@centering$
(5)

where $\rho_{m}$ and $P_{m}$ express the density and pressure corresponding on
matter. The matter term is considered as the total contribution from non-
relativistic matter namely, baryons, leptons and Cold Dark matter being the
most dominant contribution while relativistic matter contains photons and
neutrinos. As usual, matter can be described as a perfect fluid for any
cosmological era of interest. In late-time cosmological era, the perfect fluid
consists of Cold Dark Matter and radiation hence, the corresponding density is
given by the following equation,

$\centering\rho_{m}=\rho_{nr}^{(0)}\frac{1}{\alpha^{3}}\left(1+\frac{\chi}{\alpha}\right),\@add@centering$
(6)

where $\chi=\frac{\rho_{r}^{(0)}}{\rho_{nr}^{(0)}}$. The term $\rho_{r}^{(0)}$
represents the present density of relativistic matter while, $\rho_{nr}^{(0)}$
being the current density of non relativistic matter. Such parameters shall be
specified subsequently in the numerical studies. In addition, the pressure of
the perfect fluid is defined as,

$P=\sum_{i}\omega_{i}\rho_{i},$ (7)

where $\omega_{i}$ is the equation of state parameter for both relativistic
and non-relativistic matter. Here, barotropic matter is considered where
$\omega_{nr}=0$ and $\omega_{r}=\frac{1}{3}$. The continuity equation for
every component of matter has the following simple form,

$\dot{\rho}_{i}+3H\rho_{i}(1+\omega_{i})=0,$ (8)

where as shown, no mixing between matter components is observed. A convenient
way for the investigation of the dynamics of late-time Universe is the
replacement of the cosmic time with the redshift as a dynamical parameter.
Redshift is defined as following,

$\centering 1+z=\frac{1}{a(t)}\,,\@add@centering$ (9)

where the current scale factor is considered to be equal to unity for
simplicity. Hence, the present value of the redshift is equal to zero.
According to the previous variable change, the involved time derivatives can
be expressed in terms of the redshift introducing the following differential
operator,

$\centering\frac{d}{dt}=-H(z)(1+z)\frac{d}{dz}\,.\@add@centering$ (10)

Based on the aforementioned assumption, the gravitational terms in the
equations of motion and the time derivatives of the scalar field can be
rewritten as:

$\dot{H}=-H(1+z)H^{\prime},$ (11) $R=12H^{2}-6H(1+z)H^{\prime},$ (12)
$\mathcal{G}=24H^{2}(H^{2}-H(1+z)H^{\prime}),$ (13)
$\dot{\mathcal{G}}=24(1+z)^{2}H^{3}\left(3H^{\prime
2}+HH^{\prime\prime}-\frac{3HH^{\prime}}{1+z}\right),$ (14)
$\dot{\phi}=-H(1+z)\phi^{\prime},$ (15)
$\ddot{\phi}=H^{2}(1+z)^{2}\phi^{\prime\prime}+H^{2}(1+z)\phi^{\prime}+HH^{\prime}(1+z)^{2}\phi^{\prime},$
(16)

where the prime denotes the differentiation with respect to the redshift.
Furthermore, for the time being, in order to investigate the dynamics of the
cosmological evolution the Hubble rate and its derivatives are replaced by the
dimensionless statefinder $y_{H}$ which is defined as,

$\centering y_{H}(z)=\frac{\rho_{DE}}{\rho_{nr}^{(0)}}\,,\@add@centering$ (17)

where $\rho_{DE}$ is the dark energy density which involves all the geometric
terms appearing on the right hand side of Friedmann’s equation hence, it is
considered as,

$\centering\rho_{DE}=\frac{1}{2}\dot{\phi}^{2}+V(\phi)+\frac{\mathcal{G}f_{\mathcal{G}}-f}{2}-12\dot{f}_{\mathcal{G}}H^{3}\,,\@add@centering$
(18)

while according to Raychaudhuri equation, the pressure which originates from
the dark energy is:

$P_{DE}=\dot{\phi}^{2}+4H^{2}(\ddot{f}_{\mathcal{G}}-H\dot{f}_{\mathcal{G}})+8\dot{f}_{\mathcal{G}}H\dot{H}-\rho_{DE}.$
(19)

This designation of dark energy density and pressure is performed in order to
rewrite the Friedmann and Raychaudhuri equations in the usual form as
showcased below,

$\centering
3H^{2}=\kappa^{2}\left(\rho_{m}+\rho_{DE}\right)\,,\@add@centering$ (20)
$\centering-2\dot{H}=\kappa^{2}\left(\rho_{m}+P_{m}+\rho_{DE}+P_{DE}\right).\@add@centering$
(21)

Obviously, instead of a cosmological constant now a function of redshift
appears in both equations, therefore is rightfully deserves the name dark
energy as it can be a formal definition of it. Another interesting conclusion
which can easily be inferred from equations (20) and (21) is that dark energy
must behaves as a perfect fluid since Hubble and matter components are also
considered as perfect fluids. As a result the continuity equation of dark
energy has the usual form,

$\centering\dot{\rho}_{DE}+3H(\rho_{DE}+P_{DE})=0,\@add@centering$ (22)

where as was the case with the matter components, no coupling between dark
energy and matter is present. Note also that a constant dark energy density,
meaning the cosmological constant $\Lambda$, requires an equation of state
$\omega_{DE}=-1$. Therefore, one could extract similar results with the
$\Lambda$CDM model by producing a dark energy equation of state which is
approximately equal to $-1$ in the late era, thus requiring it to evolve
infinitesimally for small redshifts, however it could change in principle
drastically in previous cosmological eras of interest. In the end, only a
numerical solution of the Friedmann equation can answer this but at the very
least compatibility with data can be achieved.

At this stage, let us define the statefinder function in terms of redshift as
follows,

$\centering
y_{H}=\frac{H^{2}}{m_{s}^{2}}-\frac{\rho_{DE}}{\rho_{nr}^{(0)}}\,,\@add@centering$
(23)

where the physical quantity $m_{s}^{2}$ is the mass scale, which is given by
the following expression,
$m_{s}^{2}=\frac{\kappa^{2}\rho_{nr}^{(0)}}{3}=H_{0}^{2}\Omega^{(0)}_{m}$.
According to the latest Planck data the present value of the Hubble’s
parameter is $H_{0}=67.4\pm 0.5\frac{km}{Mpc\times sec}$ while, the current
value of the matter density parameter is $\Omega^{(0)}_{m}=0.3153$. Here, the
value based on the CMB shall be used but in general one can follow the same
steps using the Cepheid value and obtain different results. Based on the
statefinder function, the Hubble’s parameter and its derivatives can be
written as,

$\centering
H^{2}=m_{s}^{2}\left(y_{H}+\frac{\rho_{m}}{\rho_{nr}^{(0)}}\right)\,,\@add@centering$
(24) $\centering
HH^{\prime}=\frac{m_{s}^{2}}{2}\left(y_{H}^{\prime}+\frac{\rho_{m}^{\prime}}{\rho_{nr}^{(0)}}\right)\,,\@add@centering$
(25) $\centering H^{\prime
2}+HH^{\prime\prime}=\frac{m_{s}^{2}}{2}\left(y_{H}^{\prime\prime}+\frac{\rho_{m}^{\prime\prime}}{\rho_{nr}^{(0)}}\right).\@add@centering$
(26)

These terms are involved in the expressions of the Ricci and Gauss-Bonnet
topological invariant equations. Although the individual expressions seem to
be quite perplexed, the Friedmann equation can easily be solved numerically.
In the following, we present the current matter density in terms of redshift
which reads,

$\centering\rho=\rho_{nr}^{(0)}(1+z)^{3}(1+\chi(1+z)).\@add@centering$ (27)

One can ascertain whether a particular model is in agreement with the standard
cosmological model $\Lambda$CDM by comparing of the numerical values of the
the extra statefinder functions. These statefinder functions are, the
deceleration parameter $q$, the jerk $j$, the snap $s$ and $Om$ which are
given by the following expressions,

$\displaystyle\centering q\@add@centering$
$\displaystyle=-1-\frac{\dot{H}}{H^{2}}$ $\displaystyle j$
$\displaystyle=\frac{\ddot{H}}{H^{3}}-3q-2$ $\displaystyle s$
$\displaystyle=\frac{j-1}{3\left(q-\frac{1}{2}\right)}$ $\displaystyle Om$
$\displaystyle=\frac{\left(\frac{H}{H_{0}}\right)^{2}-1}{(1+z)^{3}-1}.$ (28)

where $H_{0}$ in this context should be considered as the current value of the
Hubble rate as derived from the numerical analysis of a certain model and not
the parameter taken from the CMB. From a more detailed perspective, the
deceleration parameter and the jerk have the following forms in terms of
redshift,

$\centering q(z)=-1+\frac{(1+z)H^{\prime}}{H}\,,\@add@centering$ (29)
$\centering
j=(1+z)^{2}\left(\left(\frac{H^{\prime}}{H}\right)^{2}+\frac{H^{\prime\prime}}{H}-\frac{2H^{\prime}}{(1+z)H}\right)+1\,.\@add@centering$
(30)

For the case of constant dark energy density, the jerk is trivially equal to
unity in the late era while snap should be zero by definition. Reproducing the
results of the $\Lambda$CDM suggests that the jerk and the snap should evolve
infinitesimally around the aforementioned values. Due to the numerical
solution and the definition of the snap, a more complex form should be
expected and is indeed derived as shown in the following sections. In
addition, $Om$ should currently have a value close to the matter density
parameter $\Omega_{m}$ for the sake of consistency. Finally, the equation of
state parameter $\omega_{DE}$ and the density parameter of the dark energy can
be expressed in terms of the auxiliary function $y_{H}$. These equations read,

$\omega_{DE}=-1+\frac{1+z}{3}\frac{dlny_{H}}{dz},$ (31)
$\Omega_{DE}=\frac{y_{H}}{y_{H}+\rho^{\prime}},$ (32)

where $\rho^{\prime}=\frac{\rho}{\rho^{(0)}_{nr}}$. Taking into consideration
the above equations, the validity of the model requires the numerical values
of the aforementioned parameters for $z=0$ to be consistent with their current
observational values. If this happens then, the model at hand can be in
agreement with the $\Lambda$CDM. For simplicity, the area of study shall focus
on te area of [-0.9,10] for redshift in order to test the model up to the
matter dominated era and also make certain predictions for the evolution of
the universe and see whether the description matches our current
understanding.

##  III Late-time dynamics of scalar-$f(\mathcal{G})$ gravity

In this section the late-time phenomenology of a given scalar-tensor model is
investigated. The function $f(\mathcal{G})$ is defined as,

$\centering
f(\mathcal{G})=\frac{\alpha_{1}}{\mathcal{G}^{2}}+\alpha_{2}\sqrt{\mathcal{G}}\,,\@add@centering$
(33)

while the scalar potential has the following form,

$\centering
V(\phi)=\Lambda^{2}(1-e^{-\sqrt{\frac{2}{3}}\frac{\phi}{M_{P}}})\,,\@add@centering$
(34)

where, $\Lambda$ denotes the cosmological constant. Parameters $\alpha_{i}$
are used for dimensional purposes and take the values $\alpha_{1}=1ev^{12}$
and $\alpha_{2}=1ev^{2}$. Here, two separate contributions are considered, one
being singular however both drive the late time. The scalar potential is used
for the inflationary era. For the scalar potential, an exponential function is
considered. In order to proceed with the overall phenomenology the numerical
values of the free parameters must be determined. Hereafter, the cosmological
constant, the Hubble’s rate and the mass Planck take the following numerical
values $\Lambda=11.895\cdot 10^{-67}ev^{2}$, $H_{0}=1.37187\cdot 10^{-33}ev$,
$M_{P}=1.2\cdot 10^{28}ev$ respectively. It is worth to remind the reader that
the present value of the current Hubble rate is nothing but the Hubble rate
taken from CMB data and expressed in units of energy.

Figure 1: Numerical solutions for the statefinder and the scalar field with
respect to the Planck mass in the area [-0.9,10]. As shown, both are
increasing with time. Figure 2: Evolution of Hubble’s parameter in the context
of scalar-$f(\mathcal{G})$ gravity (orange curve) compared with the
$\Lambda$CDM model (green curve) for the interval [-0.9,10]. In the
aforementioned area of integration, the two functions seem to be in a
relatively good agreement.

Figure 3: Deceleration parameter (black), the jerk (red), the snap (green) and
$Om$ (blue) as functions of redshift.

The first step of our analysis is to set the initial conditions, which are
necessary in order to solve numerically the system of differential equations
of motion (3-5) with respect to $y_{H}(z)$ and $\phi(z)$. For this purpose,
the redshift interval has been chosen as follows z=[-0.9,10] where, the
initial value of redshift is considered as $z=10$. According to this
assumption, the initial value of statefinder function and its first derivative
with respect to redshift are given by the expressions,

$y_{H}(z=10)=\frac{\Lambda}{3m_{s}^{2}}\left(1+\frac{11}{1000}\right),$ (35)
$y_{H}^{\prime}(z=10)=\frac{\Lambda}{3000m_{s}^{2}}.$ (36)

Furthermore, the initial conditions for the scalar field and its first
derivative are respectively, $\phi(z=10)=10^{-16}M_{P}$,
$\phi^{\prime}(z=10)=-10^{-17}M_{P}$.

At this point let us justify the dimensional analysis which is used for the
involved physical quantities in cosmological evolution. According to the
equation (15) the first time derivative of the scalar field has dimensions
[$\dot{\phi}]=ev^{2}$, while the Hubble parameter has dimensions $[H]=ev$.
This assumption implies that the initial value of the time derivative with
respect to the redshift has ev units.

Given the aforementioned initial conditions, the numerical solution of the
system of differential equations of motion leads to plots which are presented
in Fig. 1. It can be interpreted easily that as the redshift increases, the
numerical values of the statefinder $y_{H}$ and of the scalar field decrease.
Additionally in Fig. 2 we present the evolution of the Hubble’s parameter in
the context of scalar-$f(\mathcal{G})$ gravity compared with the $\Lambda$
model. According to the following equation,

$H_{\Lambda}(z)=H_{0}\sqrt{\Omega_{\Lambda}+\Omega_{nr}(1+z)^{3}+\Omega_{r}(1+z)^{4}},$
(37)

the Hubble’s parameter depends on the value of redshift and on the density
parameters $\Omega_{\Lambda}$, $\Omega_{nr}$ and $\Omega_{r}$ for any
cosmological era of interest. In late-time era the numerical values of the
aforementioned parameters are $\Omega_{\Lambda}\simeq 0.681369$ and
$\Omega_{nr}\simeq 0.3153$ based on the latest Planck data [60] while
$\Omega_{r}$ is irrelevant, as suggested by $\chi$ which is of order
$\mathcal{O}(10^{-4})$.

Statefinders | Numerical Results |  $\Lambda$CDM Value
---|---|---
q(0) | -0.522536 | -0.535
j(0) | 1.00087 | 1
s(0) | -0.000284313 | 0
Om(0) | 0.318309 | 0.3153$\pm$0.007
$\omega_{DE}(0)$ | -1.00033 | -1.018$\pm$0.31
$\Omega_{DE}(0)$ | 0.681499 | 0.6847$\pm 0.0073$
Table 1: Numerical comparison of the current values of the extra statefinder
functions. The current values suggest that the model is indeed in agreement
with the $\Lambda$CDM.

Figure 4: Dark energy EoS (yellow) and density parameter (orange) as functions
of redshift. While the dark energy density is smoothly increasing and reaches
unity in the future, the dark energy EoS as expected is infinitesimally
evolving thus rendering the model compatible with $\Lambda$CDM.

One can ascertain if a particular model of interest can be considered as
viable, comparing the numerical statefinder parameters and the dark energy
parameters with the $\Lambda$CDM model. In the context of
scalar-$f(\mathcal{G})$ gravity the current value of the deceleration
parameter is $q(0)\simeq-0.52$ which indicates good agreement with the
$\Lambda$CDM as presented in Table I. In addition, the jerk and the snap
parameters take also numerical values which are in consistency with the
$\Lambda$CDM as the former is almost equal to unity, while the latter can be
considered approximately as zero. The Om statefinder parameter which indicates
the present value of matter density is almost equal to $\Omega_{m}$ as
expected. The dark energy density today is predicted $\Omega_{DE}(0)\simeq
0.68$ which means that the condition $\Omega_{DE}(0)+\Omega m(0)\simeq 1$ is
satisfied given that the radiation today can be considered that has negligible
contribution in late-time cosmological era. In Figures 3 and 4 it is
demonstrated the evolution of the aforementioned statefinder parameters with
respect to redshift for the interval [-0.9,10]. Here, it is showcased that the
deceleration parameter reaches the value of $q=0.5$ for large redshifts and in
the future seems to reach the value of $q=-1$. This implies that the model
smoothly unifies the late-time era with the matter dominated era while in the
future a de Sitter expansion is achieved. Similarly, jerk is infinitesimally
increasing however it resides near the value of $j=1$ as expected. Now, the
snap, as mentioned before, seems to be currently near the expected value of
$s=0$ but decreases with redshift. This decrease is a direct consequence of
the infinitesimal increase of the jerk. It is expected that the snap will
diverge at some value of redshift where the numerical value of the
deceleration parameter shall be exactly $q=0.5$. Contrary to the snap, $Om$
seems to decrease with time while the dark energy density parameter
$\Omega_{DE}$ increases. As a result, the Friedmann constraint
$\Omega_{nr}+\Omega_{r}+\Omega_{DE}=1$ seems to be satisfied approximately by
replacing the non relativistic density parameter with $Om$. Finally, the dark
energy EoS does not currently coincide with the $\Lambda$CDM but in the future
where dark energy will be the only dominant contribution to the Friedmann
equation, the EoS shall obtain the value $\omega_{DE}=-1$ during the de Sitter
expansion. Once again, the infinitesimal evolution is indicative of the
compatibility of the current model with the $\Lambda$CDM.

##  IV Theoretical Framework of $f(\mathcal{G})$ gravity accompanied by a
k-essence scalar field

In the following section it is demonstrated briefly the theoretical framework
of late-time cosmological evolution in the presence of a k-essence scalar
field assisted by $f(G)$ gravity in the absence of scalar potential. K-essence
scalar field involves a nonstandard kinetic term which consists of the
combination of the usual part
$X=\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi$ and a higher
order kinetic term, in this case a squared contribution $X^{2}$ shall be
assumed. The gravitational action for such a model is written as,

$\centering
S=\int{d^{4}x\sqrt{-g}\left(\frac{R}{2\kappa^{2}}-c_{1}X-c_{2}X^{2}+\frac{f(\mathcal{G})}{2}+\mathcal{L}_{matter}\right)}\,,\@add@centering$
(38)

where the auxiliary parameters $c_{1}$ and $c_{2}$ are also arbitrary
constants with the first being dimensionless while the latter has mass
dimensions $[m]^{-4}$ for the sake of consistency. Here, $c_{1}$ is used in
order to have the phantom case corresponding to $c_{1}=-1$ available for the
reader however in the following numerical analysis only a canonical
contribution shall be assumed. By implementing the variation principle with
respect to the metric tensor and the scalar field respectively, the equations
of motion are derived which in this case read,

$\centering\frac{3H^{2}}{\kappa^{2}}=\rho_{m}+\frac{c_{1}}{2}\dot{\phi}^{2}-\frac{3c_{2}}{4}\dot{\phi}^{4}+\frac{\mathcal{G}f_{\mathcal{G}}-f}{2}-12\dot{f}_{\mathcal{G}}H^{3}\,,\@add@centering$
(39)
$\centering-\frac{2\dot{H}}{\kappa^{2}}=\rho_{m}+P_{m}+2c_{1}\dot{\phi}^{2}-4c_{2}\kappa^{4}\dot{\phi}^{4}+4H^{2}(\ddot{f}_{\mathcal{G}}-H\dot{f}_{\mathcal{G}})+8\dot{f}_{\mathcal{G}}H\dot{H}\,,\@add@centering$
(40) $\centering
c_{1}(\ddot{\phi}+3H\dot{\phi})-3c_{2}(\ddot{\phi}+H\dot{\phi})\kappa^{4}\dot{\phi}^{2}=0\,,\@add@centering$
(41)

Our methodology for the investigation of the late-time era is based on the
analysis which is performed in Section II. However, as the system of equations
of motion has been altered due to the inclusion of the quadratic kinetic term
$X^{2}$, the density and the pressure of the dark energy are given now by the
following relations respectively,

$\centering\rho_{DE}=\frac{c_{1}}{2}\dot{\phi}^{2}-\frac{3c_{2}}{4}\dot{\phi}^{4}+\frac{\mathcal{G}f_{\mathcal{G}}-f}{2}-12\dot{f}_{\mathcal{G}}H^{3},\@add@centering$
(42) $\centering
P_{DE}=2c_{1}\dot{\phi}^{2}-4c_{2}\kappa^{4}\dot{\phi}^{4}+4H^{2}(\ddot{f}_{\mathcal{G}}-H\dot{f}_{\mathcal{G}})+8\dot{f}_{\mathcal{G}}H\dot{H}-\rho_{DE}.\@add@centering$
(43)

The rest expressions remain unaffected by the new inclusion in the
gravitational action. This is really convenient as one can use exactly the
same equations as before for a variety of gravitational actions and just
specify properly what is dark matter density and pressure in this context.
Compatibility with observations, if it can be achieved, is studied by solving
the equations of motion numerically no matter the inclusion of additional
terms in the action.

##  V Numerical Analysis and Validity of the model

We present the late-time phenomenology in the presence of a k-essence scalar
field assisted by $f(\mathcal{G})$ gravity. For the shake of simplicity the
case of absent scalar potential is considered, which is a usual assumption for
noncanonical scalar fields with higher order kinetic terms. First of all, the
arbitrary function $f(\mathcal{G})$ which appears in the gravitational action
(38), is now specified as follows,

$f(\mathcal{G})=\frac{\alpha_{1}}{\mathcal{G}}+\alpha_{2}\mathcal{G}^{1/4}$
(44)

where $\alpha_{1}$ and $\alpha_{2}$ are auxiliary parameters introduced for
dimensional purposes but are both replaced with the value $\alpha_{i}=1$ for
simplicity. Once again, our analysis begins from the imposed initial
conditions. For this reason it is assumed that the statefinder function as
well as its first derivative with respect to redshift are given by the
following expressions $y_{H}(z=10)=\frac{1.05\Lambda}{3m_{s}^{2}}$ ,
$y_{H}^{\prime}(z=10)=\frac{\Lambda}{3000m_{s}^{2}}$ respectively.
Furthermore, the initial value of the scalar field and its first derivative
is, $\phi(z=10)=10^{-16}M_{P}$ and $\phi^{\prime}(z=10)=-10^{-17}M_{P}$. In
addition, the cosmological constant and the Planck mass remain exactly the
same as in the previous model. The numerical solution of the system of the
differential equations (39-41) with respect to the statefinder function is
presented in Figure 5. The current value of the statefinder function is
$y_{H}(0)=2.20394$. According to the plot it is abundantly clear the absence
of dark energy oscillations during the cosmological evolution.

Figure 5: Numerical solution of the statefinder function for the interval
[-0.9,10]. In contrast to the previous model, the statefinder function
decreases with time.

In this case, the solutions of the equations of motion seem to decrease with
time compared to the previous model. The difference lies with the initial
condition for statefinder $y_{H}$ for the current model which as shown manages
to alter the properties of dark energy. In Table II, the present numerical
values of the statefinder functions of the given model are demonstrated and as
shown are in consistency with the $\Lambda$CDM standard cosmological model.

Figure 6: Deceleration parameter (black), the jerk (red), the snap (green) and
$Om$ (blue) as functions of redshift in the late-time era. The overall
behaviour seems quite similar between the models however differences in
parameters containing higher derivatives of Hubble are spotted for small
values of redshift.

Figure 7: Dark energy EoS (yellow) and density parameter (orange) as functions of redshift. In the current model, the dark energy EoS is decreasing for large values of redshift. Statefinders | Numerical Results |  $\Lambda$CDM Value
---|---|---
q(0) | -0.531348 | -0.535
j(0) | 0.999532 | 1
s(0) | 0.000151217 | 0
Om(0) | 0.312435 | 0.3153$\pm$0.007
$\omega_{DE}(0)$ | -0.999679 | -1.018$\pm$0.31
$\Omega_{DE}(0)$ | 0.687818 | 0.6847$\pm 0.0073$
Table 2: Numerical values of the statefinder parameters in
k-essence-$f(\mathcal{G})$ gravity compared with the $\Lambda$CDM mdoel.

Here, the deceleration parameter seems to behave exactly as the previous
model. This would be expected since the form of $f(\mathcal{G})$ in the
current model shares similar characteristics as the previous, meaning that
both have a singular term given by an inverse power-law and a power-law part
with exponent ranging in the area $0<n<1$. Due to the fact that the functions
are inherently different, higher derivatives of Hubble such as the jerk are
affected. In Fig.6 the jerk has a different behavior for small redshifts. The
same applies to $Om$ which as shown starts increasing again in the future. The
initial conditions used in both cases are quite similar therefore the
different choices in models, both the $f(\mathcal{G})$ and the k-essence part
over the scalar potential affect this result. The major impact seems to lie
with the EoS of dark energy in Fig.7 where it behaves opposite to the previous
model. Now the eoS increases with the increase of redshift with the change
once again being spotted in the third decimal. Infinitesimal evolution is in
an essence what manages to replicate the $\Lambda$CDM results to a good
extend. The behavior of the dark energy density parameter as expected is
similar to the previous case where it increases with time. Overall, while
similar models are studied, different properties for dark energy may be
obtained if one studies scalar field assisted $f(\mathcal{G})$ gravity models.

##  VI Conclusions

In this paper we investigated the late-time cosmological era in the context of
a given scalar-tensor gravitational theory. Specifically the phenomenological
implications of the inclusion of a scalar field and a function
$f(\mathcal{G})$ in the cosmological evolution are examined. Two models of
interest were presented, the first involved the case of a canonical scalar
field, while the second a scalar field with a standard and a higher order
kinetic term without a scalar potential. In order to simplify our analysis we
performed two variable changes. Firstly, we replaced the cosmic time with
redshift and secondly we introduced the statefinder function in terms of dark
energy instead of Hubble’s parameter. Based on the aforementioned
transformations the system of equations of motions was extracted and it was
solved numerically with to the statefinder function and of the scalar field.
For both models which are investigated, it is proved that compatibility with
the latest Planck data and agreement with the $\Lambda$CDM can be achieved.
Comparing with cosmological models, which arise in the context of $f(R)$
gravitational theories the main benefits of the theory worth to be mentioned.
In this context, dark energy oscillations during the early stages of matter
dominated era are absent even though power-law models with exponents $0<n<1$
are considered. The choice of curvature invariant seems to matter. The scalar
field present is mainly used in order to unify inflation with the late-time.
The interesting part lies with the choice of a linear Ricci scalar, which was
mainly used in order to obtain simler equations of motion and also avoid the
appearance of ghost instabilities which are are typically present in
$f(R,\mathcal{G})$ models. The idea of such models however which are assisted
by scalar fields is quite interesting and could result in interesting
characteristics for dark energy. We hope to address this interesting and
rather challenging topic in a future work.

## References

  * [1] A.G.Riess et al. [Supernova Search Team], Astron. J. 116 (1998), 1009-1038 doi:10.1086/300499 [arXiv:astro-ph/9805201 [astro-ph]].
  * [2] V.K.Oikonomou, Phys. Rev. D 103 (2021) no.12, 124028 doi:10.1103/PhysRevD.103.124028 [arXiv:2012.01312 [gr-qc]].
  * [3] V. K. Oikonomou, Phys. Rev. D 103 (2021) no.4, 044036 doi:10.1103/PhysRevD.103.044036 [arXiv:2012.00586 [astro-ph.CO]].
  * [4] S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 99 (2019) no.10, 104070 doi:10.1103/PhysRevD.99.104070 [arXiv:1905.03496 [gr-qc]].
  * [5] V. K. Oikonomou, Gen. Rel. Grav. 45 (2013), 2467-2481 doi:10.1007/s10714-013-1597-7 [arXiv:1304.4089 [gr-qc]].
  * [6] S. D. Odintsov, V. K. Oikonomou and L. Sebastiani, Nucl. Phys. B 923 (2017), 608-632 doi:10.1016/j.nuclphysb.2017.08.018 [arXiv:1708.08346 [gr-qc]].
  * [7] S. Nojiri, S. D. Odintsov and V. K. Oikonomou, arXiv:1912.13128 [gr-qc].
  * [8] V. Sahni, A. Shafieloo and A. A. Starobinsky, Astrophys. J. 793 (2014) no.2, L40 doi:10.1088/2041-8205/793/2/L40 [arXiv:1406.2209 [astro-ph.CO]].
  * [9] S. Nojiri, S. D. Odintsov and V. K. Oikonomou, Phys. Rept. 692 (2017) 1 doi:10.1016/j.physrep.2017.06.001 [arXiv:1705.11098 [gr-qc]].
  * [10] S. Nojiri, S. D. Odintsov and D. Saez-Gomez, Phys. Lett. B 681 (2009) 74 doi:10.1016/j.physletb.2009.09.045 [arXiv:0908.1269 [hep-th]].
  * [11] S.Capozziello and M.De Laurentis, Phys. Rept. 509 (2011) 167 doi:10.1016/j.physrep.2011.09.003 [arXiv:1108.6266 [gr-qc]].
  * [12] S. Nojiri and S. D. Odintsov, eConf C 0602061 (2006) 06 [Int. J. Geom. Meth. Mod. Phys. 4 (2007) 115] doi:10.1142/S0219887807001928 [hep-th/0601213].
  * [13] S. Nojiri and S. D. Odintsov, Phys. Rept. 505 (2011) 59 doi:10.1016/j.physrep.2011.04.001 [arXiv:1011.0544 [gr-qc]].
  * [14] G.J.Olmo, Int. J. Mod. Phys. D 20 (2011) 413 doi:10.1142/S0218271811018925 [arXiv:1101.3864 [gr-qc]].
  * [15] S. Nojiri and S. D. Odintsov, Phys. Rev. D 68 (2003) 123512 doi:10.1103/PhysRevD.68.123512 [hep-th/0307288].
  * [16] S. Nojiri and S. D. Odintsov, Phys. Lett. B 657 (2007) 238 doi:10.1016/j.physletb.2007.10.027 [arXiv:0707.1941 [hep-th]].
  * [17] S. Nojiri and S. D. Odintsov, Phys. Rev. D 77 (2008) 026007 doi:10.1103/PhysRevD.77.026007 [arXiv:0710.1738 [hep-th]].
  * [18] G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov, L. Sebastiani and S. Zerbini, Phys. Rev. D 77 (2008) 046009 doi:10.1103/PhysRevD.77.046009 [arXiv:0712.4017 [hep-th]].
  * [19] S. Nojiri and S. D. Odintsov, Phys. Rev. D 74 (2006) 086005 doi:10.1103/PhysRevD.74.086005 [hep-th/0608008].
  * [20] S. A. Appleby and R. A. Battye, Phys. Lett. B 654 (2007) 7 doi:10.1016/j.physletb.2007.08.037 [arXiv:0705.3199 [astro-ph]].
  * [21] Y. Zhong and D. Sáez-Chillón Gómez, Symmetry 10 (2018) no.5, 170 doi:10.3390/sym10050170 [arXiv:1805.03467 [gr-qc]].
  * [22] B. Li, J. D. Barrow and D. F. Mota, “The Cosmology of Modified Gauss-Bonnet Gravity,” Phys. Rev. D 76 (2007) 044027 doi:10.1103/PhysRevD.76.044027 [arXiv:0705.3795 [gr-qc]].
  * [23] S. Nojiri and S. D. Odintsov, “Modified Gauss-Bonnet theory as gravitational alternative for dark energy,” Phys. Lett. B 631 (2005) 1 [hep-th/0508049].
  * [24] S. Nojiri, S. D. Odintsov and O. G. Gorbunova, “Dark energy problem: From phantom theory to modified Gauss-Bonnet gravity,” J. Phys. A 39 (2006) 6627 [hep-th/0510183].
  * [25] G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov and S. Zerbini, “Dark energy in modified Gauss-Bonnet gravity: Late-time acceleration and the hierarchy problem,” Phys. Rev. D 73 (2006) 084007 [hep-th/0601008].
  * [26] E. Elizalde, R. Myrzakulov, V. V. Obukhov and D. Saez-Gomez, “LambdaCDM epoch reconstruction from $F(R,\mathcal{G})$ and modified Gauss-Bonnet gravities,” Class. Quant. Grav. 27 (2010) 095007 [arXiv:1001.3636 [gr-qc]].
  * [27] K. Izumi, “Causal Structures in Gauss-Bonnet gravity,” Phys. Rev. D 90 (2014) no.4, 044037 [arXiv:1406.0677 [gr-qc]].
  * [28] V. K. Oikonomou, “Gauss-Bonnet Cosmology Unifying Late and Early-time Acceleration Eras with Intermediate Eras,” Astrophys. Space Sci. 361 (2016) no.7, 211 doi:10.1007/s10509-016-2800-6 [arXiv:1606.02164 [gr-qc]].
  * [29] K. Kleidis and V. K. Oikonomou, Int. J. Geom. Meth. Mod. Phys. 15 (2017) no.04, 1850064 doi:10.1142/S0219887818500640 [arXiv:1711.09270 [gr-qc]].
  * [30] A. Escofet and E. Elizalde, “Gauss-Bonnet modified gravity models with bouncing behavior,” Mod. Phys. Lett. A 31 (2016) no.17, 1650108 doi:10.1142/S021773231650108X [arXiv:1510.05848 [gr-qc]].
  * [31] A. N. Makarenko and A. N. Myagky, “The asymptotic behavior of bouncing cosmological models in $F(\mathcal{G})$ gravity theory,” Int. J. Geom. Meth. Mod. Phys. 14 (2017) no.10, 1750148 doi:10.1142/S0219887817501481 [arXiv:1708.03592 [gr-qc]].
  * [32] A. N. Makarenko, “The role of Lagrange multiplier in Gauss-Bonnet dark energy,” Int. J. Geom. Meth. Mod. Phys. 13 (2016) no.05, 1630006. doi:10.1142/S0219887816300063
  * [33] G. Navo and E. Elizalde, doi:10.1142/S0219887820501625 [arXiv:2007.11507 [gr-qc]].
  * [34] F. Bajardi and S. Capozziello, Eur. Phys. J. C 80 (2020) no.8, 704 doi:10.1140/epjc/s10052-020-8258-2 [arXiv:2005.08313 [gr-qc]].
  * [35] S. Capozziello, C. A. Mantica and L. G. Molinari, Int. J. Geom. Meth. Mod. Phys. 16 (2019) no.09, 1950133 doi:10.1142/S0219887819501330 [arXiv:1906.05693 [gr-qc]].
  * [36] M. Benetti, S. Santos da Costa, S. Capozziello, J. S. Alcaniz and M. De Laurentis, Int. J. Mod. Phys. D 27 (2018) no.08, 1850084 doi:10.1142/S0218271818500840 [arXiv:1803.00895 [gr-qc]].
  * [37] T. Clifton and J. D. Barrow, Class. Quant. Grav. 23 (2006) 2951 doi:10.1088/0264-9381/23/9/011 [gr-qc/0601118].
  * [38] J. D. Barrow and S. Cotsakis, Phys. Lett. B 214 (1988) 515. doi:10.1016/0370-2693(88)90110-4
  * [39] K. Bamba, S. D. Odintsov, L. Sebastiani and S. Zerbini, Eur. Phys. J. C 67 (2010) 295 doi:10.1140/epjc/s10052-010-1292-8 [arXiv:0911.4390 [hep-th]].
  * [40] M. De Laurentis, M. Paolella and S. Capozziello, Phys. Rev. D 91 (2015) no.8, 083531 doi:10.1103/PhysRevD.91.083531 [arXiv:1503.04659 [gr-qc]].
  * [41] A. De Felice, J. M. Gerard and T. Suyama, Phys. Rev. D 82 (2010) 063526 doi:10.1103/PhysRevD.82.063526 [arXiv:1005.1958 [astro-ph.CO]].
  * [42] A. de la Cruz-Dombriz and D. Saez-Gomez, Class. Quant. Grav. 29 (2012) 245014 doi:10.1088/0264-9381/29/24/245014 [arXiv:1112.4481 [gr-qc]].
  * [43] E. Elizalde, S. Nojiri, S. D. Odintsov, L. Sebastiani and S. Zerbini, Phys. Rev. D 83 (2011) 086006 doi:10.1103/PhysRevD.83.086006 [arXiv:1012.2280 [hep-th]].
  * [44] S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 101 (2020) no.4, 044009 doi:10.1103/PhysRevD.101.044009 [arXiv:2001.06830 [gr-qc]].
  * [45] K. Bamba, A. Lopez-Revelles, R. Myrzakulov, S. D. Odintsov and L. Sebastiani, Class. Quant. Grav. 30 (2013) 015008 doi:10.1088/0264-9381/30/1/015008 [arXiv:1207.1009 [gr-qc]].
  * [46] F. Fronimos, Eur. Phys. J. Plus 136 (2021) no.10, 1014 doi:10.1140/epjp/s13360-021-02014-6 [arXiv:2110.00353 [gr-qc]].
  * [47] V. K. Oikonomou, Class. Quant. Grav. 38 (2021) no.19, 195025 doi:10.1088/1361-6382/ac2168 [arXiv:2108.10460 [gr-qc]].
  * [48] S. D. Odintsov, V. K. Oikonomou and F. P. Fronimos, Nucl. Phys. B 958 (2020), 115135 doi:10.1016/j.nuclphysb.2020.115135 [arXiv:2003.13724 [gr-qc]].
  * [49] V. K. Oikonomou and F. P. Fronimos, Class. Quant. Grav. 38 (2021) no.3, 035013 doi:10.1088/1361-6382/abce47 [arXiv:2006.05512 [gr-qc]].
  * [50] V. K. Oikonomou and F. P. Fronimos, Eur. Phys. J. Plus 135 (2020) no.11, 917 doi:10.1140/epjp/s13360-020-00926-3 [arXiv:2011.03828 [gr-qc]].
  * [51] S. D. Odintsov, V. K. Oikonomou and F. P. Fronimos, Annals Phys. 424 (2021), 168359 doi:10.1016/j.aop.2020.168359 [arXiv:2011.08680 [gr-qc]].
  * [52] S. A. Venikoudis and F. P. Fronimos, Gen. Rel. Grav. 53 (2021) no.8, 75 doi:10.1007/s10714-021-02846-8 [arXiv:2107.09457 [gr-qc]].
  * [53] S. A. Venikoudis and F. P. Fronimos, Eur. Phys. J. Plus 136 (2021) no.3, 308 doi:10.1140/epjp/s13360-021-01298-y [arXiv:2103.01875 [gr-qc]].
  * [54] B. P. Abbott et al. [LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR and SKA South Africa/MeerKAT], Astrophys. J. Lett. 848 (2017) no.2, L12 doi:10.3847/2041-8213/aa91c9 [arXiv:1710.05833 [astro-ph.HE]].
  * [55] S. D. Odintsov, V. K. Oikonomou and F. P. Fronimos, Phys. Dark Univ. 29 (2020), 100563 doi:10.1016/j.dark.2020.100563 [arXiv:2004.08884 [gr-qc]].
  * [56] F. P. Fronimos and S. A. Venikoudis, [arXiv:2110.12457 [gr-qc]].
  * [57] A. Munyeshyaka, J. Ntahompagaze and T. Mutabazi, Int. J. Mod. Phys. D 30 (2021) no.07, 2150053 doi:10.1142/S021827182150053X [arXiv:2103.09597 [gr-qc]].
  * [58] S. D. Odintsov, V. K. Oikonomou and F. P. Fronimos, Class. Quant. Grav. 38 (2021) no.7, 075009 doi:10.1088/1361-6382/abe24f [arXiv:2102.02239 [gr-qc]].
  * [59] S.D.Odintsov, V.K.Oikonomou, F.P.Fronimos and K.V.Fasoulakos, Phys. Rev. D 102 (2020) no.10, 104042 doi:10.1103/PhysRevD.102.104042 [arXiv:2010.13580 [gr-qc]].
  * [60] N. Aghanim et al. [Planck], Astron. Astrophys. 641 (2020), A6 [erratum: Astron. Astrophys. 652 (2021), C4] doi:10.1051/0004-6361/201833910 [arXiv:1807.06209 [astro-ph.CO]].