11institutetext: Siemens AG, Munich, Germany11institutetext: Siemens AG,
Munich, Germany
11email: {anna.himmelhuber, stephan.grimm, sonja.zillner, mitchell.joblin,
martin.ringsquandl<EMAIL_ADDRESS>22institutetext: Technical
University of Munich, Munich, Germany

# Combining Sub-Symbolic and Symbolic Methods for Explainability

Anna Himmelhuber 1122 Stephan Grimm 11 Sonja Zillner 11 Mitchell Joblin 11
Martin Ringsquandl 11 Thomas Runkler 1122

###### Abstract

Similarly to other connectionist models, Graph Neural Networks (GNNs) lack
transparency in their decision-making. A number of sub-symbolic approaches
have been developed to provide insights into the GNN decision making process.
These are first important steps on the way to explainability, but the
generated explanations are often hard to understand for users that are not AI
experts. To overcome this problem, we introduce a conceptual approach
combining sub-symbolic and symbolic methods for human-centric explanations,
that incorporate domain knowledge and causality. We furthermore introduce the
notion of fidelity as a metric for evaluating how close the explanation is to
the GNN’s internal decision making process. The evaluation with a chemical
dataset and ontology shows the explanatory value and reliability of our
method.

###### Keywords:

Graph Neural Networks XAI Symbolic Methods Inductive Logic Learning.

##  1 Introduction

Many important real-world data sets come in the form of graphs or networks,
including social networks, knowledge graphs, protein-interaction networks, the
World Wide Web and many more. Graph neural networks are connectionist models
that capture the dependence structure induced by links via message passing
between the nodes of graphs. Unlike standard neural networks, GNNs retain a
state that can represent information from its neighborhood with arbitrary
depth as well as incorporate node feature information [10]. Similarly to other
connectionist models, GNNs lack transparency in their decision-making. Since
the unprecedented levels of performance of such AI methods lead to increasing
use in the daily life of humans, there is an emerging need to understand the
decision-making process of such systems [1]. While symbolic methods such as
inductive logic learning come with explainability, they perform best when
dealing with relatively small and precise data. Sub-symbolic methods such as
graph neural networks are able to handle large datasets, have a higher
tolerance to noise in real world data, generally have high computing
performance and are easier to scale up [4].

Through the increasing popularity and need for explainability in AI, a variety
of explainable models for neural networks are being developed [18]. These
include surrogate models which are interpretable models that are trained to
approximate the predictions of a black box model [8]. Other approaches include
the identification of the most relevant features [18] [11]. The explainer
methods named above, allow the user to relate properties of the inputs to
their output. However, the user is responsible for compiling and comprehending
the explanations, relying on their own implicit form of knowledge and
reasoning about them. Since humans are depending on their background knowledge
and therefore also their biases about the data and its domain, different
explanations about why a model makes a decision may be deduced. Since such
sub-symbolic models are often built for AI researchers, it can make them hard
to understand for non-experts. We strive to go beyond that by justifying
predictions with background or common sense knowledge in a human
understandable way [3]. This is of increased importance, as explainable AI and
with it the widespread application of AI models are more likely to succeed if
the evaluation of these explainer models is focused more on the user’s needs
[17].

We aim to develop a hybrid method by combining GNNs, sub-symbolic explainer
methods and inductive logic learning. This enables human-centric and causal
explanations through extracting symbolic explanations from identified decision
drivers and enriching them with available background knowledge. These are
generated for individual predictions, and are therefore instance-level
explanations. With this method, high-accuracy sub-symbolic predictions come
with symbolic-level explanations, and provide an effective solution for the
performance vs. explainability trade-off.
As far as we know, this is the first work to study integrating a sub-symbolic
explainer with symbolic methods for more human-centric instance-level
explanations. Our fidelity metric indicates how close an explanation is to the
GNN’s internal decision making process. Additionally, the employment of
justifications in our method provides causality that makes use of the
background knowledge.

##  2 Background and Problem Definition

For incorporating explicit domain knowledge into our explanation method on the
side of symbolic representation, we use ontologies expressed in the W3C OWL 2
standard111https://www.w3.org/TR/owl2-overview/ [6] based on the description
logic formalism. In this section we first introduce semantic web ontology,
then revisit the notions of entailment, inductive logic learning and
justifications, followed by graph neural networks and a sub-symbolic explainer
method. Eventually, we define the problem of learning explainer classes by
combining a GNN’s output with inductive logic learning 222For better
readability we will denominate variables represented in ontology form in greek
letters and sub-symbolic graph representations in latin letters..

Semantic Web Ontology The basic constituents for representing knowledge in OWL
are individuals, classes and properties. They are used to form axioms, i.e.
statements within the target domain, and an ontology $\mathcal{O}$ is a set of
axioms to describe what holds true in this domain. The most relevant axioms
for our work are class assertions $\tau(\sigma)$ assigning an individual
$\sigma$ to a class $\tau$, property assertions $\rho(\sigma_{1},\sigma_{2})$
connecting two individuals $\sigma_{1},\sigma_{2}$ by property $\rho$, and
subclass axioms $\tau_{1}\sqsubseteq\tau_{2}$ expressing that class $\tau_{1}$
is a subclass of class $\tau_{2}$. Classes can be either atomic class names,
such as ’Compound‘ or ’Bond‘, or they can be composed by means of complex
class expressions. An example for a complex class expression noted in
Manchester syntax is ’Compound and hasStructure some Nitrogen_Dioxide‘, which
refers to all molecule compounds having some nitrogen dioxide compound in
their structure. For details about all types of axioms and the way complex
concepts are constructed we refer to [6].

(1) |  Carbon $\sqsubseteq$ Atom |  carbons are atoms
---|---|---
(2) |  Hetero_aromatic_5_ring $\sqsubseteq$ Ring_Size_5 $\sqsubseteq$ RingStructure |  hetero-aromatic rings of size 5 are rings of size 5, which are ring structures
(3) |  Nitrogen(feature_100_5) |  feature_100_5 is a nitrogen
(4) |  Compound(graph_100) |  graph_100 is a compound
(5) |  hasAtom(graph_100, feature_100_5) |  graph_100 has atom feature_100_5
Table 1: Example excerpt of $\delta^{Mutag}$.

Example 1. (Mutagenesis Ontology)
As we are combining GNNs and ontologies, graph data has to be available as
triples as well as background knowledge. We chose a chemical domain to test
our method, as it comes with structured background knowledge. The domain
knowledge used in our approach is given by the Mutagenesis ontology
$\delta^{Mutag}$ 333 https://github.com/SmartDataAnalytics/DL-
Learner/tree/develop/examples
/mutagenesis, which is exemplified in Table 1.
Definition 1. (Entailment)
Given ontology $\mathcal{O}$, if axiom $\alpha$ logically follows from
$\mathcal{O}$, as can be derived by a standard OWL reasoner, then we call
$\alpha$ an _entailment_ of $\mathcal{O}$ and write $\mathcal{O}\models\alpha$
444As defined in [13].
Definition 2. (Inductive Logic Learning (ILL)): Given an ontology
$\mathcal{O}$, a set of positive instances $E^{+}$ and a set of negative
instances $E^{-}$, a resulting target predicate class expression $\varepsilon$
is constructed such that $\mathcal{O}\models\varepsilon(\sigma)$ holds for all
individuals $\sigma\in E^{+}$ and does not hold for individuals $\sigma\in
E^{-}$ 555As defined in [16].
In the context of OWL ontologies, ILL attempts to construct class expressions
from an ontology $\mathcal{O}$ and two sets $E^{+},E^{-}$ of individuals that
act as positive and negative examples for being instances of the target class,
respectively. Concretely, we use DL-Learner [16] as the key tool to derive OWL
class expressions to be used for our explanations.
Definition 3. (Justification):
Given an ontology $\mathcal{O}$ and an entailment $\alpha$, the
_justification_ $\mathcal{J}(\mathcal{O},\alpha)$ for $\alpha$ in
$\mathcal{O}$ is a set $\mathcal{J}\subseteq\mathcal{O}$, such that
$\mathcal{J}\models\alpha$ and $\mathcal{J}^{\prime}\not\models\alpha$ for all
proper subsets $\mathcal{J}^{\prime}\subset\mathcal{J}$
Graph Neural Network (GNN)
For a GNN, the goal is to learn a function of features on a graph $G=(V,E)$
with edges $E$ and nodes $V$. The input is comprised of a feature vector
$x_{i}$ for every node $i$, summarized in a feature matrix
$X\in\mathbb{R}^{n\times d_{in}}$ and a representative description of the link
structure in the form of an adjacency matrix $A$. The output of one layer is a
node-level latent representation matrix $Z\in\mathbb{R}^{n\times d_{out}}$,
where $d_{out}$ is the number of output latent dimensions per node. Therefore,
every layer can be written as a non-linear function: $H^{(l+1)}=f(H^{(l)},A)$,
with $H^{(0)}=X$ and $H^{(L)}=Z$, $L$ being the number of stacked layers. The
vanilla GNN model employed in our framework, uses the propagation rule [9]:

$f(H^{(l)},A)=\hat{s}(\hat{D}^{-\frac{1}{2}}\hat{A}\hat{D}^{-\frac{1}{2}}H^{(l)}W^{(l)}),$

with $\hat{A}=A+I$, $I$ being the identity matrix. $\hat{D}$ is the diagonal
node degree matrix of $\hat{A}$, $W^{(l)}$ is a weight matrix for the $l-th$
neural network layer and $\hat{s}$ is a non-linear activation function. Taking
the latent node representations $Z$ of the last layer we define the logits of
node $v_{i}$ for classification task as
$\hat{y}_{i}=\text{softmax}(z_{i}W^{\top}_{c})$, where
$W_{c}\in\mathbb{R}^{d_{out}\times k}$ projects the node representations into
the $k$ dimensional classification space.
Example 2. (GNN Classifications for Mutag Dataset)
For GNN predictions, the dataset Mutag is utilized, which is from a different
source and therefore independent of the Mutagenesis ontology. It contains
molecule graphs and is classified through a 3-layer vanilla Graph
Convolutional Network with 85% accuracy [11]. The molecule graphs $G_{i}$ =
($A_{i}$, $X_{i}$), which are compounds existing out of atoms and bonds, with
certain structures such as carbon rings, can be classified as mutagenic (m) or
nonmutagenic (n) depending on their mutagenic effect on the Gram-negative
bacterium S. typhimurium [12].
Sub-symbolic Explainer Method
The sub-symbolic explainer method takes a trained GNN and its prediction(s),
and it returns an explanation in the form of a small subgraph of the input
graph together with a small subset of node features that are most influential
for the prediction. For their selection, the mutual information between the
GNN prediction and the distribution of possible subgraph structures is
maximized through optimizing the conditional entropy. The explainer method
output is comprised of edge masks $M_{E_{i}}\in\\{0,1\\}^{n\times n}$
$\subset$ $A_{i}$ and node feature masks $M_{X_{i}}\in\\{0,1\\}^{n\times d}$
$\subset$ $X_{i}$, which is used as input to our framework.
Since it is the state-of-the-art method, which outperforms alternative
baseline approaches by 43.0% in explanation accuracy [11], we chose the
GNNExplainer for our framework, but our approach will work with any other
explainer subgraph generation method.
Example 3. (GNNExplainer Output for Mutag Dataset Classifications)
The GNNExplainer is applied to identify the most influential parts of the
respective graph for the classification decision. Figure 1 shows the original
graph, its edge mask $M_{E}$ as identified by the GNNExplainer and the ground
truth for a mutagenic (left) and nonmutagenic (middle) molecule as well as the
identified node feature mask $M_{X}$ (right). It can be seen that the
identified important graph motifs and node features align with some of ground
truth mutagenic properties, as given by [12]. These include ring structures
and the node features $C$, $O$, $N$ and $H$. However, the fact that these
results represent a carbon ring as well as the chemical group $NO_{2}$
(Nitrogen dioxide) is left up to the user for interpretation.

Figure 1: GNNExplainer results for Mutag dataset classifications. Figure
adapted from Figure 4 and Figure 5 in [11].

Definition 4. (Explainer Class Learning)
Given ontology $\bigO$ and a set of graph individuals $\\{\eta_{i}\\}\in\bigO$
666Mapping sub-symbolic graph representations ($X_{i}$, $A_{i}$) $\mapsto$
$\bigO$, resulting in individuals $\eta_{i}$ is specified in Section 3.1 with
their respective classifications $\\{y_{1},y_{2},...y_{i}\\}$ provided by a
GNN for a certain category, we define explainer class learning as inductive
logic learning such that $\eta_{j}|y_{j}=category\ \in E^{+}$ and
$\eta_{k}|y_{k}\not=category\ \in E^{-}$.
$\bigO$ provides the background knowledge for inductive logic learning, and
the classification decision by the GNN provides the positive and negative
examples in order to learn explainer classes. We also define a metric called
fidelity metric (to be specified in Section 3) for quantitative measurement.
The higher the fidelity metric, the higher the reliability of the entailed
explainer class.

##  3 Combining Sub-Symbolic and Symbolic Methods

###  3.1 Explainer Class Learning

Figure 2: Learning Explainer Classes process flow.

We are proposing a hybrid method, within which the coupling of the sub-
symbolic explainer method GNNExplainer with the symbolic DL-Learner is used to
explain GNN instance-level predictions. Our approach is shown for a graph
classification task, but would equally apply to node classification or link
prediction. The process flow of learning explainer classes can be seen in
Figure 2. Firstly, a GNN is trained on and applied to training and testing
data and subsequently the sub-symbolic explainer method GNNExplainer is
applied to all generated predictions, as can be seen in in Figure 2 (Step 1
and Step 2). Secondly, to create explainer classes for the GNN decision making
process, DL-Learner is applied for a specific predicted category, with
positive and negative examples labelled accordingly through $y_{i}$ (Step 4).
The background knowledge used by the DL-Learner to learn explainer classes is
comprised of the adjacency matrices $A_{i}$ and node feature matrices $X_{i}$,
edge masks $M_{E_{i}}$ and node feature masks $M_{X_{i}}$ and domain knowledge
$\delta$. As the DL-Learner can only process ontologies, the matrices are
mapped to an ontology (Step 3) through $\lambda$ as detailed below:
Extraction and Mapping Step
A set of graphs detailed in their associated matrices $A_{i}$ and $X_{i}$
777Their size is dependent on the number of layers used by the GNN, to keep
the consistency in coupling the sub-symbolic with the symbolic method. are
modelled as set of individuals $\\{\eta_{i}\\}$. Their edges and node features
are extracted from $A_{i}^{\prime}s$ and $X_{i}^{\prime}s$ edge and feature
lists and modelled as set of individuals $\\{\upsilon_{j}\\}$ and
$\\{\chi_{k}\\}$. If there are graph-specific structures common in the
respective domain, such as certain motifs, e.g. a ring structure, the set of
possible structures $\\{structure_{z}\\}$ along with their extraction
functions $\\{\gamma_{z}(A_{i},X_{i})\\}$ is defined and mapped through
mapping function $S:\\{structure_{z}\\}$ $\mapsto$
$\\{\gamma_{z}(A_{i},X_{i})\\}$.
If $structure_{1}$ is contained in $(A_{i},X_{i})$, extraction function
$\gamma_{1}(A_{i},X_{i})$ returns all individuals contained in the structure.
The found structures are modelled as a set of individuals $\\{\psi_{g}\\}$. To
assign all individuals their type declarations and roles, a set of roles
$\\{\rho_{v}\\}$ and type declarations $\\{\tau_{w}\\}$ as well as further
mapping functions based on domain knowledge $\delta$ are needed. Defining
these sets and mapping functions has been done as a one-time manual step, with
their complexity depending on the domain.
$P:\\{\eta_{i}\\}\times(\\{\upsilon_{j}\\}\cup\\{\chi_{k}\\}\cup\\{\psi_{g}\\})\mapsto\\{\rho_{v}\\}$,
maps a pair of individuals to their role.
$T:(\\{\eta_{i}\\}\cup\\{\upsilon_{j}\\}\cup\\{\chi_{k}\\}\cup\\{\psi_{g}\\}\mapsto\\{\tau_{w}\\}$
maps individuals to their types. All extracted individuals, roles and type
declarations are added as axioms to ontology $\bigO$ through function
$AddAxiom(\bigO,axiom)$ as is shown in Algorithm 1. Therefore, $\lambda$ is
defined as $\lambda(A_{i},X_{i},T,P,S)\mapsto\bigO$. Equivalently, $\lambda$
is carried out for all corresponding sub-symbolic explainer subgraphs with
their associated edge masks $M_{E_{i}}$ and node feature masks $M_{X_{i}}$,
with the set of explainer graphs modelled as individuals $\eta\\_sub_{i}$.
Additionally, mapping function $\mu$ is defined as bijective function, as is
shown in Algorithm 1. This function is needed for the fidelity calculation.
Function $\mu$ is defined in such a way, that if the input, e.g.
${\sigma_{1}}$, doesn’t map to anything, ${\sigma_{1}}$ will be returned as
output.
Example 4. (Mapping Mutag Dataset with Mutagenesis Ontology)
The mapping functions $S^{Mutag}=\\{Azanide:\gamma_{Azanide}$ ,
$Methyl:\gamma_{Methyl},...\\}$,
$R^{Mutag}=\\{(\eta_{i},\upsilon_{j}):hasBond,(\eta_{i},\chi_{k}):hasAtom,...\\}$
and
$T^{Mutag}=\\{\eta_{i}:Compound,\upsilon_{j}:Bond,\chi_{k}:Carbon,...\\}$ are
defined based on domain terminology $\delta^{Mutag}$. For example, from
molecule graph $G_{1}$ with associated matrices $X_{1}$ and $A_{1}$, the edge
individuals edge_1_2, edge_1_3, etc., are modelled. For extracting structure
$Methyl$ ($CH_{3}$), which is defined as containing one carbon atom bonded to
three hydrogen atom, function $\gamma_{Methyl}(A_{1},X_{1})$ is employed. All
accruing axioms are added to the ontology $\bigO^{Mutag}$. Through $\mu$, the
set of edges forming the identified structure, e.g. { edge_1_2, edge_1_3,
edge_1_4 } is mapped to the individual structure_1_1_1.

Data: Set of graphs with adjacency matrices $A_{i}$, feature matrices $X_{i}$,
mapping functions for type declarations $T(\sigma)$, roles
$P(\sigma_{1},\sigma_{2})$ and structures $S(x)$

Result: $\bigO$, $\mu$

$\bigO:\\{\\}$

foreach  _graph in range(i)_ do

AddAxiom($\bigO$, $T(graph)$ ($\eta_{graph}$))

foreach  _edge in Edgelist( $A_{graph})$_ do

AddAxiom($\bigO$, $T(edge)$($\upsilon_{edge\\_graph}$) )

AddAxiom($\bigO$, $P(graph,edge)(\eta_{graph},\upsilon_{edge\\_graph})$ )

end foreach

foreach  _feature in Featurelist( $X_{graph}$)_ do

AddAxiom($\bigO$, $T(feature)(\chi_{feature\\_graph})$ )

AddAxiom($\bigO$, $P(graph,feature)(\eta_{graph},\chi_{feature\\_graph})$ )

end foreach

foreach  _structure in $\\{structure_{z}\\}$_ do

if  _$S(structure)(A_{graph},X_{graph})\ not\ \textbf{None}$_ then

foreach  _number in range(count( $S(structure)(A_{graph},X_{graph})$))_ do

AddAxiom($\bigO$, $T(structure)(\psi_{graph\\_structure\\_number})$ )

AddAxiom($\bigO$,
$P(graph,structure)(\eta_{graph},\psi_{graph\\_structure\\_number}))$

$\mu:S(structure)(A_{graph},X_{graph})\mapsto\psi_{graph\\_structure\\_number}$

end foreach

end if

end foreach

end foreach

Algorithm 1 Graph Structure Extraction $\lambda$

According to the GNN’s classifications positive and negative examples of
graphs are distinguished and explainer classes are learned. The background
knowledge is the ontology $\bigO$ = $\delta$ $\cup$
$\lambda(A_{i},X_{i},T,P,S)$ $\cup$ $\lambda(M_{E_{i}},M_{X_{i}},T,P,S)$. We
differentiate between two types of explainer classes:

Input-Output Explainer Classes
Given Def. 4, background knowledge $\delta$ $\cup$
$\lambda(A_{i},X_{i},T,P,S)$, $\eta_{i}|y_{i}=category\in E^{+}$ and
$\eta_{i}|y_{i}\not=category\in E^{-}$, a set of Input-Output Explainer
Classes $\\{\phi_{n}^{category}\\}$ are learned. Input-output explainer
classes are candidate explanations, that capture the global behavior of a GNN
through investigating what input patterns can lead to a specific class
prediction, comparable to the input-output mapping approach in [15].
Importance Explainer Classes
Given Def.4, background knowledge $\delta$ $\cup$
$\lambda(M_{E_{i}},M_{X_{i}},T,P,S)$, $\eta\\_sub_{i}|y_{i}=category\in E^{+}$
and $\eta\\_sub_{i}|y_{i}\not=category\in E^{-}$, a set of Importance
Explainer Classes $\\{\varphi_{m}^{category}\\}$ are learned. Importance
Explainer classes show which edges, nodes, features and motifs are important
for the GNN to predict a certain class. These class expressions represent the
inner workings of a GNN, by incorporating the output of the sub-symbolic
explainer.

###  3.2 Explainer Class Application for Instance-Level Explanations

The pool of possible explainer classes for all categories as learned in
Section 3.1, consisting of $\\{\phi_{n}\\}$ and $\\{\varphi_{n}\\}$, are used
in the application step to generate instance-level explanations through
explainer class entailment and justification steps.
Explainer Class Entailment
Given Def. 1., a set of explainer classes $\\{\phi_{n}^{category}\\}$ and
$\\{\varphi^{category}\\}$, ontology $\bigO$ and individual $\eta_{j}$
classified as category, entailments for $\eta_{j}$ are generated. By doing so,
we check if the learned overall decision-making pattern of the GNN applies to
a specific instance. For all available explainer classes, entailments for a
specific individual $\eta_{j}$ are generated. It is possible, that several
entailments hold, just as it is possible that a classification decision of
$G_{j}$ is based on several different factors. The set of entailments for
$\eta_{j}$ is given by $C_{Exp}(\eta_{j})=\\{\phi\ |\ \bigO\models\
\phi^{category}(\eta_{j})\\}\cup\\{\varphi\ |\ \bigO\models\
\varphi^{category}(\eta_{j})\\}$.
Definition 5. (Entailment Frequency)
Given an ontology $\bigO$, explainer class $\phi^{category}_{i}$ and a set of
indivdiuals $\\{\eta_{i}\\}$, we define the entailment frequency as the number
of entailments for $|\\{\eta\in\\{\eta_{i}\\}:\bigO\models\
\phi_{i}^{category}(\eta)\\}|$ over the number of instances
$|\\{\eta_{i}\\}|$.
The entailment frequency gives insight over the generality or specificity of
explainer classes and representing the average frequency with which a certain
explainer class is entailed.
Explainer Class Entailment Justification
Given $\bigO$ and entailment $\bigO\models\ \phi_{i}^{category}(\eta_{j})$,
justification $\mathcal{J}$($\bigO$, $\phi_{i}^{category}$($\eta_{j}$)) is
generated. The number of generated axioms gives some insight about the level
of domain knowledge employed. As there can be several justifications for an
entailment, we limit them to only one. It is not in the scope of this paper to
determine which justification would provide the best explanation, but since a
shorter justification tends to be more efficient, the justification with the
minimum number of axioms is chosen.
Example 5. (Justification for Mutag Explainer Class)
Table 2 shows an example justification for the entailment $\bigO^{Mutag}\
\models\ \phi_{8}^{m}(\eta_{1})$, which contributes to a meaningful
explanation, as it carries causal information present in expert knowledge
about the conclusion.

(1) |  $\eta_{1}$ hasStructure structure_1_1_1
---|---
(2) |  structure_1_1_1 Type Hetero_aromatic_5_ring
(3) |  Hetero_aromatic_5_ring SubClassOf Ring_size_5
(4) |  $\phi_{8}^{m}$ EquivalentTo hasStructure some Ring_size_5
Table 2: Example justification $\mathcal{J}$($\bigO^{Mutag}$,
$\phi_{8}^{m}$($\eta_{1}$)).

Fidelity Calculation
Fidelity is defined as the measure of the accuracy of the student model (DL-
Learner) with respect to the teacher model (GNN). High fidelity is therefore
fundamental, whenever a student model is to be claimed to offer a good
explanation for a teacher model. Without high fidelity, an apparently
perfectly good explanation produced by an explainable system is likely not to
be an explanation of the underlying sub-symbolic system which it is expected
to explain [21]. We calculate Fidelity as follows:

$Fidelity(\phi_{i},\eta_{j})=\frac{|\mu^{-1}(ind(\mathcal{J}(\bigO,\phi_{i}(\eta_{j}))))\cap\eta\\_sub_{j}|}{|\mu^{-1}(ind(\mathcal{J}(\bigO,\phi_{i}(\eta_{j}))))|},$

where $ind()$ is a function that collects all individuals that are provable
instances of a set of axioms. The denominator equals the count of the set of
edges or node features that have to be part of $\eta_{i}$, for the entailment
of explainer class $\phi_{i}$ to hold. The fidelity metric is defined as the
overlap of the sub-symbolic explainer output with the entailed explainer
classes, as can be seen in Figure 3, which means that the effectiveness of the
sub-symbolic explainer method in representing the GNN decision making is
therefore assumed.
Example 6. (Fidelity for Explainer Class hasStructure some Methyl)
As the explainer classes are represented through axioms, e.g. $\phi_{2}^{n}$ =
hasStructure some Methyl, we apply the justification mechanism to arrive at
the axioms containing the corresponding individual(s) for the specific example
$\eta_{1}$, such as $\eta_{1}$ hasStructure structure_1_2_1 $\in$
$\mathcal{J}(\phi_{2}^{n},\eta_{1})$. Since there might be a multiplicity of
individuals, function $ind(\mathcal{J}(\bigO,\phi_{2}(\eta_{1}))$ is applied,
which collects all individuals that are provable instances of the
justification. These individuals are then inversely mapped ($\mu^{-1}$) to
their corresponding set of individuals, in this example { edge_1_2, edge_1_3,
edge_1_4 }. In case there is no corresponding set of individuals, the inverse
mapping simply returns the given individual. For the numerator, we count the
overlap of the identified set of individuals with the individuals in
$\eta\\_sub_{i}$, the subgraph identified by the GNNExplainer.
Definition 6. (Final Explanation)
Given the set of entailments, that hold for $\eta_{j}$, we define the final
explanation E($\eta_{j}$) as the set of the respective justifications
$E(\eta_{j})={\\{\mathcal{J}(\bigO,C(\eta_{j}))\\}\mid C\in
C_{Exp(\eta_{j})}}$.
Example 7. (Molecule Graph $G_{1}$) In Figure 3, the final explanation for the
classification of molecule graph $G_{1}$ as mutagenic can be seen, complete
with justifications and fidelity score.

Figure 3: Final explanation for molecule graph $G_{1}$, which has been
classified as mutagenic.

##  4 Evaluation

Experiment Setting. We used a subset of 530 molecule graphs as training data
to learn explainer classes, and 800 molecule graphs as testing data. The
graphs have been classified by a 3-layer vanilla Graph Convolutional Network.
All molecule graphs come with adjacency matrices $A_{i}^{Mutag}$, and feature
matrices $X_{i}^{Mutag}$ and their corresponding GNNExplainer importance masks
($M_{E_{i}}^{Mutag}$ and $X_{E_{i}}^{Mutag}$), equally split between mutagenic
and nonmutagenic classifications. The DL-Learner can create arbitrarily many
class expressions, functioning as explainer classes, which are ordered by
predictive accuracy (number of correctly classified examples divided by the
number of all examples). We are taking a cut-off point of $>50\%$ predictive
accuracy, as an explainer class with less than $50\%$ predictive accuracy,
wouldn’t represent a pattern for mutagenic classification decisions but rather
the opposite, and v.v. for nonmutagenic classification decisions 888All
experimental data, code and results are available from https://github.com/XAI-
sub-symbolic/Combining-Sub-Symbolic-Explainer-Methods-with-SWT..
Explainer Classes
The generated pool of explainer classes provides a total of 14 explainer
classes for mutagenic and 12 explainer classes for nonmutagenic
classifications. All the comprehensible explanation for mutagenic
classification decisions that can be identified and interpreted from the
GNNExplainer output (see Section 2), have been learnt by the DL-Learner. These
include

$\phi_{2}^{m}$ = hasStructure some Carbon_6_ring,
$\phi_{7}^{m}$ = hasStructure some Nitrogen_dioxide,
$\varphi_{1}^{m}$ = hasAtom some Carbon,
$\varphi_{2}^{m}$ = hasAtom some Hydrogen,
$\varphi_{3}^{m}$ = hasAtom some Nitrogen,
$\varphi_{4}^{m}$ = hasAtom some Oxygen,

along with several others, which have not been identified by the GNNExplainer.
The explainer class $\phi_{6}^{m}$ = hasStructure some Phenanthrene is a
compelling example for the effectiveness of our hybrid approach, as
Phenanthrene is a strong indicator for mutagenic potency [12], but isn’t
identifiable in the GNNExplainer output. This shows that our hybrid method can
identify and verbalize decision-making processes of the GNN, which a
comprehensible sub-symbolic explainer system, whose output might not be easily
understood and interpreted by a user, is missing.

Entailment Frequency

Explainer Class Type |  Number |  Avg. Pred. Acc. (SD) |  Avg. Entailment Rate (SD) |  Avg. Fidelity (SD)
---|---|---|---|---
$\phi_{n}^{m}$ |  1,…,10 |  0.56 (0.04) |  0.64 (0.3) |  0.88 (0.12)
$\phi_{n}^{n}$ |  1,…,5 |  0.59 (0.03) |  0.09 (0.04) |  0.82 (0.12)
$\varphi_{n}^{m}$ |  1,…,4 |  0.77 (0.06) |  0.86 (0.15) |  0.99 (0.01)
$\varphi_{n}^{n}$ |  1,…,7 |  0.56 (0.01) |  0.41 (0.25) |  0.81 (0.05)
Table 3: Input-output and importance explainer classes with avg. pred.
accuracy (DL-Learner), entailment rate and fidelity with their respective
standard deviations (SD).

The entailment frequency gives us insight over the generality or specificity
of explainer classes. As can be seen in Table 3 (Avg. Entailment Rate), there
is a wide range of entailment rates. Some explainer classes, e.g.
$\phi_{4}^{m}$ = hasAtom some Carbon always apply, while others are quite
rare, such as $\phi_{4}^{n}$ = hasAtom some Phosphorus, that comes with only a
4% entailment rate. As expected, we have an overall lower entailment rate for
nonmutagenic explainer classes, as the there are also less distinct factors
indicating nonmutagenicity [12]. Most nonmutagenic classifications come with
about 3 entailments, while mutagenic classifications come with more than 5
entailments on average. This is due to a lower generality of the explainer
classes, which implies that such an explainer class only applies to specific
instances. This notion is also confirmed by the lower average predictive
accuracy of the DL-Learner results for nonmutagenic (57%) as opposed to
mutagenic (63%) explainer classes, as can be seen in Table 3 (Avg. Pred. Acc).
The predictive accuracy of the DL-Learner is defined as the number of
correctly classified examples divided by the number of all examples [14].

Explanation Fidelity
Fidelity gives the user a measure of reliability of the explanation, with the
average fidelity ranging from 64% for $\phi_{5}^{n}$ = hasStructure some
Carbon_5_ring to 100% for e.g. $\varphi_{2}^{m}$ = hasAtom some Hydrogen.
While an explainer class with an average fidelity of 64% might still give the
user some insight, its explanatory value cannot be considered as reliable as
for an explainer class with a higher fidelity. An explainer class, that has a
low generality, meaning it is rarely applied to explain a classification, can
nonetheless come with a high fidelity such as $\phi_{4}^{n}$ (100%). This
suggests that also low generality explainer classes can be valuable for
specific instances.

We can observe a positive correlation of 88% between the average fidelity and
predictive accuracy for $\\{\varphi_{n}\\}$ and of 50% between the average
fidelity and $\\{\varphi_{n}\\}$ $\cup$ $\\{\phi_{n}\\}$, signalising the
effectiveness of representing the sub-symbolic decision-making process with
the DL-Learner. As the predictive accuracy of the output given by the DL-
Learner is the metric on which we base our choice of explainer classes
included in the pool, the correlation with the fidelity indicates that this
approach leads to reliable explanations.

Explainability of sub-symbolic methods is desirable not only to justify
actions taken based on the predictions made by the system, but also to
identify false predictions. Therefore, it is also important to evaluate our
method based on its ability to not generate explanation for wrong predictions
and therefore validating them. Table 4 shows the difference in entailments for
the correctly classified (true positives TP) and incorrectly classified graphs
(false positives FP). We can see, that the average fidelity for entailments is
30 percentage points lower for mutagenic FP than mutagenic TP, and 38
percentage points for nonmutagenic FP. While this might not be sufficient to
clearly identify a wrong classification, it indicates the validity of the
fidelity metric, as it is significantly lower for explainer classes applied to
incorrect classification.

|  $TP^{m}$ |  $FP^{m}$ |  $TP^{n}$ |  $FP^{n}$
---|---|---|---|---
Number of instances |  371 |  29 |  374 |  26
Average fidelity |  0.96 |  0.66 |  0.82 |  0.44
Table 4: Average fidelity for true positives and false positives.

Justification Axioms
Through justifications we provide causality for explanations, based on domain
knowledge. The ontology $\delta^{Mutag}$ utilized has little structural depth
as can be seen in the example excerpt in Table 1. Nonetheless, there is a
minimum of 3 axioms for all entailments. For 20% of explainer classes, 4 axiom
justifications and for 8% of explainer classes, 5 axiom justifications are
generated. This means, that for all explanations generated, the explanations
carry some causal information about the conclusion, supported by expert
knowledge.

###  4.1 Comparison of our Hybrid Method with DL-Learner Explanations and
Input-Output Explanations

DL-Learner: Classifications along with corresponding explanations can be
generated by only using a symbolic classifier such as the DL-Learner. When
comparing this purely symbolic approach with our hybrid method, we find that
using only the DL-Learner comes with significantly lower prediction accuracy
and also explanatory value. The predictive accuracy of the GNN using the same
subset of training data is $78\%$, so considerably above the the DL-Learner
result, as shown below. When applying the DL-Learner to carry out
classifications, we are restricted to only one classifier. This means, even if
we allow more complex class expressions, we only have one explanation for the
target predicate mutagenic:

hasStructure some Nitrogen_dioxide or hasThreeOrMoreFusedRings
value true (pred. acc.: 65.76%)

GNN with Input-Output Explanations: We want to look at the benefits of
integrating a sub-symbolic explainer into our framework, as opposed to
explaining GNN predictions with only the input-output matching method as done
in e.g. [15]. We can see, that for some explainer classes such as
$\phi_{3}^{m}=\varphi_{3}^{m}$ = hasAtom some Nitrogen, we have overlap of the
importance explainer classes with the input-output explainer classes. However,
the importance explainer classes come with a significantly higher predictive
accuracy of 77% as can be seen in Table 3, indicating their significance for
the classification decision. For the nonmutagenic classifications, explainer
class $\varphi_{2}^{n}$ = hasStructure some Carbon_6_ring, which is equivalent
with the ground truth as shown in Figure 1, wouldn’t have been included in
$\phi_{n}^{n}$. Here, we can clearly see the added benefit of generating
explainer classes from the GNNExplainer as opposed to only observing the
input-output behaviour of a GNN. The main benefit of including such a sub-
symbolic explainer, however, is the provision of the fidelity metric. Without
such a metric there is no means to quantify the reliability of the
explanation. These results justify the strategy of using a hybrid method.

###  4.2 Deeper Integration of GNNs with Domain Knowledge

We carried out an initial integration of sub-symbolic and symbolic methods, by
mapping and integrating the GNN input and GNNExplainer output to and with the
available domain knowledge. A deeper integration could be reached through
integrating available domain knowledge $\delta$ into the GNN before training.
As the domain knowledge $\delta$ and the input graphs $G_{i}$ are from
different sources, they are independent. It was therefore not known if their
integration could significantly worsen the GNN classification results. Initial
results, where we included common molecule structure from domain knowledge
$\delta^{Mutag}$ as a simple binary vector into the feature matrices
$X_{i}^{Mutag}$, show that the overall prediction accuracy of the GNN only
decreases insignificantly by 2 percentage points, which is a promising first
result. It indicates that the domain knowledge, and the explanations generated
with it, don’t contradict the decision making process of the GNN.

##  5 Related Work

Explainable AI including model-level interpretation and instance-level
explanations have been the focus of research for years [3]. In this section we
first give an overview for explainable AI for Graph Neural Networks and then
for using symbolic methods to explain sub-symbolic models.
Sub-Symbolic Explainer Methods Current work towards explainable GNNs attempts
to convert approaches initially designed for Convolutional Neural Networks
(CNNs) into graph domain [18]. The drawback of reusing explanation methods
previously applied to CNNs are their inability to incorporate graph-specific
data such as the edge structure. Another method, a graph attention model,
augments interpretability via an attention mechanism by indicating influential
graph structures through learned edge attention weights [7]. It cannot,
however, take node feature information into account and is limited to a
specific GNN architecture. To overcome these problems, [11] created the model-
agnostic approach GNNExplainer, that finds a subgraph of input data which
influence GNNs predictions in the most significant way by maximizing the
subgraph’s mutual information with the model’s prediction.
Explanations with symbolic methods A different type of explainability method
tries to integrate ML with symbolic methods. The symbolic methods utilized
alongside Neural Networks are quite agnostic of the underlying algorithms and
mainly harness ontologies and knowledge graphs [19]. One approach is to map
network inputs or neurons to classes of an ontology or entities of a knowledge
graph. For example, [15] map scene objects within images to classes of an
ontology. Based on the image classification outputted by the Neural Network,
the authors run ILP on the ontology to create class expressions that act as
model-level explanations. Furthermore, [20] learn a mapping between individual
neurons and domain knowledge. This enables the linking of a neuron’s weight to
semantically grounded domain knowledge. A ontology-based approach for human-
centric explanation of transfer learning is proposed by [2]. While there is
some explanatory value to these input-output methods, they fail to give
insights into the inner workings of a graph neural network and cannot identify
which type of information was influential in making a prediction. This work
bridges this gap by combining the advantages of both approaches is among the
first to study the coupling of a sub-symbolic explanation method with symbolic
methods.

##  6 Conclusion

In this paper, we addressed the problem of grounding explanations in domain
knowledge while keeping them close to the decision making process of a GNN. We
showed that combining sub-symbolic with symbolic methods can generate reliable
instance-level explanations, that don’t rely on the user for correct
interpretation. We tested our hybrid framework on the Mutag dataset mapped to
the Mutagenesis ontology, to evaluate its explanatory value, its
practicability and the validity of the idea. We used data from a chemical
domain, as it comes with complex domain knowledge that is universally accepted
and can therefore be considered as ground truth when evaluating explanations.
Our results show, that there are significant advantages of our hybrid
framework over only using the sub-symbolic explainer, where the output is
susceptible to biased or faulty interpretations by the user. Equally, there
are advantages of our hybrid method over a purely symbolic method such as ILL,
as it comes with significantly higher accuracy, while for an input-output
method, the decision-making process of the neural network isn’t considered and
there are no means to validate the reliability of the explanations. In future,
we will evaluate how our hybrid framework compares for different datasets.
Furthermore, we will analyze the effect on explanations when the coupling of
available domain knowledge with GNNs is deepened before training.

## References

  * [1] Arrieta, Alejandro Barredo, et al. “Explainable Artificial Intelligence (XAI): Concepts, taxonomies, opportunities and challenges toward responsible AI.” Information Fusion 58 (2020): 82-115
  * [2] Chen, Jiaoyan, et al. “Knowledge-based transfer learning explanation.” Sixteenth International Conference on Principles of Knowledge Representation and Reasoning. 2018.
  * [3] Biran, Or, and Courtenay Cotton. “Explanation and justification in machine learning: A survey.” IJCAI-17 workshop on explainable AI (XAI). Vol. 8. No. 1. 2017.
  * [4] Ilkou, Eleni, and Maria Koutraki. “Symbolic Vs Sub-symbolic AI Methods: Friends or Enemies?.” CIKM (Workshops). 2020.
  * [5] Tiddi, I. “Foundations of explainable knowledge-enabled systems.” Knowledge Graphs for eXplainable Artificial Intelligence: Foundations, Applications and Challenges 47 (2020): 23.
  * [6] McGuinness, Deborah L., and Frank Van Harmelen. “OWL web ontology language overview.” W3C recommendation 10.10 (2004): 2004.
  * [7] Veličković, Petar, et al. “Graph attention networks.” ICLR (2018).
  * [8] Ribeiro, Marco Tulio, Sameer Singh, and Carlos Guestrin. “Why should i trust you? Explaining the predictions of any classifier.” Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining. 2016.
  * [9] Kipf, Thomas N., and Max Welling. “Semi-supervised classification with graph convolutional networks.” ICLR (2017).
  * [10] Zhou, Jie, et al. “Graph neural networks: A review of methods and applications.” AI Open 1 (2020): 57-81.
  * [11] Ying, Rex, et al. “Gnnexplainer: Generating explanations for graph neural networks.” Advances in neural information processing systems 32 (2019): 9240.
  * [12] Debnath, Asim Kumar, et al. “Structure-activity relationship of mutagenic aromatic and heteroaromatic nitro compounds. correlation with molecular orbital energies and hydrophobicity.” Journal of medicinal chemistry 34.2 (1991): 786-797.
  * [13] Horridge, Matthew, et al. “Understanding Entailments in OWL.” OWLED. 2008.
  * [14] Lehmann, Jens, et al. “Class expression learning for ontology engineering.” Journal of Web Semantics 9.1 (2011): 71-81.
  * [15] Sarker, Md Kamruzzaman, et al. “Explaining trained neural networks with semantic web technologies: First steps.” arXiv preprint arXiv:1710.04324 (2017).
  * [16] Lehmann, Jens. “DL-Learner: learning concepts in description logics.” The Journal of Machine Learning Research 10 (2009): 2639-2642.
  * [17] Miller, Tim, Piers Howe, and Liz Sonenberg. “Explainable AI: Beware of inmates running the asylum or: How I learnt to stop worrying and love the social and behavioural sciences.” IJCAI Workshop on Explainable Artificial Intelligence (XAI). 2017.
  * [18] Pope, Phillip E., et al. “Explainability methods for graph convolutional neural networks.” Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. 2019.
  * [19] Seeliger, Arne, Matthias Pfaff, and Helmut Krcmar. “Semantic Web Technologies for Explainable Machine Learning Models: A Literature Review.” PROFILES/SEMEX@ ISWC 2465 (2019): 1-16.
  * [20] Selvaraju, Ramprasaath R., et al. “Choose your neuron: Incorporating domain knowledge through neuron-importance.“ Proceedings of the European conference on computer vision (ECCV). 2018.
  * [21] Garcez, Artur d’Avila, and Luis C. Lamb. “Neurosymbolic AI: the 3rd Wave.” arXiv preprint arXiv:2012.05876 (2020).


# Attractors in STGR with Boundary Correction

Andronikos Paliathanasis<EMAIL_ADDRESS>Institute of Systems Science,
Durban University of Technology, Durban 4000, South Africa  Departamento de
Matemáticas, Universidad Católica del Norte, Avda. Angamos 0610, Casilla 1280
Antofagasta, Chile

###### Abstract

We investigate the asymptotic behavior of the cosmological field equations in
Symmetric Teleparallel General Relativity, where a nonlinear function of the
boundary term is introduced instead of the cosmological constant to describe
the acceleration phase of the universe. Our analysis reveals constraints on
the free parameters necessary for the existence of an attractor that
accurately represents acceleration. However, we also identify asymptotic
solutions depicting Big Rip and Big Crunch singularities. To avoid these
solutions, we must impose constraints on the phase-space, requiring specific
initial conditions.

Cosmology; Symmetric teleparallel; boundary term; $f\left(Q,C\right)$-theory;
attractors

##  I Introduction

Recent cosmological observations od1 ; od2 ; od3  challenge Einstein’s General
Relativity (GR). In recent years, many cosmologists have proposed various
modified gravitational models to explain these observations Buda ; Ferraro ;
mod1 ; rt11 ; rt12 ; f6 ; rt1 ; rt2 ; rt4 ; rt5 ; rt7 ; rt10 ; ff3 ; ff4 ; ff5
; rr5 ; rr6 . Although these models can be examined using numerical
techniques, analytical treatment is necessary to derive constraints on the
models’ free parameters and draw conclusions about their cosmological
viability.

In gravitational physics, the field equations are nonlinear differential
equations, making the derivation of analytic solutions a challenging task.
Certain approaches employed in the literature for constructing analytic
solutions rely on symmetry analysis ns1 ; ns2 ; ns3 ; ns4  and the Painlevé
algorithm p1 ; p2 ; p3 . However, in order to understand the global behaviour
of the physical properties in a given gravitational model, we can study the
behaviour of the field equations in the long run. The analysis of the
asymptotics has been widely used in gravitational theories, yielding many
interesting results dn1 ; dn2 ; dn3 ; dn4 ; dn5 ; dn6 ; dn7 ; dn8 ; dn9 ; dn10
.

In this work, we investigate the impact of nonlinear boundary corrections in
Symmetric Teleparallel General Relativity (STGR) nester  on cosmological
solutions. In particular we assume the gravitational Lagrangian to be
$L=Q+f\left(B\right)$, where $Q$ is the nonmetricity scalar, which is the
fundamental scalar of STGR, and $f\left(B\right)$ is an arbitrary function of
the boundary term $B$ which gives the difference of $Q$ with the Ricci scalar
of the background spacetime. This model belongs to the the family of
$f(Q,B)$-gravity, which proposed recently in ftc0 ; ftc1 ; ftc2  and it is a
generalization of $f\left(Q\right)$ theory lav2 ; lav3 . The boundary term has
been introduced before in the gravitational Action Integral previously within
the framework of teleparallelism bah  with various applications in
cosmological studies, see for instance ftb1 ; ftb2 ; ftb3 ; ftb4 ; ftb5 .

The boundary term introduces higher-order derivatives into the field
equations, which can be attributed to scalar fields anbound . These newly
introduced scalar fields play a role in the field equations, imparting dynamic
behavior to dark energy. In the following sections, we utilize asymptotic
analysis to establish constraints on the free parameters of the gravitational
theory and discuss the model’s viability based on initial conditions dn9 ;
dn10 . The paper is structured as follows.

In Section II we briefly discuss the basic definitions of STGR and we
introduce the boundary correction term. We consider a
Friedmann–Lemaître–Robertson–Walker (FLRW) geometry and in Section III we
present the field equations for the gravitational model of our consideration.
Section IV includes the main results of this study where we present a detailed
analysis of the dynamics for the cosmological field equations. Finally, in
Section V we summarize our results.

##  II Symmetric Teleparallel General Relativity

Consider a four-dimensional (non-Riemannian) manifold characterized by the
metric tensor $g_{\mu\nu}$ and the symmetric and flat connection
$\Gamma_{\mu\nu}^{\kappa}$, defining the covariant derivative $\nabla_{\mu}$.
Furthermore, we assume that the metric tensor and the connection
$\Gamma_{\mu\nu}^{\kappa}$ possess identical symmetries.

Because $\Gamma_{\mu\nu}^{\kappa}$ differs from the Levi-Civita connection it
follows $\nabla_{\lambda}g_{\mu\nu}\neq 0$. The tensor field

$Q_{\lambda\mu\nu}\equiv\nabla_{\lambda}g_{\mu\nu}=\frac{\partial
g_{\mu\nu}}{\partial
x^{\lambda}}-\Gamma_{\;\lambda\mu}^{\sigma}g_{\sigma\nu}-\Gamma_{\;\lambda\nu}^{\sigma}g_{\mu\sigma}$

is called the nonmetricity tensor and it is the essential for the STGR.

By definition, the connection is symmetric and flat, implying that both the
curvature tensor and the torsion tensor vanish.

From the nonmetricity tensor we can construct the scalar $Q$ as nester

$Q=Q_{\lambda\mu\nu}P^{\lambda\mu\nu},$ (1)

which is the Lagrangian function of STGR.

In particular, in STGR the Action Integral is defined as

$S_{STGR}=\int d^{4}x\sqrt{-g}Q.$ (2)

Tensor $\ P^{\lambda\mu\nu}$ is given by the following expression nester

$P_{\;\mu\nu}^{\lambda}=-\frac{1}{4}Q_{\;\mu\nu}^{\lambda}+\frac{1}{2}Q_{(\mu\phantom{\lambda}\nu)}^{\phantom{(\mu}\lambda\phantom{\nu)}}+\frac{1}{4}\left(Q^{\lambda}-\bar{Q}^{\lambda}\right)g_{\mu\nu}-\frac{1}{4}\delta_{\;(\mu}^{\lambda}Q_{\nu)},$
(3)

and

$Q_{\mu}=Q_{\mu\nu}^{\phantom{\mu\nu}\nu}~{},~{}\bar{Q}_{\mu}=Q_{\phantom{\nu}\mu\nu}^{\nu\phantom{\mu}\phantom{\mu}}.$
(4)

Let $\tilde{\Gamma}_{\mu\nu}^{\kappa}$ be the Levi-Civita connection for the
metric tensor $g_{\mu\nu}$, that is $\tilde{\nabla}_{\lambda}g_{\mu\nu}=0$,
and $\tilde{R}$ is the Ricciscalar defined by the Levi-Civita connection. Then
by definition lav3  $R-Q=B_{Q}$where $B$ is a boundary given by the expression
lav3  $B_{Q}=-\tilde{\nabla}_{\mu}\left(Q^{\mu}-\bar{Q}^{\mu}\right)$.

Consequently it follows

$\int d^{4}x\sqrt{-g}Q\simeq\int d^{4}x\sqrt{-g}R+\text{boundary terms.}$ (5)

and the theory is equivalent to GR.

###  II.1 Boundary corrections

The influence of the boundary term on gravitational phenomena has been
previously investigated, particularly in the context of teleparallel gravity
bah .

In this study we consider the modified gravitational Action Integral

$\hat{S}_{STGR}=\int d^{4}x\sqrt{-g}\left(Q+f\left(B\right)\right),$ (6)

where we introduce a nonlinear function $f$ which depend on the boundary $B$.
Action $\hat{S}_{STGR}$ belongs to the family of $f(Q,B)$ (also known as
$f\left(Q,C\right)$) theory ftc0 ; ftc1 ; ftc2 . In our consideration we
assume the dynamical degrees of freedom provided by the correction term
$f\left(B\right)$ to play the role of a dynamical dark energy.

The gravitational field equations are

$\displaystyle 0$
$\displaystyle=\frac{2}{\sqrt{-g}}\nabla_{\lambda}\left(\sqrt{-g}P_{\;\mu\nu}^{\lambda}\right)-\frac{1}{2}\left(Q+B\right)g_{\mu\nu}+\left(P_{\mu\rho\sigma}Q_{\nu}^{\;\rho\sigma}-2Q_{\rho\sigma\mu}P_{\phantom{\rho\sigma}\nu}^{\rho\sigma}\right)$
$\displaystyle+\left(\frac{B}{2}g_{\mu\nu}-\nabla_{\mu}\nabla_{\nu}+g_{\mu\nu}g^{\kappa\lambda}\nabla_{\kappa}\nabla_{\lambda}-2P_{~{}~{}\mu\nu}^{\lambda}\nabla_{\lambda}\right)f_{,B},$
(7)

or equivalent

$G_{\mu\nu}=T_{\mu\nu}^{f\left(B\right)},$ (8)

where now $G_{\mu\nu}$ is the Einstein tensor and
$T_{\mu\nu}^{f\left(B\right)}$ attributes the dynamical degrees of freedom of
the boundary term, that is,

$T_{\mu\nu}^{f\left(B\right)}=\left(\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}g^{\kappa\lambda}\nabla_{\kappa}\nabla_{\lambda}+2P_{~{}~{}\mu\nu}^{\lambda}\nabla_{\lambda}\right)f_{,B}+\frac{1}{2}\left(f\left(B\right)-Bf_{,B}\right)g_{\mu\nu}.$
(9)

We introduce the scalar field $\zeta=f_{,B}$; thus, the energy momentum tensor
$T_{\mu\nu}^{f\left(B\right)}$ becomes

$T_{\mu\nu}^{f\left(B\right)}=\left(\nabla_{\mu}\nabla_{\nu}\zeta-
g_{\mu\nu}g^{\kappa\lambda}\nabla_{\kappa}\nabla_{\lambda}\zeta+2P_{~{}~{}\mu\nu}^{\lambda}\nabla_{\lambda}\zeta\right)+\frac{1}{2}V\left(\zeta\right)g_{\mu\nu},$
(10)

where $V\left(\zeta\right)=\left(f\left(B\right)-Bf_{,B}\right).$

##  III FLRW Cosmology

On very large scales, the universe is isotropic and homogeneous, described by
the spatially flat FLRW geometry with the line element

$ds^{2}=-N\left(t\right)dt^{2}+a^{2}\left(t\right)\left(dx^{2}+dy^{2}+dz^{2}\right),$
(11)

where $a\left(t\right)$ is the scale factor and
$H\left(t\right)=\frac{1}{N}\frac{\dot{a}}{a}$ is the Hubble function and
$N\left(t\right)$ is the lapse function, where without loss of generality we
assume $N\left(t\right)=1$.

The FLRW geometry admits six isometries consisted by the three translation
symmetries

$\partial_{x}~{},~{}\partial_{y}~{},~{}\partial_{z}~{},$

and the three rotations

$y\partial_{x}-x\partial_{y}~{},~{}z\partial_{x}-x\partial_{z}~{},~{}z\partial_{y}-y\partial_{z}\text{.}$

The requirement for the connection $\Gamma_{\mu\nu}^{\kappa}$ to be symmetric,
flat, and to inherit the symmetries of the background geometry leads to three
distinct families of connections Heis2 ; Zhao . The cosmological field
equations for these three families of connections were derived previously in
anbound .

For the first connection, namely $\Gamma_{1}$, the modified Friedmann
equations are

$\displaystyle-3H^{2}$
$\displaystyle=3\dot{\zeta}H+\frac{1}{2}V\left(\zeta\right),$ (12)
$\displaystyle-2\dot{H}-3H^{2}$
$\displaystyle=\ddot{\zeta}+\frac{1}{2}V\left(\zeta\right),$ (13)

in which the scalar field $\zeta$ satisfy the Klein-Gordon equation

$\dot{H}+3H^{2}+\frac{1}{6}V_{,\zeta}=0.$ (14)

For connection $\Gamma_{2}$ the cosmological field equations are

$\displaystyle-3H^{2}$
$\displaystyle=3\dot{\zeta}H-\frac{3}{2}\dot{\zeta}\dot{\psi}+\frac{1}{2}V\left(\zeta\right),$
(15) $\displaystyle-2\dot{H}-3H^{2}$
$\displaystyle=\frac{3}{2}\dot{\zeta}\dot{\psi}+\ddot{\zeta}+\frac{1}{2}V\left(\zeta\right),$
(16)

where the scalar fields satisfy the equations motion

$\displaystyle\ddot{\zeta}+3H\dot{\zeta}$ $\displaystyle=0,$ (17)
$\displaystyle 6\dot{H}+18H^{2}-9H\dot{\psi}-3\ddot{\psi}+V_{,\zeta}$
$\displaystyle=0.$ (18)

Finally, for the third connection, i.e. $\Gamma_{3}$, the field equations are

$\displaystyle-3H^{2}$
$\displaystyle=3H\dot{\zeta}+\frac{3}{2}\frac{\dot{\zeta}}{a^{2}\dot{\Psi}}+\frac{1}{2}V\left(\zeta\right),$
(19) $\displaystyle-2\dot{H}-3H^{2}$
$\displaystyle=\frac{V\left(\zeta\right)}{2}+\frac{\dot{\zeta}}{8a^{2}\dot{\Psi}},$
(20)

and the equations of motion for the scalar fields are

$\displaystyle
3\dot{\Psi}\ddot{\zeta}+\dot{\zeta}\left(\dot{\Psi}H-2\ddot{\Psi}\right)$
$\displaystyle=0,$ (21) $\displaystyle
6\dot{H}+18H^{2}+V_{,\zeta}-\frac{3}{a^{2}\dot{\Psi}^{2}}\left(H\dot{\Psi}-\ddot{\Psi}\right)$
$\displaystyle=0.$ (22)

We emphasize that for connections $\Gamma_{2}$ and $\Gamma_{3}$, scalars
$\psi$ and $\Psi$ play crucial roles in the dynamics’ evolution. In contrast,
in the case of connection $\Gamma_{1}$, scalar $\gamma$ serves as a gauge
function. Connection $\Gamma_{1}$ is defined in the coincidence gauge, whereas
connections $\Gamma_{2}$ and $\Gamma_{3}$ are defined in the noncoincidence
gauge.

We continue our study with the phase-space analysis of the field equations
corresponding to the three connections. For the potential function
$V\left(\zeta\right)$ we consider the exponential function
$V\left(\zeta\right)=V_{0}e^{\lambda\zeta}$ which corresponds to the function

$f\left(B\right)=-\frac{B}{\lambda}\ln\left(-\frac{B}{\lambda
V_{0}}\right)-\frac{B}{\lambda}\text{.}$ (23)

##  IV Analysis of asymptotics

In the following, we conduct an analysis of the asymptotics for the considered
cosmological model. Specifically, we introduce dimensionless variables and
express the field equations as a set of algebraic-differential equations. We
calculate the stationary points and investigate their stability properties.
Each stationary point corresponds to an asymptotic solution with specific
physical properties. Finally, based on the stability properties, we can
establish constraints on the free parameters of the models and discuss the
initial value problem. This analysis is applied to the three different
families of connections.

###  IV.1 Connection $\Gamma_{1}$

To examine the asymptotic evolution of the field equations for the first
connection, we introduce the new variables

$~{}z=\frac{\dot{\zeta}}{H}~{},~{}y=\frac{V\left(\zeta\right)}{6H^{2}}~{},~{}\lambda=\frac{V_{,\zeta}}{V}~{},~{}\tau=\ln
a,$ (24)

with inverse transformation

$\dot{\zeta}=zH~{},~{}V\left(\zeta\right)=6yH^{2}~{},~{}V_{,\zeta}=\lambda
V~{},~{}a=e^{\tau}.$ (25)

Thus, the field equations transform into

$\displaystyle\frac{dz}{d\tau}$
$\displaystyle=3\left(1+z\right)+y\left(\lambda\left(2+z\right)-3\right),$
(26) $\displaystyle\frac{dy}{d\tau}$
$\displaystyle=y\left(6+\lambda\left(2y+z\right)\right),$ (27)
$\displaystyle\frac{d\lambda}{d\tau}$
$\displaystyle=\lambda^{2}z\left(\Gamma\left(\lambda\right)-1\right)~{},~{}\Gamma\left(\lambda\left(\zeta\right)\right)=\frac{V_{,\zeta\zeta}V}{\left(V_{,\zeta}\right)^{2}},$
(28)

and

$1+y+z=0.$ (29)

Moreover, the equation of state parameter $w_{eff}^{\Gamma_{1}}$ is expressed
as follows

$w_{eff}^{\Gamma_{1}}=1+\frac{2}{3}\lambda y.$ (30)

For the exponential potential, where $\lambda$ is always a constant and with
the application of the latter constraint equation we reduce the field
equations to the single differential equation

$\frac{dz}{d\tau}=\left(1+z\right)\left(6-\lambda\left(2+z\right)\right).$
(31)

The stationary points $A=\left(z\left(A\right)\right)$ of the latter equation
are two (in the finite and infinity regimes), point $A_{1}=-1$ and point
$A_{2}=\frac{6}{\lambda}-2$. Point $A_{1}$ corresponds to a stiff fluid
solution with $w_{eff}^{\Gamma_{1}}\left(A_{1}\right)=1$, while for point
$A_{2}$ we calculate
$w_{eff}^{\Gamma_{1}}\left(A_{2}\right)=-3+\frac{2\lambda}{3}$, where it
follows that the de Sitter universe is recovered for $\lambda=3$, and the
solution describes acceleration for $\lambda<4$. Finally, for $\lambda>6$
point $A_{1}$ is an attractor while for $\lambda<6$, the attractor is point
$A_{2}$.

###  IV.2 Connection $\Gamma_{2}$

We introduce the dimensionless variables

$x=\frac{\dot{\psi}}{2H}~{},~{}z=\frac{\dot{\zeta}}{H}~{},~{}y=\frac{V\left(\zeta\right)}{6H^{2}}~{},~{}\lambda=\frac{V_{,\zeta}}{V}~{},~{}\tau=\ln
a,$ (32)

that is,

$\dot{\psi}=2xH~{},~{}\dot{\zeta}=zH~{},~{}V\left(\zeta\right)=6yH^{2}~{},~{}V_{,\zeta}=\lambda
V~{},~{}a=e^{\tau}.$ (33)

In terms of the new variables the field equations are

$\displaystyle\frac{dx}{d\tau}$
$\displaystyle=\frac{1}{2}\left(3x\left(y-2z-1\right)+3x^{2}z+3\left(1+z\right)+\left(2\lambda-3\right)y\right),$
(34) $\displaystyle\frac{dz}{d\tau}$
$\displaystyle=\frac{3}{2}z\left(y-1+z\left(x-1\right)\right),$ (35)
$\displaystyle\frac{dy}{d\tau}$
$\displaystyle=y\left(3\left(1+y\right)+\left(\lambda-3\left(1-x\right)\right)z\right),$
(36) $\displaystyle\frac{d\lambda}{d\tau}$
$\displaystyle=\lambda^{2}z\left(\Gamma\left(\lambda\right)-1\right)~{},~{}\Gamma\left(\lambda\left(\zeta\right)\right)=\frac{V_{,\zeta\zeta}V}{\left(V_{,\zeta}\right)^{2}}.$
(37)

Friedmann’s first equation yields the constraint

$1+z\left(1-x\right)+y=0,$ (38)

while the equation of state parameter reads

$w_{eff}^{\Gamma_{2}}=y-z\left(1-x\right).$ (39)

For the exponential potential, i.e. $\lambda$ is a constant, the stationary
points $B=\left(B\left(x\right),B\left(z\right),B\left(y\right)\right)$ are

$B_{1}=\left(x,-\frac{1}{1-x},0\right)~{},~{}B_{2}=\left(1-\frac{\lambda}{3},0,-1\right).$
(40)

$B_{1}$ describes a family of points with correspond to stiff fluid solutions,
that is, $w_{eff}^{\Gamma_{2}}\left(B_{1}\right)=1$. On the other hand,
$B_{2}$ describes the de Sitter universe with
$w_{eff}^{\Gamma_{2}}\left(B_{2}\right)=-1$.

As far as the stability is concerned the eigenvalues of the two-dimensional
system in the space of variables $\left\\{x,z\right\\}$, around the stationary
points $B_{1}$ are $\left\\{0,6-\frac{\lambda}{1-x}\right\\}$, while around
the point $B_{2}$ are $\left\\{-3,-3\right\\}$. Thus, point $B_{2}$ is always
an attractor, while for the stability properties of $B_{1}$ we should employ
the center manifold theorem (CMT).

We introduce the new variable $\bar{z}=z+\frac{1}{1-x}$, such that the
coordinates of points $B_{1}$ to be $B_{1}=\left(x,0,0\right)$. Then, we
assume $\bar{z}=h\left(x\right)$ in order to determine the center manifold. In
order a stable manifold to exist it should hold $h\left(x_{1}\right)=0$ and
$\frac{dh}{dx}|_{x\rightarrow x_{1}}=0$. We calculate
$h\left(x\right)=\frac{1}{1-x}-h_{0}\left(1-x-\lambda\right)$; thus, points
$B_{1}$ describe always unstable solutions.

####  IV.2.1 Poincare variables

Because the dynamical variables are not constraint, they can take values at
the infinity. Hence, in order to study the analysis at the infinity we
introduce the Poincare variables

$x=\frac{X}{\sqrt{1-X^{2}-Z^{2}}}~{},~{}z=\frac{Z}{\sqrt{1-X^{2}-Z^{2}}}~{},~{}dT=\sqrt{1-X^{2}-Z^{2}}d\tau,$

where $\left\\{X^{2},Z^{2}\right\\}\leq 1$.

In terms of the new variables the field equations are expressed as

$\frac{dX}{dT}=F_{1}\left(X,Z\right)~{},~{}\frac{dZ}{dT}=F_{2}\left(X,Z\right),$
(41)

while the equation of state parameter is

$w_{eff}^{\Gamma_{2}}=-1-2Z\frac{X-\sqrt{1-X^{2}-Z^{2}}}{1-X^{2}-Z^{2}}.$ (42)

The stationary points
$B^{\infty}=\left(B^{\infty}\left(X\right),B^{\infty}\left(Z\right)\right)$ at
the infinity are

$B_{1\pm}^{\infty}=\left(\pm 1,0\right)~{},~{}B_{2\pm}^{\infty}=\left(0,\pm
1\right),$ (43)

and for $\lambda=3$, there exist the family of points

$B_{3\pm}^{\infty}=\left(X,\pm\sqrt{1-X^{2}}\right).$ (44)

Stationary points $B_{1\pm}^{\infty}$ describe de Sitter solutions, i.e.
$w_{eff}^{\left(\Gamma_{2}\right)}\left(B_{1\pm}^{\infty}\right)=-1$, while
the asymptotic solutions at points $B_{2\pm}^{\infty}$ describe Big Crunch or
Big Rip singularities, that is
$w_{eff}^{\left(\Gamma_{2}\right)}\left(B_{2\pm}^{\infty}\right)=\mp\infty$.
Similarly, the family of points $B_{3\pm}^{\infty}$ describe Big Crunch and
Big Rip singularities.

As far as the stability is concerned, the eigenvalues for the linearized
system around points $B_{1\pm}^{\infty}$ are $\left\\{0,0\right\\}$, from
where we infer that the stationary points describe unstable solutions. For
points $B_{2\pm}^{\infty}$ the eigenvalues are $\left\\{\pm
6,\pm\left(\lambda-3\right)\right\\}$, which means that the Big Crunch
solution $B_{2-}^{\infty}$ is an attractor for $\lambda>3$. Finally, for
$\lambda=3$ the stability of the points $B_{3\pm}^{\infty}$ depend on the sing
of the dynamical variable $X$.

In Fig. 1 we present phase-space portraits for the dynamical system, where it
is clear for $\lambda<3$, the unique attractor is the de Sitter solution
described by point $B_{2}$. Moreover, in Fig. 2 we present qualitative
evolution of the equation of state parameter for various sets of initial
conditions.

Figure 1: Phase-space portraits for the cosmological field equations of
connection $\Gamma_{2}$ in the Poincare variables (41). The phase-space
portraits are for $\lambda=1$, $\lambda=3$ and $\lambda=4$. With dots are the
stationary points and red lines correspond to the family of points $B_{1}$. We
observe that for $\lambda<3$, the unique attractor is the de Sitter solution
described by point $B_{2}$.  Figure 2: Qualitative evolution of the equation
of state parameter $w_{eff}^{\left(\Gamma_{2}\right)}$ for different values of
$\lambda=\left(1,3,4\right)$ and for various initial conditions
$\left(X_{0},Z_{0}\right)$. Solid lines are for $\left(0.95,0.1\right)$,
dashed lines are for $\left(0.9,0.19\right)$, dotted lines are for
$\left(0.8,0.3\right)$ and dash-dotted lines are for $\left(-0.8,0.3\right)$.

###  IV.3 Connection $\Gamma_{3}$

For the set of field equations related to the connection $\Gamma_{3}$ we
introduce the dimensionless variables

$\bar{x}=\frac{1}{2a^{2}H\dot{\Psi}}~{},~{}z=\frac{\dot{\zeta}}{H}~{},~{}y=\frac{V\left(\zeta\right)}{6H^{2}}~{},~{}\lambda=\frac{V_{,\zeta}}{V}~{},~{}\tau=\ln
a,$ (45)

that is,

$\frac{1}{\dot{\Psi}}=2a^{2}\bar{x}H~{},~{}\dot{\zeta}=zH~{},~{}V\left(\zeta\right)=6yH^{2}~{},~{}V_{,\zeta}=\lambda
V~{},~{}a=e^{\tau}.$ (46)

The field equations read

$\displaystyle\frac{d\bar{x}}{d\tau}$
$\displaystyle=\frac{\bar{x}}{2z}\left(z\left(1+\bar{x}\right)-3+y\left(3-2\lambda\left(1-z\right)\right)\right),$
(47) $\displaystyle\frac{dz}{d\tau}$
$\displaystyle=3+\left(3-\bar{x}\right)z+y\left(\lambda\left(2+z\right)-3\right),$
(48) $\displaystyle\frac{dy}{d\tau}$
$\displaystyle=y\left(6+\lambda\left(2y+z\right)\right),$ (49)
$\displaystyle\frac{d\lambda}{d\tau}$
$\displaystyle=\lambda^{2}z\left(\Gamma\left(\lambda\right)-1\right)~{},~{}\Gamma\left(\lambda\left(\zeta\right)\right)=\frac{V_{,\zeta\zeta}V}{\left(V_{,\zeta}\right)^{2}},$
(50)

and constraint

$1+y+\left(1+\bar{x}\right)z=0.$ (51)

Moreover, the equation of state parameter is expressed as

$w_{eff}^{\Gamma_{3}}=1+\frac{2}{3}\lambda y.$ (52)

We remark that for the exponential potential $\lambda$ is always constant.
Hence the stationary points
$C=\left(C\left(\bar{x}\right),C\left(z\right),C\left(y\right)\right)$ for the
latter algebraic-differential system are

$C_{1}=\left(0,-1,0\right)~{},~{}C_{2}=\left(0,\frac{2}{\lambda}\left(3-\lambda\right),1-\frac{6}{\lambda}\right).$
(53)

Point $C_{1}$ describes a stiff fluid solution, with
$w_{eff}^{\Gamma_{3}}\left(C_{1}\right)=1$. On the other hand, the asymptotic
solution at point $C_{2}$ describes an ideal gas with equation of state
parameter $w_{eff}^{\Gamma_{3}}\left(C_{2}\right)=-3+\frac{2}{3}\lambda$. The
stationary point describes acceleration for $\lambda<4$, while the
cosmological constant is recovered for $\lambda=3$.

We make use of the constraint equation(51) and we reduce by one the dimension
of the dynamical system. The eigenvalues of the linearized system around the
stationary points $C_{1}$ and $C_{2}$ are $\left\\{2,6-\lambda\right\\}$;
$\left\\{\lambda-6,\frac{1}{2}\left(3\lambda-14\right)\right\\}$ respectively.
Therefore, point $C_{1}$ is a saddle point when $\lambda>6$, otherwise is a
source; while point $C_{2}$ is an attractor for $\lambda<\frac{14}{3}$. We
remark that when $C_{2}$ describes acceleration it is always attractor.

####  IV.3.1 Poincare variables

For the analysis at the infinity regime we work in the two-dimensional space
defined by the dynamical variables $\left\\{\bar{x},z\right\\}$.

The Poincare variables are defined as

$\bar{x}=\frac{\bar{X}}{\sqrt{1-\bar{X}^{2}-Z^{2}}}~{},~{}z=\frac{Z}{\sqrt{1-\bar{X}^{2}-Z^{2}}}~{},~{}dT=\sqrt{1-\bar{X}^{2}-Z^{2}}d\tau,$

where $\left\\{\bar{X}^{2},Z^{2}\right\\}\leq 1$.

The field equations are reduced to the system of the form

$\frac{dX}{dT}=F_{1}\left(X,Z\right)~{},~{}\frac{dZ}{dT}=F_{2}\left(X,Z\right),$
(54)

and the equation of state parameter is expressed as

$w_{eff}^{\Gamma_{3}}=1-\frac{2}{3}\lambda\left(1+\frac{Z\left(\bar{X}+\sqrt{1-X^{2}-Z^{2}}\right)}{1-X^{2}-Z^{2}}\right).$
(55)

The stationary points
$C^{\infty}=\left(X\left(C^{\infty}\right),Z\left(C^{\infty}\right)\right)$ at
the infinity; that is,
$1-\left(X\left(C^{\infty}\right)\right)^{2}-\left(Z\left(C^{\infty}\right)\right)^{2}=0$,
are

$C_{1\pm}^{\infty}=\left(0,\pm 1\right).$

We calculate
$w_{eff}^{\Gamma_{3}}\left(C_{1\pm}^{\infty}\right)=\mp\lambda\infty$. Hence,
for $\lambda>0$, $C_{1+}^{\infty}$ corresponds to a Big Rip singularity, and
$C_{1-}^{\infty}$ to a Big Crunch; nevertheless for $\lambda<0$,
$C_{1+}^{\infty}$ corresponds to a Big Crunch singularity, and
$C_{1-}^{\infty}$ to a Big Rip singularity.

The eigenvalues of the linearized system around the stationary points are
$\left\\{\pm 2\lambda,0\right\\}$. Because the second eigenvalue is zero, we
employ the CMT and we found that the stationary points does not posses any
submanifold where the solutions are stable. Thus, the stationary points are
saddle points or sources.

In Fig. 3 we present phase-space portraits for the dynamical system in
Poincare variables. Furthermore, in Fig. 2 we present qualitative evolution of
the equation of state parameter for various sets of initial conditions.

Figure 3: Phase-space portraits for the cosmological field equations of
connection $\Gamma_{3}$ in the Poincare variables (41). The phase-space
portraits are for $\lambda=2$, $\lambda=3$ and $\lambda=4$. With dots are the
stationary points. We observe that the unique attractor is point $C_{2}$,
while the two stationary points at the infinity regime are saddle points.
Figure 4: Qualitative evolution of the equation of state parameter
$w_{eff}^{\left(\Gamma_{3}\right)}$ for different values of
$\lambda=\left(2,2.5,3\right)$ and for various initial conditions
$\left(X_{0},Z_{0}\right)$. Solid lines are for $\left(-0.95,0.05\right)$,
dashed lines are for $\left(-0.8,0.2\right)$, dotted lines are for
$\left(0.8,-0.2\right)$ and dash-dotted lines are for $\left(0.5,0.7\right)$.

##  V Conclusions

We conducted a detailed analysis of the asymptotic dynamics for an extension
of STGR in which nonlinear components of the boundary term are introduced in
the gravitational integral. In STGR, the definition of the connection is not
unique, and for the spatially flat FLRW, there are three families of different
connections. Although the selection of the connection does not affect the
gravitational model in STGR when nonlinear terms of the boundary scalar are
introduced, new dynamical degrees of freedom appear. The new degrees of
freedom can be attributed to scalar fields, leading to three different sets of
gravitational field equations, corresponding to the families of connections.

For the three different models, we employed dimensionless variables and
determined the stationary points. We reconstructed the asymptotic solutions at
the stationary points and investigated their stability properties. The results
are summarized in Table 1.

The three gravitational models admit asymptotic solutions that can describe
late-time acceleration. However, for connections $\Gamma_{2}$ and
$\Gamma_{3}$, asymptotic solutions appear that can describe Big Rip or Big
Crunch singularities. For connection $\Gamma_{2}$, a Big Crunch singularity is
stable when the parameter $\lambda>3$. Thus, in this case, for the unique
attractor to be the de Sitter solution, it follows that $\lambda\leq 3$.
However, for these two specific models, the phase-space analysis leads to the
derivation of sets for the initial conditions to avoid the appearance of such
types of singularities.

The $Q+f\left(B\right)$ theory provides a mechanism to explain the late-time
acceleration phase of the universe. Nevertheless, it is not obvious from this
analysis if this theory can describe other eras of the cosmological history
and if it solves the $H_{0}$-tension. In a future study, we plan to
investigate this specific problem in the context of
$f\left(Q,B\right)$-gravity.

Table 1: Asymptotic solutions in STGR with boundary correctons Point | $\mathbf{w}_{eff}$ | Acceleration | Stability
---|---|---|---
Connection $\mathbf{\Gamma}_{1}$
$A_{1}$ | $1$ | No | Stable$~{}\lambda>6$
$A_{2}$ | $-3+\frac{2}{3}\lambda$ | $\lambda<4$ | Stable$~{}\lambda<6$
Connection $\mathbf{\Gamma}_{2}$
$B_{1}$ | $1$ | No | Unstable
$B_{2}$ | $-1$ | Yes | Stable
$B_{1\pm}^{\infty}$ | $-1$ | No | Unstable
$B_{2\pm}^{\infty}$ | $\mp\infty$ | Big Rip $B_{2+}^{\infty}$ |  $B_{2-}^{\infty}~{}$Stable $\lambda>3$
$B_{3\pm}^{\infty}~{}\left(\lambda=3\right)$ | $\mp\infty$ | Yes | Stable
Connection $\mathbf{\Gamma}_{3}$
$C_{1}$ | $1$ | No | Unstable
$C_{2}$ | $-3+\frac{2}{3}\lambda$ | $\lambda<4$ | Stable $\lambda<\frac{14}{3}$
$C_{1\pm}^{\infty}$ | $\mp\lambda\infty$ | Big Rip | Unstable

###### Acknowledgements.

AP thanks the support of VRIDT through Resolución VRIDT No. 096/2022 and
Resolución VRIDT No. 098/2022.

## References

  * (1) A.G. Riess et al., Astron. J. 116, 1009 (1998)
  * (2) M. Tegmark et al., Astrophys. J. 606, 702 (2004)
  * (3) E. Komatsu et al., Astrophys. J. Suppl. Ser. 180, 330 (2009)
  * (4) S. Nojiri, S.D. Odintsov and V.K. Oikonomou, Phys. Rept. 692, 1 (2017)
  * (5) T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82 451, (2010)
  * (6) S. Bahamonde, K. F. Dialektopoulos, C. Escamilla-Rivera, G. Farrugia, V. Gakis, M. Hohmann, J. L. Said, J. Mifsud and E. Di Valentino, Rep. Prog. Phys. 86, 026901 (2023)
  * (7) H.A. Buchdahl, Mon. Not. Roy. Astron. Soc. 150, 1 (1970)
  * (8) R. Ferraro and F. Fiorini, Phys. Rev. D 75, 084031 (2007)
  * (9) J. B. Jimenez, L. Heisenberg and T. Koivisto, Phys. Rev. D 98, 044048 (2018)
  * (10) S. M. Carrol, V. Duvvuri, M. Trodden and M. S. Turner, Phys. Rev. D., 70, 043528 (2004);
  * (11) L. Amendola, D. Polarski and S. Tsujikawa, Phys. Rev. Lett. 98, 131302 (2007)
  * (12) T. Clifton and J.D. Barrow, Class. Quant. Grav. 23, 2951 (2006)
  * (13) R.C. Nunes, A. Bonilla, S. Pan and E.N. Saridakis, EPJC 77, 230 (2017)
  * (14) A. de la Cruz-Dombriz, A. Dobado and A. L. Maroto, Phys. Rev. D 80, 124011 (2009)
  * (15) J. Aftergood and A. DeBenedictis, Phys. Rev. D 90, 124006 (2014)
  * (16) A. Lymperis, JCAP 11, 018 (2022)
  * (17) S.H. Shekh, Phys. Dark Univ. 33, 100850 (2021)
  * (18) F.K. Anagnastopoulos, S. Basilakos and E.N. Saridakis, Phys. Lett. B 822, 136634 (2021)
  * (19) W. Hu and I. Sawicki, Phys. Rev. D. 76, 064004 (2007)
  * (20) M. Krssak, R. J. van den Hoogen, J. G. Pereira, C. G. Boehmer and A. A. Coley, Class. Quantum Grav. 36, 183001 (2019)
  * (21) M. Tsamparlis and A. Paliathanasis, Symmetry 10, 233 (2018)
  * (22) K. Dialektopoulos, J.L. Said, Z. Oikonomopoulou, Eur. Phys. J. C 82, 259 (2022)
  * (23) A.R. Akbarieh, P.S. Ilkhchi and Y. Kucukakca, Int. J. Geom. Meth. Mod. Phys. 20, 2350106 (2023)
  * (24) S. Hembron, R. Bhaumik, S. Dutta and S. Chakraborty, Eur. Phys. J. C 84, 110 (2024)
  * (25) S. Cotsakis and P.G.L. Leach, J. Phys. A: Math. Gen. 27, 1625 (1994)
  * (26) J. Demaret and C. Scheen, J. Math. Phys. A: Math. Gen. 29, 59 (1996)
  * (27) A. Paliathanasis, J.D. Barrow and P.G.L. Leach, Phys. Rev. D 94, 023525 (2016)
  * (28) C.R Fadragas and G. Leon, Class. Quantum Grav. 31, 195011 (2014)
  * (29) R. Lazkoz and G. Leon, Phys. Lett. B 638, 303 (2006)
  * (30) T. Gonzales, G. Leon and I. Quiros, Class. Quantum Grav. 23, 3165 (2006)
  * (31) A.A. Coley, Phys. Rev. D 62, 023517 (2000)
  * (32) H. Farajollahi and A. Salehi, JCAP 07, 036 (2011)
  * (33) E.J. Copeland, A.R. Liddle and D. Wands, Phys. Rev. D 57, 4686 (1998)
  * (34) P. Christodoulidis, D. Roest and E.I. Sfakianakis, Scaling attractors in multi-field inflation, JCAP 12, 059 (2019)
  * (35) L. Amendola, D. Polarski and S. Tsujikawa, IJMPD 16, 1555 (2007)
  * (36) L. Amendola, D. Polarski and S. Tsujikawa, Phys. Rev. Lett. 98, 131302 (2007)
  * (37) L. Amendola, R. Gannouji, D. Polarski and S. Tsujikawa, Phys. Rev. D 75, 083504 (2007)
  * (38) J.M. Nester and H-J Yo, Chinese Journal of Physics 37, 113 (1999)
  * (39) C.G. Boehmer and E. Jensko, Phys. Rev. D 104, 024010 (2021)
  * (40) S. Capozziello, V. De Falco and C. Ferrara, Eur. Phys. J. C. 83, 915 (2023)
  * (41) A. De, T.-H. Loo and E.N. Saridakis, Non-metricity with bounday terms: f(Q,C) gravity and cosmology, (2023) [arXiv:2308.00652]
  * (42) L. Heisenberg, Phys. Reports 796, 1 (2019) lav2,lav3
  * (43) J.B. Jimenez, L. Heisenberg, T. Koivisto and S. Pekar, Phys. Rev. D 101, 103507 (2020)
  * (44) S. Bahamonde, C.G. Bohmer and M. Wright, Phys. Rev. D 10, 104042 (2015)
  * (45) M. Wright, Phys. Rev. D 93, 103002 (2016)
  * (46) M. Caruana, G. Farrugia and J.L. Said, EPJC 80, 640 (2020)
  * (47) S. Bahamonde, A. Golovnev, M.-J. Guzmán, J.L. Said and C. Pfeifer, JCAP 01, 037 (2022)
  * (48) A. Paliathanasis, JCAP 1708, 027 (2017)
  * (49) A. Paliathanasis and G. Leon, EPJC 81, 718 (2021)
  * (50) A. Paliathanasis, Phys. Dark Univ. 43, 101388 (2024)
  * (51) F. D’ Ambrosio, L. Heisenberg and S. Kuhn, Class. Quantum Grav. 39 025013 (2022)
  * (52) D. Zhao, Eur. Phys. J. C 82, 303 (2022)