\section{Introduction}
This paper examines whether using graph embeddings\cite{Goyal2018} with a dimensionality reducing algorithm, t-SNE\cite{vandermaaten2208}, to visualise argumentation graphs and their structural properties through unsupervised learning can lead to interesting results. This approach has been effective at data visualisation in Deep Neural Network research for graph-based methods\cite{Pezzotti2019}. The contribution of this paper is to show that this technique also holds promise in the visualisation of argumentation graphs by applying it at both the node and the graph level using both a standard and a custom built embedding approach. In particular, we show that it is possible to clearly visualise the functional partitions of arguments in the Sembuster domain\cite{Rodrigues2017} using this method and to separate argumentation domains into visual clusters at the graph level using a custom GCN embedding, raising the possibility that both differences between argument graphs and the function of arguments within argumentation graphs can be clustered and shown in a visually intuitive way using unsupervised methods.

\section{Graph Embeddings}
Graph Embeddings are used in a range of graph analysis applications including node classification, link prediction, clustering, and visualisation either directly or as additional input features to machine learning algorithms. There are several different algorithmic approaches to generating graph embeddings that are suited to different use cases. Goyal and Ferrara define three main approaches in their 2018 survey of the field\cite{Goyal2018}.

The first category is factorisation-based approaches that share the property of working with a matrix representation (e.g. adjacency or Laplacian matrix) and a proximity measure to calculate the node embeddings. The next category is based on random walks through the graph to generate the embedding. Deepwalk, for instance, maximises the probability of seeing the last \textit{k} and the next \textit{k} nodes in the random walk centred at a given node. The last major category of approaches is based on Deep Neural Networks. For instance, Graph Convolutional Networks (GCN) generate an embedding by iteratively aggregating neighbourhood embeddings\cite{Kipf2016}. In this paper, we will use HOPE\cite{Ou2016}, a factorisation based approach, and a GCN based neural network approach to obtain the embeddings that we will visualise for argumentation graphs.

\section{Visualising Embeddings using t-SNE}
Embeddings are usually of low dimensionality, but even so, they are not easy to visualise. However, once an embedding has been generated for complex data types such as graphs, words, or images, it becomes possible to use dimensionality reduction techniques to visualise them effectively.

In the Deep Neural Network community, the dimensionality reduction technique of choice is t-distributed stochastic neighbor embedding (t-SNE). In contrast to other common dimensionality reduction techniques such as Principal Component Analysis (PCA), t-SNE does not rely on a linear projection but uses local relationships between points. It uses Student's t-distribution to model the relationship between points in the higher dimensional space and then recreates those relationships in the lower dimensional space by way of a gradient descent based algorithm\cite{vandermaaten2208}.

A simple example that demonstrates the way this visualisation approach works can be found by applying it to the MNIST dataset of handwritten digits\cite{lecun-mnisthandwrittendigit-2010}\footnotetext[1]{Reproduced from \cite{Pezzotti2019}}. We see that t-SNE 2-D embeddings are able to cluster the 10 digits from MNIST dataset into 10 distinct clusters \cite{Pezzotti2019}. This shows that graph embeddings can extract visual similarity from the raw data without supervision.

\begin{figure}[h!]
	\centering
	\begin{subfigure}[b]{0.3\linewidth}
		\includegraphics[width=\linewidth]{mnist_tsne.png}
		\caption{MNIST Visualised with t-SNE.}
	\end{subfigure}
	\begin{subfigure}[b]{0.3\linewidth}
		\includegraphics[width=\linewidth]{cora_tsne.png}
		\caption{CORA Visualised with t-SNE.}
	\end{subfigure}
	\par
	\begin{subfigure}[b]{0.3\linewidth}
		\includegraphics[width=\linewidth]{word2vec_tsne.png}
		\caption{Word2Vec Visualised with t-SNE.}
	\end{subfigure}
	\caption{Examples of t-SNE Visualisation}
	\label{fig:mnist_cora2}
\end{figure}

A more practical example can be found when using t-SNE with word embeddings. Word embeddings project words into an embedding space in an analogous way to graph embeddings. A good word embedding would, therefore, place words with similar meanings closer to each other in the embedding space. Using Word2Vec\cite{Mikolov2013}, a word embedding trained by trying to predict individual words given a context, the following chart of cities was generated. In general, cities from the same countries are placed closer together, although there are some anomalies\footnotetext[2]{Reproduced from https://nlpforhackers.io/word-embeddings/}.

This approach also works well for graph-structured data. Once a graph embedding has been created for a graph it can also be visualised using t-SNE and the results will also tend to preserve a visual representation of the proximity measure that the embedding optimises. An example of this can be found by applying the GraphSAGE mode\cite{Hamilton2017}, a Deep Learning approach, to the CORA dataset of academic papers and citations\footnotetext[3]{Reproduced from https://towardsdatascience.com/using-graphsage-to-learn-paper-embeddings-in-cora-a94bb1e9dc9d}. This results in clusters on the class of the papers in the dataset, represented by the point colours, but with several exceptions that invite further investigation. In this case, we see that the information in the graph embeddings enables a clustering close to the underlying distribution of document classes without prior knowledge.

\section{Visualising Argumentation Graphs}
\subsection{Introduction}
Argumentation graphs are to some extent themselves a way of visualising arguments. However, for large argumentation structures with hundreds or thousands of arguments looking at the unadorned graph is of little utility. In the following sections, we will give two examples of how one can use graph embeddings to visualise information about the structure of arguments, one using node-level information to analyse a single argument and one to visualise properties of whole argumentation graphs. Both examples are based on the Abstract Argumentation formalism\cite{Dung1995}, but there is nothing inherent in this approach that limits it to one type of graph-based representation.

\subsection{Node-Level Visualisation}
To examine the possibilities of using Graph Embeddings and t-SNE to visualise argumentation graphs, we have started with a formal domain, Sembuster, that has the interesting property of having three types of arguments with distinct structural characteristics. The Sembuster domain, originally proposed by Caminada and Verheij\cite{Caminada2010}, is composed of unique graphs generated for each cardinality, $k$, that can be partitioned into three different types: A, B, and C.
\begin{wrapfigure}{r}{0.4\linewidth}
	\begin{center}
		\includegraphics[width=\linewidth]{sembuster.png}
	\end{center}
	\caption{Sembuster Scheme}
\end{wrapfigure}
Arguments in the A partition only attack themselves. An argument $B_i$ attacks all arguments $A_j$ where $i  >= j$, arguments $B_j$ where $i  > j$ and the argument $C_i$. An argument $C_i$ attacks the corresponding argument $B_i$.

To visualise this configuration, we trained a graph embedding using the HOPE\cite{Ou2016} method on a set of Sembuster graphs with a cardinality from 300 to 4500. HOPE is a directionality preserving embedding that works well on directed graphs such as argumentation graphs. The graphs were trained with a dimensionality of 128 features using a single K80 GPU. Training time was from 48 seconds to 5 minutes 12 seconds for a single embedding. These graphs were then visualised using 2-dimensional t-SNE using a colour coding corresponding to the argument partitions and subjected to visual inspection.
\begin{figure}[h!]
	\centering
	\begin{subfigure}[b]{0.3\linewidth}
		\includegraphics[width=\linewidth]{sembuster_600_vis.pdf}
		\caption{k=600}
	\end{subfigure}
	\begin{subfigure}[b]{0.3\linewidth}
		\includegraphics[width=\linewidth]{sembuster_1200_vis.pdf}
		\caption{k=1200}
	\end{subfigure}
	\begin{subfigure}[b]{0.3\linewidth}
		\includegraphics[width=\linewidth]{sembuster_1800_vis.pdf}
		\caption{k=1800}
	\end{subfigure}
	\caption{T-SNE Visualisation of Sembuster Graphs. Cyan represents partition A, Red partition B, and Blue partition C.}
	\label{fig:mnist_cora}
\end{figure}
If we examine the three examples, the three types of arguments stand out clearly. The self-referential A arguments are curled in on themselves in the visualisation and the related B arguments are closer to the A than the unrelated C arguments.
The string-like construction of the Sembuster graphs seems at least superficially to also be represented in the string-like nature of the visualisation. While this is a highly formal example, it does indicate that the functional structure of arguments can in some cases be visualised using this type of technique, although its application to arguments of a less formal nature would need further investigation. Applying similar techniques on other argumentation graphs that are under analysis would potentially be able to show functional clusterings of arguments based on the type of graph embedding that has been applied. It may even be possible to design specific embedding approaches that allow training based on the particular properties we are interested in analysing for.

\subsection{Graph-Level Visualisation}
We also trained a GCN model\cite{Kipf2016} to classify variants of graphs from different abstract argumentation graph domains. We trained the embeddings on 10 separate domains using a custom built model. The domains are defined in Rodrigues et al.\cite{Rodrigues2017} and the table below give a description of their basis. The dataset was based on a subset of the tasks from the 2nd International Competition on Computational Models of Argumentation using an even sample drawn from these 10 domains, which all form part of the competition corpus.

The GCN model used 4 convolutional layers followed by 2 fully connected layers to try to predict the class of a given input graph. The implication would be that the visualisation should show the graphs in a given domain following in a recognisable pattern (e.g. clustering or equal spacing). This could for instance be useful in situations where you have an argumentation corpus that is not clustered or classified and you want to group them by similarity and be able to show that similarity in a visual mode. After the model had been trained for 4 hours, the training was stopped and the output of the last convolutional layer was used as a graph embedding for visualisation by extracting the raw features of the last convolutional layer.

\begin{table}[H]
	\caption{Domains for Graph-Level Visualisation}
	\begin{tabularx}{\textwidth}{|l|l|X|}
		\textbf{Identifier} & \textbf{Domain}  & \textbf{Description}\\
		afinput & ABA2AF & Assumption-Based Argumentation translated to abstract argumentation frameworks \\
		admbuster & AdmBuster & AdmBuster graphs, based on Caminada and Podlaszewski\cite{Caminada2017} \\
		BA & Barabasi-Albert & Barabasi-Albert graphs, randomly generated \\
		ER & Erdös-Rényi & Erdös-Rényi graphs, randomly generated  \\
		grd & GroundedGenerator & Randomly generated argumentation frameworks containing only a grounded extension \\
		.cnf & Planning2AF & Planning problems transformed to abstract argumentation problems \\
		sembuster & SemBuster & SemBuster graphs, see section 4.2\\
		scc & SccGenerator & Randomly generated argumentation frameworks containing multiplestrongly connected components  \\
		.gml. & Traffic  & Traffic networks converted to abstract argumentation frameworks\\
		WS & Watts-Strogatz & Watts-Strogatz graphs, randomly generated\\
	\end{tabularx}
	\label{tab:addlabel}%
\end{table}

\begin{figure}[h!]
	\centering
	\begin{subfigure}[b]{0.45\linewidth}
		\includegraphics[width=\linewidth]{graph_embedding.pdf}
	\end{subfigure}
	\begin{subfigure}[b]{0.45\linewidth}
		\includegraphics[width=\linewidth]{sem_scc.pdf}
	\end{subfigure}
	\begin{subfigure}[b]{0.45\linewidth}
		\includegraphics[width=\linewidth]{afp_adm_ba.pdf}
	\end{subfigure}
	\begin{subfigure}[b]{0.45\linewidth}
		\includegraphics[width=\linewidth]{gml_WS.pdf}
	\end{subfigure}

	\caption{Graph-level Visualisation using t-SNE.}
	\label{fig:graphlevel}
\end{figure}
This had more ambiguous results. For most domains it was possible to identify a clear separation between classes either by clustering or by equal spacing. Examples of these are shown in Figure~\ref{fig:graphlevel}. First, for the \textit{ER}, \textit{.cnf}, and \textit{grd} domains we can see a clear separation between the three classes through three spacing patterns that are distinct by class. For \textit{sembuster} and \textit{scc} this pattern is even more distinct. However, for the afinput, admbuster, and BA domains it is less clear with less separation between the classes. For .gml and WS there is some overlap showing that there is a separation between the classes, but that some graphs are difficult to distinguish.

While this approach does demonstrate the feasibility of training and visualising graph-level embeddings of argumentation graphs, it is some way from being practically useful in its present form. What it does demonstrate is that a properly trained unsupervised embedding will be able to accurately separate argumentation graphs based on a training task that will then be available for visualisation. That means potentially being able to cluster corpora of argumentation graphs based on a variety of training tasks in order to discover new ways of classifying and categorising the arguments. In this case a simple supervised training task is used to generate a separation into already known classes, but both other well-known graph embeddings or a custom designed embedding for the task at hand might give improved results.

\section{Conclusion}
In this paper, we have shown that it is possible to use graph embeddings combined with t-SNE to visualise properties of argumentation graphs at both the node and the graph level. The node-level experiments were more successful, showing three distinct clusters according to the function of the arguments in the graph. However, the graph level example did shown that visual clustering of argument graphs by using a graph embedding is a least possible, although it needs some refinement to be applicable in practice. Further development of both the training task and the network architecture should yield improved results.

Wider exploration of this space would require the analysis of argumentation graphs of diverging provenance across many different domains. This would also involve testing a wider range of graph embeddings and eventually exploring the creation of embeddings from scratch with the kind of properties that would make them particularly useful for visualising argumentation graphs. The aim would be to enable the clustering argumentation graphs in an unsupervised manner to enable visual analysis of similarity at the graph level and commensurately to be able to see visual clustering of the individual arguments in graphs to be able to visually group arguments within a graph by some measure of similarity generated by the graph embedding.

\section{Introduction}
This document provides instructions for style and layout and how to submit the final
version. Although it was written for individual authors contributing to IOS Press books,
it can also be used by the author/editor preparing a monograph or an edited volume.

Authors should realize that the manuscript submitted by the volume editor to IOS
Press will be almost identical to the final, published version that appears in the book,
except for the pagination and the insertion of running headlines. Proofreading as
regards technical content and English usage is the responsibility of the author.

A template file for \LaTeX2e is available from http://www.iospress.nl/service/authors/\allowbreak latex-and-word-tools-for-book-authors/. \LaTeX{} styles required for the \LaTeX{} template are also
available.\footnote{For authors using MS Word separate Instructions as well
as a Word template are available from http://www.iospress.nl/service/authors/latex-and-word-tools-for-book-authors/.}

\section{Typographical Style and Layout}

\subsection{Type Area}
The text output area is automatically set within an area 12.4 cm
horizontally and 20 cm vertically. Please do not use any
\LaTeX{} or \TeX{} commands that affect the layout or formatting of
your document (i.e. commands like \verb|\textheight|,
\verb|\textwidth|, etc.).

\subsection{Font}

The font type for running text (body text) is 10~point Times New Roman.
There is no need to code normal type (roman text). For literal text, please use
\texttt{type\-writer} (\verb|\texttt{}|)
or \textsf{sans serif} (\verb|\textsf{}|). \emph{Italic} (\verb|\emph{}|)
or \textbf{boldface} (\verb|\textbf{}|) should be used for emphasis.

\subsection{General Layout}
Use single (1.0) line spacing throughout the document. For the main
body of the paper use the commands of the standard \LaTeX{}
``article'' class. You can add packages or declare new \LaTeX{}
functions if and only if there is no conflict between your packages
and the \texttt{IOS-Book-Article}.

Always give a \verb|\label| where possible and use \verb|\ref| for cross-referencing.

\subsection{(Sub-)Section Headings}
Use the standard \LaTeX{} commands for headings: {\small \verb|\section|, \verb|\subsection|, \verb|\subsubsection|, \verb|\paragraph|}.
Headings will be automatically numbered.

Use initial capitals in the headings, except for articles (a, an, the), coordinate
conjunctions (and, or, nor), and prepositions, unless they appear at the beginning
of the heading.

\subsection{Footnotes and Endnotes}
Please keep footnotes to a minimum. If they take up more space than roughly 10\% of
the type area, list them as endnotes, before the References. Footnotes and endnotes
should both be numbered in arabic numerals and, in the case of endnotes, preceded by
the heading ``Endnotes''.

\subsection{References}

Please use the Vancouver citing \& reference system, and the National Library of
Medicine (NLM) style.

Place citations as numbers in square brackets in the text. All publications cited in
the text should be presented in a list of references at the end of the manuscript.
List the references in the order in which they appear in the text. Some examples of
the NLM style:

\medskip
\noindent\ul{Journal article:}\par\noindent
Petitti DB, Crooks VC, Buckwalter JG, Chiu V. Blood pressure levels before dementia.
Arch Neurol. 2005 Jan;62(1):112-6.

\medskip
\noindent\ul{Paper from a proceedings:}\par\noindent
Rice AS, Farquhar-Smith WP, Bridges D, Brooks JW. Canabinoids and pain. In: Dostorovsky
JO, Carr DB, Koltzenburg M, editors. Proceedings of the 10th World Congress on Pain;
2002 Aug 17-22; San Diego, CA. Seattle (WA): IASP Press; c2003. p. 437-68.

\medskip
\noindent\ul{Contributed chapter in a book:}\par\noindent
Whiteside TL, Heberman RB. Effectors of immunity and rationale for immunotherapy. In:
Kufe DW, Pollock RE, Weichselbaum RR, Bast RC Jr, Gansler TS, Holland JF, Frei~E~3rd,
editors. Cancer medicine 6. Hamilton (ON): BC Decker Inc; 2003. p. 221-8.

\medskip
\noindent\ul{Book by author(s):}\par\noindent
Jenkins PF. Making sense of the chest x-ray: a hands-on guide. New York: Oxford
University Press; 2005. 194 p.

\medskip
\noindent\ul{Edited book:}\par\noindent
Izzo JL Jr, Black HR, editors. Hypertension primer: the essentials of high blood pressure.
3rd ed. Philadelphia: Lippincott Williams \& Wilkins; c2003. 532 p.

\medskip
\noindent\ul{Proceedings:}\par\noindent
Ferreira de Oliveira MJ, editor. Accessibility and quality of health services. Proceedings of
the 28th Meeting of the European Working Group on Operational Research Applied to Health
Services (ORAHS); 2002 Jul 28-Aug 2; Rio de Janeiro, Brazil. Frankfurt (Germany): Peter Lang;
c2004. 287 p.

\section{Illustrations}

\subsection{General Remarks on Illustrations}
The text should include references to all illustrations. Refer to illustrations in the
text as Table~1, Table~2, Figure~1, Figure~2, etc., not with the section or chapter number
included, e.g. Table 3.2, Figure 4.3, etc. Do not use the words ``below'' or ``above''
referring to the tables, figures, etc.

Do not collect illustrations at the back of your article, but incorporate them in the
text. Position tables and figures with at least 2 lines
extra space between them and the running text.

Illustrations should be centered on the page, except for small figures that can fit
side by side inside the type area. Tables and figures should not have text wrapped
alongside.

Place figure captions \textit{below} the figure, table captions \textit{above} the table.
Use bold for table/figure labels and numbers, e.g.: \textbf{Table 1.}, \textbf{Figure 2.},
and roman for the text of the caption. Keep table and figure captions justified. Center
short figure captions only.

The minimum \textit{font size} for characters in tables is 8 points, and for lettering in other
illustrations, 6 points.

On maps and other figures where a \textit{scale} is needed, use bar scales rather than
numerical ones of the type 1:10,000.

\subsection{Quality of Illustrations}
Embed the fonts used if the application provides that option.
Ensure consistency by using similar sizes and fonts for a group of small figures.
To add lettering to figures, it is best to use Helvetica or Arial (sans serif fonts)
to avoid effects such as shading, outline letters, etc.

 Avoid using illustrations
taken from the Web. The resolution of images intended for viewing on a screen is
not sufficient for the printed version of the book.

If you are incorporating screen captures, keep in mind that the text may not be
legible after reproduction.

\subsection{Color Illustrations}
Please note, that illustrations will only be printed in color if the volume editor agrees to
pay the production costs for color printing. Color in illustrations will be retained
in the online (ebook) edition.

\section{Equations}
Number equations consecutively, not section-wise. Place the numbers in parentheses at
the right-hand margin, level with the last line of the equation. Refer to equations in the
text as Eq. (1), Eqs. (3) and (5).

\section{Fine Tuning}

\subsection{Type Area}
\textbf{Check once more that all the text and illustrations are inside the type area and
that the type area is used to the maximum.} You may of course end a page with one
or more blank lines to avoid `widow' headings, or at the end of a chapter.

\subsection{Capitalization}
Use initial capitals in the title and headings, except for articles (a, an, the), coordinate
conjunctions (and, or, nor), and prepositions, unless they appear at the beginning of the
title or heading.

\subsection{Page Numbers and Running Headlines}
You do not need to include page numbers or running headlines. These elements will be
added by the publisher.

\section{Submitting the Manuscript}
Submit the following to the volume editor:

\begin{enumerate}
\item The main source file, and any other required files. Do not submit more than
one version of any item.

\item The source files should compile without errors with pdflatex or latex.

\item Figures should be submitted in EPS, PDF, PNG or JPG format.

\item A high resolution PDF file generated from the source files you submit.
\end{enumerate}


\section{\label{sec:Introduction}Introduction}
The origin of dark matter (DM) and dark energy (DE) has been one of the main pursuits for new physics for many decades. DM was proposed as an invisible mass that holds a galaxy together while DE was an explanation of the accelerating universe. Although there are many candidates from the particle physics side such as Axion-Like Particle (ALP) which could act as both DM and DE depending on their masses \cite{Arvanitaki:2009fg,Marsh:2015xka}, but none of these have been experimentally confirmed yet.
On the other hand, a natural place for which the candidate for DM and DE could rise is the gravity theory. Over the past decades, a tremendous progress has been made in generalisation of Einstein gravity especially in the class of models called the Scalar-Tensor theory \cite{Fujii:2003pa, Faraoni:2004pi} and reference therein.
One of the most studied and celebrated models is Brans-Dicke (BD) theory \cite{Brans:1961sx}. Motivated by Mach's principle, the scalar field is included to modify the gravity in a way that Weak Equivalence Principle is still intact. As a result, the gravitational constant is promoted to a scalar field which coupled to both mass and geometry. A huge number of studies for this particular model can be found here \cite{Kolitch:1994kr,Santos:1996jc,Copeland:1997et,Abdalla:2007qi,deSouza:2005ig,Hrycyna:2013hla,Hrycyna:2013yia,Hrycyna:2014cka,Garcia-Salcedo:2015naa,Papagiannopoulos:2016dqw,Felegary:2016znh,Roy:2017mnz,Ghaffarnejad_2017,Lu:2019twd,Shabani_2019,Zucca:2019ohv,Giacomini:2020grc}.

The idea of extra dimensions has also been an essential development for the unification of fundamental forces \cite{Bailin:1987jd}. The scalar fields from extra-dimensional theories such String/M theory are ubiquitous and play a crucial role in several phenomena such as the inflaton field \cite{Gasperini:1993hu,Damour:1995pd,Chamblin:1999ya}, DM \cite{Damour:1990tw}, dark radiation\cite{Svrcek:2006yi,Acharya:2015zfk} and cosmological constant \cite{Wetterich:1987fm,Shaposhnikov:2008xi}. In our previous work \cite{Pongkitivanichkul:2020txi}, we have proposed an extension of Einstein gravity using degrees of freedom coming from compactification of extra-dimensions, see \cite{Overduin:1998pn} for review and references therein. In particular, the toroidal 5-dimensional spacetime ($\mathbb{R}^4 \times S^1$) is the starting point and the 5-dimensional metric containing a scalar field and a vector field is assumed. This process of circular compactification is also known as Kaluza-Klein (KK) compactification. We also include 2 additional gauge fields in order to solve the anisotropy problem of the energy momentum tensor. In our previous paper, dynamical system analysis has been carried out and we demonstrated the potential to have both DM and DE phases of the universe.

It was shown that including a matter field into the system can sometimes changes the behaviour of the dynamical system \cite{Hrycyna:2013hla, Hrycyna:2013yia}. In this paper, we would like to provide a more extensive study of our model. In particular, we are interested in the inclusion of matter field and radiation into the existing dynamical system in order to achieve the semi-realistic universe. The paper is organised as follows. First, we give a review on the model in section \ref{sec:The Model} where the additional matter and radiation field are also included. The full list of equation of motions will be given in section \ref{sec:EOM}. In section \ref{sec:Lower Mode Cases}, the range of compactification radius is chosen such that the number of scalar fields relevant to low energy physics will be finite and the numerical analysis can be done. The dynamical system analysis is performed in section \ref{sec:Dynamical System}. The full results will be given in section \ref{sec:Result} and the discussion of the results can be found in section \ref{sec:Discussion}. Finally, the conclusion is given in section \ref{sec:Conclusion and Outlook}.
\section{\label{sec:Conclusion and Outlook}Conclusion and Outlook}

In this paper, we have revisited the KK inspired BD model where barotropic fluid is included. The UV limit of this model is the traditional 5-dimensional KK action with 1 additional massive scalar field and 2 gauge fields. The KK compactification reduces the UV theory to a 4-dimensional gravity theory with a dilaton coupling to a tower of scalar fields. Another relic of the higher dimensions is the gauge field from 5-dimensional metric which forms a Cosmic Triad solution with 2 additional gauge fields. Then the usual Einstein field equations with the barotropic fluid are calculated. Together with equations of motion, the complete autonomous system is constructed and the appropriate dimensionless parameters are defined.

With the dynamical system approach, we found that realistic DM and DE phases are readily present in many solutions. Since our universe has undergone the matter dominated phase and currently it is in the DE dominated phase, the DM solutions are expected to be saddle points with $w_{\rm eff} = 0$. In the analysis we found 1 DM solution, i.e,  \hyperlink{FB1}{$\mathcal{F}_1^{(B)}$} (and equivalently \hyperlink{FD1}{$\mathcal{F}_1^{(D)}$}). The \hyperlink{FB1}{$\mathcal{F}_1^{(B)}$} critical point consists of the oscillating $\dot{\eta}_1$ field and ($\phi$, $A^{\mu}$, $\rho_m$) which play a role of background fields. We found 4 dark energy solutions which are stable critical points with $w_{\rm eff} = -1$, i.e.,  \hyperlink{FBbar3}{$\bar{\mathcal{F}}_1^{(B)}$}, \hyperlink{FCbar1}{$\bar{\mathcal{F}}_1^{(C)}$}, \hyperlink{FCbar3}{$\bar{\mathcal{F}}_3^{(C)}$} and \hyperlink{FDbar1}{$\bar{\mathcal{F}}_1^{(D)}$}. Noted that all of these are in the slow-roll approximation as one might expected. Interestingly, we also found many solutions which resemble the phantom DE with $w_{\rm eff} < -1$. These are \hyperlink{FB2}{$\mathcal{F}_2^{(B)}$}, \hyperlink{FDbar5}{$\bar{\mathcal{F}}_5^{(D)}$}, and \hyperlink{FDbar4}{$\bar{\mathcal{F}}_4^{(D)}$}. All of these solutions contain non-zero kinetic terms which suggests that the phantom DE could change the value along the history of the universe.

Although this approach could not provide the origin of the barotropic fluid at the 5-dimensional level, we consider this approach as an effective analysis of the overall contributions of other fields. The lack of a strong connection to UV physics is compensated by the fact that the revisited model presented here has many more facets of the DM/DE behaviours comparing to the original model. One could conclude that the barotropic fluid in the KK inspired BD model plays an important role in enriching phenomena of the model tremendously.

\section{\label{sec:EOM}Equation of Motions}
In this section, we compute the equations of motion of the KK inspired BD model and we will employ the results to investigate the dynamics of the universe in this model as a candidate of DM and DE. With the flat FLRW metric, $ds^2 = dt^2 - a^2(t)(dx^2+dy^2+dz^2)$, we firstly obtain the equation of motion for the dilaton field, $\phi$ from the Euler-Lagrange equation as
\begin{align}
       \mathcal{R} =&~ \frac{4}{3\sqrt{-g}\phi} \partial^{\mu}(\sqrt{-g}\partial_{\mu}\phi) - \frac{3}{4}\phi^2 F_{\mu\nu} F^{\mu\nu} - \frac{2}{3\phi^2}\partial^{\mu}\phi\partial_{\mu}\phi \nonumber\\
        & + \sum_{n=0}^{\infty}\left[ \partial_{\mu}\eta_n \partial^{\mu}\eta_n - M_{(5)}^2 \eta^2 + \left(\frac{1}{\phi^2} + A_{\nu}A^{\nu}\right)\frac{n^2}{R_k^2}\eta_n^2 \right] \nonumber\\
        & +  \sum_{a=1}^{2}\sum_{n=0}^{\infty} \left( -\frac{3}{4}\phi^2\widetilde{F}^a_{n\mu\nu}\widetilde{F}^{a\mu\nu}_n + \left( V'' - \frac{3}{4}\frac{\phi^2 n^2}{ R_k^2} \right) \widetilde{A}^a_{\mu,n} \widetilde{A}^{\mu,a}_n \right)
        .
\end{align}
According to the action in Eq.(\ref{full-action}), the matter Lagrangian together with the barotropic fluid can defined as
\begin{eqnarray}
\label{eq:Lmatter}
     \mathcal{L}_{matter} &=& -\frac{1}{4}\phi^{3}F^{2}-\frac{2}{3}\frac{\partial_{\mu}\phi\partial^{\mu}\phi}{\phi}
     \nonumber\\
     && +\, \phi \left[ \partial_{\mu}\eta_n \partial^{\mu}\eta_n - M_{(5)}^2 \eta^2 + \left(\frac{1}{\phi^2} + A_{\nu}A^{\nu}\right)\frac{n^2}{R_k^2}\eta_n^2 \right]
          + \mathcal{L}_{fluid}
        ,
\end{eqnarray} where $\mathcal{L}_{fluid}$ represents other fields (particles) which could be existing in the 4 dimensions.
The energy-momentum tensor can be determined by
\begin{align}
\label{eq:T}
    T^{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}\mathcal{L}_{matter})}{\delta g^{\mu\nu}}.
\end{align}
Substituting Eq.(\ref{eq:Lmatter}) into Eq.(\ref{eq:T}), we find
\begin{align}
     T_{\mu\nu}=&~\phi^{3}F_{\mu\alpha}F_{\nu}\!^{\alpha}-\frac{g_{\mu\nu}}{4}\phi^3 F^2 - \frac{g_{\mu\nu}}{4}\phi^3 \sum_{a=1}^2 \Tilde{F}_{\mu\nu,a}\Tilde{F}^{\mu\nu}_a + \frac{2}{3}\frac{\partial_{\mu}\phi\partial_{\nu}\phi}{\phi} - \phi \sum_{n=0}^{\infty} \partial_{\mu}\eta_{n}\partial_{\nu}\eta_{n} \nonumber\\
        & - 2 \phi A_{\mu}A_{\nu}\sum_{n=0}^{\infty} \frac{n^2\eta_n^2}{R_k^2} + g_{\mu\nu} \phi \sum_{n=0}^{\infty}\left[\left(\frac{1}{\phi^2} +  A_{\nu}A^{\nu}\right)\frac{n^2}{R_k^2}\eta_n^2 - M_{(5)}^2 \eta_n^2\right] \nonumber \\
        & + \phi^3 \sum_{n=0}^{\infty} \widetilde{F}^a_{\mu\alpha,n} \widetilde{F}^{a,n}_{\nu} - \phi \sum_{n=0}^{\infty} \left( V'' - \frac{n^2 \phi^2}{4 R_k^2} \right) (2 \widetilde{A}^a_{\mu,n} \widetilde{A}^{a}_{\nu,n} - g_{\mu\nu} \widetilde{A}^a_{\rho,n} \widetilde{A}^{\rho,a}_{n}) + \widetilde{T}_{\mu\nu}\,.
        \label{energy-momentum-tensor}
\end{align}
The barotropic fluid energy-momentum, $\widetilde{T}_{\mu\nu}$ is given by
\begin{eqnarray}
\widetilde{T}_{\mu\nu} = (\rho_m + \rho_r)u_\mu u_\nu + p_r g_{\mu\nu}\,,
\end{eqnarray}
where $u_\mu$ is the usual comoving 4-velocity, whereas $\rho_m$, $\rho_r$ and $p_r$ are the energy density of dust matter, radiation and pressure of radiation, respectively. In order to fulfill the isotropy of the energy momentum tensor, the zero mode configuration of each species $(a =1,2)$ needs to satisfy
\begin{align}
    A^{\mu} =& (0,A,0,0)\\
   \widetilde{A}^{\mu,1}_{0} \equiv \widetilde{A}^{\mu,1}  =& (0,0,A,0)\\
   \widetilde{A}^{\mu,2}_{0} \equiv \widetilde{A}^{\mu,2} =& (0,0,0,A).
\end{align}
This specific configuration is also known as the Cosmic Triad \cite{ArmendarizPicon:2004pm,Golovnev:2008cf,Maleknejad:2012fw}.
In addition, the diagonal spatial components of the energy-momentum tensor in Eq.(\ref{energy-momentum-tensor})
require that
\begin{equation}
   V'' = \sum_{n=0}^{\infty} \frac{n^2\eta_n^2}{R_k^2} \equiv M_A^2.
\end{equation}
The Einstein field equation of the KK inspired BD gravity takes the form of
\begin{equation}
     \phi(\mathcal{R}_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\mathcal{R}) + g_{\mu\nu}\nabla_{\sigma}\nabla^{\sigma}\phi - \nabla_{\mu}\nabla_{\nu}\phi = T_{\mu\nu}.
     \label{EFE}
\end{equation}
Considering the $tt$ component, the Einstein tensor reads
\begin{equation}
G_{tt}=\frac{\dot{a}^2}{a^2},
\end{equation}
with the fact that
\begin{align}
    \widetilde{F}^{t\sigma,1} \widetilde{F}^{t,1}\!_{\sigma} &= \widetilde{F}^{t\sigma,2}\widetilde{F}^{t,2}\!_{\sigma} = F^{t\sigma}F^{t}\!_{\sigma} = g_{ij}F^{ti}F^{tj} = -a^2\left(\dot{A} + 2 H A\right)^2 \\
    \widetilde{F}^{1}_{\mu\nu} \widetilde{F}^{\mu\nu,1} &= \widetilde{F}^{2}_{\mu\nu} \widetilde{F}^{\mu\nu,2} =  F_{\mu\nu}F^{\mu\nu} = F_{ti}F^{ti} + F_{it}F^{it} = -2a^2\left(\dot{A} + 2 H A\right)^2\,.
\end{align}
Then the Friedmann equation in the KK inspired BD model can be written as
\begin{eqnarray}
        \left(\frac{\dot{a}}{a}\right)^{2} &=& -\frac{3\phi^2a^2}{2}\left(\dot{A} + 2 H A\right)^2 +\frac{2\dot{\phi}^2}{3\phi^2}
        \nonumber\\
        && -\, \sum_{n=0}^{\infty}\dot{\eta}_n^2  + \sum_{n=0}^{\infty} \left[ \left(\frac{1}{\phi^2} -3 a^2 A^2\right)\frac{n^2}{R_k^2} - M_{(5)}^2 \right]\eta_n^2 + \frac{\rho_m + \rho_r}{\phi} .
\end{eqnarray}
Recalling relation $\dot{H}=\frac{\Ddot{a}}{a}-\frac{\dot{a}^2}{a^2}$, the Reychaudhuri equation can be determined from the spatial components of the field equation in Eq.(\ref{EFE}) and it reads
\begin{align}
    - 3\ddot{\phi} =& ~ 3\phi H^2 + 6\phi \left(\dot{H} + H^2\right) + 9 H \dot{\phi} \nonumber \\
     -& \frac{3 a^2\phi^3}{2} \left(\dot{A} + 2 H A\right)^2 - 6 \phi a^2 A^2 \sum_{n=0}^{\infty} \frac{n^2\eta_n^2}{R_k^2} \nonumber\\
    -& 3\phi \sum_{n=0}^{\infty}\left[\left(\frac{1}{\phi^2} -3 a^2 A^2\right)\frac{n^2}{R_k^2}\eta_n^2 - M_{(5)}^2 \eta_n^2\right] + a^2\textcolor{black}{p_r}. \label{eq:EFEij}
\end{align}
In this section, we have prepared relevant equations of motion in the KK inspired BD gravity model to study the dynamics of the universe in the later. Since, the scalar field $\eta_n$ represents the infinite summation of the KK exited states and leads to complicated physical quantities in the model. Therefore we will consider only the lower mode of the KK excited states in the next section.
\section{Discussion}\label{sec:Discussion}
From the previous section, we obtained a number of critical points from KK inspied BD gravity model with barotropic matter. The critical points are classified by their stability. Then the effective equation of state for each critical point are calculated. In analysis, the stability and effective equation of state are used to interpret the physical meaning behind each critical point. Consequently, we acquired critical points which are both the compatible and incompatible with the observable universe.

For the DM, there are 2 critical points which responsible for it. The first one comes from $\mathcal{F}_i^{(B)}$ system where non-vanishing fields are $\phi, \eta_1, A^\mu, \rho_m$, i.e., \hyperlink{FB1}{$\mathcal{F}_1^{(B)}$}, while the other one comes from the $\mathcal{F}_i^{(D)}$ system, i.e., \hyperlink{FD1}{$\mathcal{F}_1^{(D)}$}, where all fields are not vanished. Apparently, these two critical points are the same critical point coming from different system. Since the effective equation of state of these critical points is zero, i.e. $w_\text{eff}=0$, and their stability are saddle, these critical points can be represented as matter-dominated phase in the standard cosmological model. Noted that the matter-dominated phase is contributed by both matter and DM, then DM behaviour is obtainable from our model.

For the DE, we have 4 critical points from 3 systems, namely, $\bar{\mathcal{F}}_i^{(B)}$, $\bar{\mathcal{F}}_i^{(C)}$ and $\bar{\mathcal{F}}_i^{(D)}$. The first point is \hyperlink{FBbar3}{$\bar{\mathcal{F}}_1^{(B)}$} where $\phi, \eta_1, A^\mu, \rho_m$ are non-vanishing fields. The second and third points are \hyperlink{FCbar1}{$\bar{\mathcal{F}}_1^{(C)}$} and \hyperlink{FCbar3}{$\bar{\mathcal{F}}_3^{(C)}$} where  $\phi, \eta_0, \eta_1, A^\mu, \rho_m$ are non-vanishing fields. The fourth point is \hyperlink{FDbar1}{$\bar{\mathcal{F}}_1^{(D)}$} where all fields are not vanished. Noted that for barred system we assumed that the fields $\eta_0$ and $\eta_1$ are slow-roll fields causing the dynamics of these fields to be zero, i.e. $\dot{\eta}_0=\dot{\eta}_1=0$. These 4 critical points are all stable points and their effective equations of state are simply $-1$, $w_\text{eff}=-1$ which give us accelerating universe driven by DE. This means that our model is sufficient to predict DE phase.

Moreover, we also obtain the phantom DE solutions, i.e. $w_\text{eff}<-1$. This solution gives us the highly accelerated expanding universe leading to the Big Rip. We have 3 suitable critical points for the phantom DE. The first one is \hyperlink{FB2}{$\mathcal{F}_2^{(B)}$} where $\phi, \eta_1, A, \rho_m$ are non-vanishing fields. This critical point will be stable point with effective equation of state, $w_\text{eff}<-5/3$ for $\lambda > 3/\mu$. The second one is \hyperlink{FDbar5}{$\bar{\mathcal{F}}_5^{(D)}$} where non-vanishing fields are the same as \hyperlink{FB2}{$\mathcal{F}_2^{(B)}$} excepting that $\eta_1$ is a slow-roll field in this system, i.e. $\dot{\eta}_1=0$. The critical point is stable point with effective equation of state, $w_\text{eff}=-5/3$. The third one is \hyperlink{FDbar4}{$\bar{\mathcal{F}}_4^{(D)}$} where all fields are not vanished. For this critical point, we assumed that $\dot{\eta}_0$ and $\dot{\eta}_1$ are slow-roll fields. This critical point appears to be a stable point and its effective equation of state is less than $-5/3$, i.e., $w_\text{eff}=-5/3$.

For the latter, we also obtain critical points which incompatible with the observable universe. Now, we are getting into details of those critical points. Starting from \hyperlink{FBbar4}{$\bar{\mathcal{F}}_4^{(B)}$}, \hyperlink{FCbar2}{$\bar{\mathcal{F}}_2^{(C)}$} and \hyperlink{FDbar3}{$\bar{\mathcal{F}}_3^{(D)}$}, they are all saddle point. The effective equations of state for these critical points are perfectly stand on the border of decelerated and accelerated expansion of the universe, i.e., $w_\text{eff}=-1/3$. This suggests constantly expanded universe and we will call this situation as critical DE phase. Now, for \hyperlink{FA1}{$\mathcal{F}_1^{(A)}$}, \hyperlink{FA2}{$\mathcal{F}_2^{(A)}$} and \hyperlink{FC1}{$\mathcal{F}_1^{(C)}$}, they are stable points and their effective equation of state, $w_\text{eff}$, is zero. This kind of critical point represents the universe evolving to matter-dominated phase and remaining at this phase forever.
Next, \hyperlink{FBbar2}{$\bar{\mathcal{F}}_2^{(B)}$}, \hyperlink{FBbar3}{$\bar{\mathcal{F}}_3^{(B)}$} and \hyperlink{FDbar2}{$\bar{\mathcal{F}}_2^{(D)}$}, they are saddle points with DE’s equation of state, i.e. $w_\text{eff}=-1$. The universe, according to the properties of these critical points, is evolving to DE dominated era. Then, it will transit to unpredictable phase. Lastly \hyperlink{FB3}{$\mathcal{F}_3^{(B)}$} and \hyperlink{FD2}{$\mathcal{F}_2^{(D)}$}, they are also saddle points with effective equation of state related to phantom DE, i.e. $w_\text{eff}=-5/3$. This predicts that the phantom DE is not the final phase of the universe. The universe will continually evolve to the unknown phase after.

\section{\label{sec:The Model}A Kaluza-Klein inspired Brans-Dicke Model}
In this section, we set up the crucial ingredients of the KK inspired BD model and use for studying the dynamical system analysis in the latter. First of all, the 5-dimensional action for KK gravity with free massive scalar field, $\Tilde{\eta}$ is given by \cite{Pongkitivanichkul:2020txi}
\begin{align}
    S=\int d^4 x dy \sqrt{-\Tilde{g}}\sum_{a=1}^2 \left( -\Tilde{\mathcal{R}}+\hat{\partial}_A \Tilde{\eta}^*\partial^A \Tilde{\eta}-M^2_{(5)} \Tilde{\eta}^* \Tilde{\eta} - \frac{1}{4}\phi^2\Tilde{F}^a_{AB}\Tilde{F}^{aAB} + V(\Tilde{A}^a_M) \right)
    ,
\end{align}
where $M_{(5)}$ is the mass of the scalar field $\tilde\eta$ and $\tilde{\mathcal{R}}$, $\phi$, $\tilde{F}_{AB}^a \equiv \partial_A\,A^a_B - \partial_B\,A^a_A$, $A_M^a$ are Ricci scalar, dilaton field, strength tensor and gauge field in 5 dimensions, respectively. In addition, we use natural unit system where $\sqrt{16\pi G} = 1$. The capital Latin indices represent the bulk dimensional spacetime indices, $A,B,C,\cdots = 0,1,2,3,5$ whereas the Latin indices stand for the triad gauge field indices, $a,b,c,\cdots = 1,2$. While $V( \tilde{A}_M^a)$ is the potential and it is composed of the the gauge field $\tilde{A}_M^a$ only and it will be specified in the later. Moreover, all variables with tilde symbol, $\tilde{~}$, are described the physical quantities in 5-dimensional spacetime. The metric tensor $\tilde g_{AB}$ of the bulk 5-dimensional spacetime is written by \cite{Overduin:1998pn},
\begin{eqnarray}
\widetilde{g}_{AB} = \left(
\begin{array}{cc}
g_{\mu\nu} + \phi^2 A_\mu A_\nu & \phi^2 A_\mu
\\
\phi^2 A_\nu & \phi^2
\end{array}
\right).
\end{eqnarray}
It has been shown in Ref. \cite{Pongkitivanichkul:2020txi} that the gravitational action of the KK theory can be reduced to the 4-dimensional spacetime by mean of the compactification of the $y$ dimension. The full form of the KK theory in the 4 dimensions is written in the following form,
\begin{eqnarray}
S = S_{BD} + S_{KK} + S_G\,,
\label{full-action}
\end{eqnarray}
where $S_{BD}$, $S_{KK}$ and $S_G$ are the actions of the BD gravity, the 4 dimensional KK scalar field and the gauge fields, respectively. The explicit forms of three actions above are given by,
\begin{eqnarray}
S_{BD} &=& \int d^4 x \sqrt{-g}\,\phi\left[ -\mathcal{R} - \frac{1}{4}\phi^2F_{\mu\nu}F^{\mu\nu} - \frac{2}{3}\frac{\partial^{\mu}\phi \partial_{\mu}\phi}{\phi^2}\right],
\label{BD-action}
\\
S_{KK} &=& \int d^4 x \sqrt{-g}\,\phi \sum_{n=0}^{\infty}\left[ \partial_{\mu}\eta_n \partial^{\mu}\eta_n - M_{(5)}^2 \eta^2 + \left(\frac{1}{\phi^2} + A_{\nu}A^{\nu}\right)\frac{n^2}{R_k^2}\eta_n^2 \right],
\label{4d-scalar-KK}
\\
S_G &=& \int d^4x \sqrt{-g}\,\phi \sum_{a=1}^{2}\sum_{n=0}^{\infty} \Bigg[ -\frac{1}{4}\phi^2\widetilde{F}^a_{n\mu\nu}\widetilde{F}^{a\mu\nu}_n + \left( V'' - \frac{\phi^2 n^2}{4 R_k^2} \right) \widetilde{A}^a_{\mu,n} \widetilde{A}^{\mu,a}_n
\nonumber
\\
&&\qquad\qquad\qquad\qquad\quad\quad   -\frac{1}{4}\phi^2 \partial_{\mu} \widetilde{A}_{5,n}^a \partial^{\mu}\widetilde{A}^{5,a}_n + V''_{(5)}\widetilde{A}^a_{5,n} \widetilde{A}^{5,a}_n \Bigg],
\label{gauge-action}
\end{eqnarray}
where $V'' = \frac{\delta^2 V}{\delta A_{\mu}^2}$ and $V''_{(5)} = \frac{\delta^2 V}{\delta A_{5}^2}$. To obtain the actions in Eqs.(\ref{BD-action},\ref{4d-scalar-KK},\ref{gauge-action}), we have assumed that the $\eta_n$ and $\phi$ are real scalar fields. For all detail derivations of the actions, we refer to Ref.\cite{Pongkitivanichkul:2020txi} and we will not repeat them in this work. We close this section by summarising the main result in this model. We have shown that the BD gravity with scalar and gauge fields can be generated from the dimensional compactification of the 5-dimensional KK theory. The dilaton field is coupled to both Ricci scalar and the scalar field in 4-dimensional actions. In the following sections, we will use the KK inspired BD model in the presence of the barotropic fluid matter to demonstrate that the dilaton, scalar and gauge fields might play the role of DM and DE.

\section{\label{sec:Result}Results}
Since there are 9 main dimensionless parameters, $X_1-X_9$ and 1 extra parameter, $\mu$. It is more convenient to categorise the dynamical system into many cases. First, let's consider the case that the 5-dimensional mass, $M_{(5)}$ does not exist in the universe. Therefore, the $X_8$ is decoupled from the Friedmann equation. We analyse with assumption, $X_3=X_8=0$, and found that there is no real solution in this case.
Another case is the assumption that $M_{(5)}\neq0$. We will consider this in various sets of the non-vanishing fields in the dynamical system. All possible real solutions are summarised in Tab.(\ref{tab:tab1}). Moreover, the slow-roll approximation of the $\dot\eta_{0}$ and $\dot\eta_{1}$ is taken into account for some of these solutions.

\subsection{\label{sec:M5!=0}$M_{(5)}\neq0$ case}
In the case of $M_{(5)}\neq0$, there are many solutions of the critical point and they are written in terms of the parameters $\lambda$ and $\mu$. In the analysis, the Eqs.(\ref{eq:cont},\ref{eq:2},\ref{eq:5}) prevent $X_2=X_5=0$ which mean that $\phi$ cannot be decoupled from the system. The critical points of the relevant non-vanishing fields in the autonomous system are classified below.

\subsubsection{$\mathcal{F}_i^{(A)}$ : $(\phi,\eta_1,\rho_m)$-system (A) }
First, we consider $(\phi, \eta_1, \rho_m)$-system which leads to $X_1=X_3=X_6=X_8=0$, with the following constraint equation,
\begin{equation}
    1=X_2^2 - X_4^2 +X_5^2X_7^2 + X_5X_9^2-\lambda^2X_7^2. \label{eq:contA}
\end{equation}
In order to find the critical points, we use the constraint equation, Eq.(\ref{eq:contA}), to reduce the autonomous system by replacing a parameter for the remaining differential equation.\par
$\bullet$ \hypertarget{FA1}{$\mathcal{F}_1^{(A)}$}, we eliminate $X_2$ by using the constraint equation in Eq.(\ref{eq:contA}). The critical point is given by
\allowdisplaybreaks
\begin{align}
X_4=& -\frac{3}{\sqrt{\frac{8 \lambda ^2 \mu ^2}{3}-2}},
\quad X_5= \frac{\sqrt{4 \lambda ^2 \mu ^2-9}}{2 \mu },
\nonumber\\
X_7=& \frac{\sqrt{6} \mu }{\sqrt{4 \lambda ^2 \mu ^2-3}},
\quad X_9= -2\sqrt{2} \sqrt{\frac{\mu  \left(4 \lambda ^4 \mu ^4+15 \lambda ^2 \mu ^2-54\right)}{\left(4 \lambda ^2 \mu ^2-9\right)^{3/2} \left(4 \lambda ^2 \mu ^2-3\right)}}.
\end{align}
The effective equation of state of this system then simply reduces to
\begin{equation}
    w_\text{eff}=0.
\end{equation}
The critical point exists for all positive $\lambda$ and $\mu$ with $\lambda\mu > \frac{3}{2}$. After using the definition of stability matrix in Eq.(\ref{eq:eigenv}), the real part of eigenvalues is read
\begin{equation}
\left( 0, 0, 0, -\frac{9}{5} \right).
\end{equation}
This means that this critical point is a stable point. Then, the $\mathcal{F}_1^{(A)}$ can be represented as the late time DM dominated phase. The density parameters are given by
\begin{align}
\Omega_{\phi }=\frac{3 \left(4 \lambda ^2 \mu ^2-9\right)}{8 \lambda^2 \mu ^2-6},\quad
\Omega_{\eta _1}=\frac{3 \left(4 \lambda ^2 \mu ^2+9\right)}{6-8 \lambda ^2 \mu ^2},\quad
\Omega_{\rho _m}=\frac{4 \left(\lambda ^2 \mu ^2+6\right)}{4 \lambda ^2 \mu ^2-3}.
\end{align}
\begin{figure}[H]%
    \centering
    \subfloat[\centering Stability]{{\includegraphics[height=5.5cm]{section/Python_pic/FA1_Stab.png}}}%
    \subfloat[\centering Density parameters]{{\includegraphics[height=5.5cm]{section/Python_pic/omegaFA1.png}}}%
    \caption{On the left panel, the graph (a) shows the stability of the $\mathcal{F}_1^{(A)}$ critical point. The stable critical points lie in the blue shaded region while non-shaded region is forbidden for critical points. On the right panel, the graph (b) illustrates the relation between density parameter of non-vanishing fields and free parameters, i.e., $\lambda$ and $\mu$. The black, blue and red line in graph (b) indicates $\Omega_\phi$, $\Omega_{\eta_1}$ and $\Omega_{\rho_m}$, respectively.}
    \label{fig:plotFA1}%
\end{figure}

$\bullet$ \hypertarget{FA2}{$\mathcal{F}_2^{(A)}$}, the $X_4$ is eliminated by using Eq.(\ref{eq:contA}). The critical point is determined as
\begin{align}
    X_2= 0,\quad X_5= \lambda,\quad X_7= \frac{2 \mu }{\sqrt{\frac{8 \lambda ^2 \mu ^2}{3}-3}},\quad
    X_9= \frac{\sqrt{8 \lambda ^2 \mu ^2+18}}{\sqrt{\lambda  \left(8 \lambda ^2 \mu ^2-9\right)}}.
\end{align}
Having use Eq.(\ref{eq:eigenv}), the eigenvalues of the critical point are obtained. As a result, the eigenvalues are very complicated form and cannot be written in the analytical form. We therefore study them numerically. The stability of this critical point behaves as a stable node for all positive $\lambda$ and $\lambda\mu > \frac{3}{2\sqrt{2}}$. The plot of the stable region in $\lambda-\mu$ plane is depicted in Fig.(\ref{fig:plotFA2}). While the effective equation of state of this system then simply reduces to
\begin{equation}
    w_\text{eff}=0.
\end{equation}
Therefore, the {$\mathcal{F}_2^{(A)}$} system behaves as late time dominated DM in the universe. The results of the density parameters are shown below,
\begin{align}
\Omega _{\phi }=\frac{12 \lambda ^2 \mu ^2}{8 \lambda ^2 \mu ^2-9},\quad
\Omega _{\eta _1}=\frac{12 \lambda ^2 \mu ^2+27}{9-8 \lambda ^2 \mu ^2},\quad
\Omega _{\rho _m}=\frac{2 \left(4 \lambda ^2 \mu ^2+9\right)}{8 \lambda ^2 \mu ^2-9}.
\end{align}
\begin{figure}[H]%
    \centering
    \subfloat[\centering Stability]{{\includegraphics[height=5.5cm]{section/Python_pic/FA2_Stab.png}}}%
    \subfloat[\centering Density parameters]{{\includegraphics[height=5.5cm]{section/Python_pic/omegaFA2.png}}}%
    \caption{
    On the left panel, the graph (a) shows the stability of the $\mathcal{F}_2^{(A)}$ critical point. The stable critical points lie in the blue shaded region while non-shaded region is forbidden for critical points. On the right panel, the graph (b) illustrates the relation between density parameter of non-vanishing fields and free parameters, i.e., $\lambda$ and $\mu$. The black, blue and red line in graph (b) indicates $\Omega_\phi$, $\Omega_{\eta_1}$ and $\Omega_{\rho_m}$, respectively.}

    \label{fig:plotFA2}%
\end{figure}
\newpage
\subsubsection{$\mathcal{F}_i^{(B)}$ : $(\phi,\eta_1,A,\rho_m)$-system (B)}
In this sector, we consider $(\phi,\eta_1,A,\rho_m)$-system case i.e., $X_3=X_8=0$. The constraint equation becomes
\begin{equation}
1= X_2^2 - X_4^2 - \frac{X_1^2}{2 X_5^2} + X_5^2 X_7^2 - X_6^2 X_7^2 + X_5 X_9^2 - \lambda^2X_7^2. \label{eq:contB}
\end{equation}
The critical points of this system are listed below.\\

$\bullet$ \hypertarget{FB1}{$\mathcal{F}_1^{(B)}$}, we get rid of the $X_6$ variable by using constraint equation Eq.(\ref{eq:contB}). By doing this, the system provides 2 solutions of critical point. The first solution is given by
\begin{align}
    X_1=&0,\quad X_2=0,\quad X_4=-\frac{3 \sqrt{3}}{\sqrt{8 \lambda ^2 \mu ^2+9}},\quad X_5=\frac{\sqrt{8 \lambda ^2 \mu ^2+9}}{2 \sqrt{2} \mu }\nonumber \\
     X_7=&\frac{2 \sqrt{3} \mu }{\sqrt{8 \lambda ^2 \mu ^2+9}},\quad X_9=2^{3/4} \sqrt{\frac{\mu  \left(8 \lambda ^2 \mu ^2+63\right)}{\left(8 \lambda ^2 \mu ^2+9\right)^{3/2}}}.
\end{align}
We note that the existence of the critical point exists for all positive $\lambda$ and $\mu$. The eigenvalues of this critical points cannot write in the close form and numerical study is required. The stability analysis shows that this critical point is always saddle point. The effective equation of state of this system is given by
\begin{equation}
    w_\text{eff}=0.
\end{equation}
According to the dynamical system analysis, $\mathcal{F}_1^{(B)}$ can behave like DM phase or matter dominated era of the universe. Finally, the density parameters of the system are read,
\begin{align}
\Omega _A=&-\frac{81}{16 \lambda ^2 \mu ^2+18},\quad \Omega _{\phi }=\frac{3}{2},\nonumber\\
\Omega _{\eta _1}=&-\frac{3 \left(4 \lambda ^2 \mu ^2+9\right)}{8 \lambda ^2 \mu ^2+9}, \quad \Omega _{\rho _m}=\frac{8 \lambda ^2 \mu ^2+63}{8 \lambda ^2 \mu ^2+9}.
\end{align}
\newpage
\begin{figure}[H]%
    \centering
    \includegraphics[height=5.5cm]{section/Python_pic/omegaFB1.png}%
    \caption{
    The figure illustrates the relation between density parameter of non-vanishing fields and free parameters, i.e., $\lambda$ and $\mu$, for $\mathcal{F}_1^{(B)}$ critical point. Noted that the green, black, blue and red line in the graph indicates $\Omega_A$, $\Omega_\phi$, $\Omega_{\eta_1}$ and $\Omega_{\rho_m}$, respectively.
    }

    \label{fig:plotFB1}%
\end{figure}
$\bullet$ \hypertarget{FB2}{$\mathcal{F}_2^{(B)}$}, this is the second solution of the $(\phi,\eta_1,A,\rho_m)$ system where we eliminate $X_6$. It is written by
\begin{align}
    X_1=&0,\quad X_2=0,\quad X_9=0,\nonumber \\
    X_4=&\frac{1}{12 \lambda  \mu }\left(5 \sqrt{3}-\sqrt{8 \lambda ^2 \mu ^2+75}\right) \sqrt{\sqrt{3} \sqrt{8 \lambda ^2 \mu ^2+75}+15},\nonumber \\
    X_5=&-\frac{1}{2 \sqrt{3}}\sqrt{\frac{80 \lambda ^2 \mu ^2-61 \sqrt{3} \sqrt{8 \lambda ^2 \mu ^2+75}+915}{\mu ^2}},\nonumber \\
    X_7=&-\frac{1}{2 \lambda }\sqrt{\sqrt{\frac{8 \lambda ^2 \mu ^2}{3}+25}+5}.
\end{align}
The eigenvalues of the $\mathcal{F}_2^{(B)}$ system are very complicated. By numerical study, the stability of this critical point can be both stable and saddle points for all positive $\lambda$ and $\mu$. In addition, their existences of the stable node in the $\lambda-\mu$ plane are shown in Fig.(\ref{fig:plotFB2}). The effective equation of state of this system can be written as
\begin{equation}
    w_\text{eff}=\frac{1}{9} \left(6-\sqrt{3} \sqrt{8 \lambda ^2 \mu ^2+75}\right).
\end{equation}
The density parameters of the system are shown below
\begin{align}
\Omega_A=&\frac{1}{12} \left(39-5 \sqrt{3} \sqrt{8 \lambda ^2 \mu ^2+75}\right),\nonumber\\
\Omega_{\phi }=&\frac{5}{3} \sqrt{\frac{8 \lambda ^2 \mu ^2}{3}+25}-\frac{11}{6},\nonumber\\
\Omega_{\eta _1}=&-\frac{5}{36} \left(\sqrt{3} \sqrt{8 \lambda ^2 \mu ^2+75}+3\right),\nonumber\\
\Omega_{\rho _m}=&0.
\end{align}
\begin{figure}[H]%
    \centering
    \subfloat[\centering Stability]{{\includegraphics[height=5.5cm]{section/Python_pic/FB2_Stab.png}}}%
    \subfloat[\centering Effective equation of state]{{\includegraphics[height=5.5cm]{section/Python_pic/FB2_Weff.png}}}%
    \caption{
    On the left panel, the graph (a) shows the stability of the $\mathcal{F}_2^{(B)}$ critical point. The saddle critical points lie in the red shaded region while blue shaded region, $\lambda > 3/\mu$, is the region where critical points are stable. On the right panel, the graph (b) illustrates the changes of effective equation of state over $\lambda$ and $\mu$. Noted that the black line in graph (b) indicates the border between saddle and stable critical points.}

    \label{fig:plotFB2}%
\end{figure}

\begin{figure}[H]%
    \centering
    \includegraphics[height=5.5cm]{section/Python_pic/omegaFB2.png}%
    \caption{
    The figure illustrates the relation between density parameter of non-vanishing fields and free parameters, i.e., $\lambda$ and $\mu$, for $\mathcal{F}_2^{(B)}$ critical point. Noted that the green, black, blue and red line in the graph indicates $\Omega_A$, $\Omega_\phi$, $\Omega_{\eta_1}$ and $\Omega_{\rho_m}$, respectively.
    }

    \label{fig:plotFB2_2}%
\end{figure}

$\bullet$ \hypertarget{FB3}{$\mathcal{F}_3^{(B)}$}, the $X_9$ is excluded from the system by using Eq.(\ref{eq:contB}). We obtain the critical point as
    \begin{align}
    X_2=&0,\quad X_1=\frac{\sqrt{46 \lambda ^2 \mu ^2-\Delta-429}}{2 \sqrt{10} \mu }, \quad X_4= \frac{\sqrt{\mu ^2 \left(234-56 \lambda ^2 \mu ^2 - \Delta\right)}}{\mu  \sqrt{12 \lambda ^4 \mu ^4-68 \lambda ^2 \mu ^2+63}},\nonumber \\
    X_5=& \frac{1}{2 \sqrt{10}}\sqrt{\frac{86 \lambda ^2 \mu ^2-\Delta-269}{\mu ^2}},\quad X_6=\frac{\sqrt{46 \lambda ^2 \mu ^2-\Delta-429}}{2 \sqrt{10} \mu },\nonumber \\
    X_7=&\frac{\sqrt{\mu ^2 \left(234-56 \lambda ^2 \mu ^2-\Delta\right)}}{\sqrt{12 \lambda ^4 \mu ^4-68 \lambda ^2 \mu ^2+63}},
\end{align}
where
\begin{align}
    \Delta=\sqrt{6916 \lambda ^4 \mu ^4-47628 \lambda ^2 \mu ^2+74601}. \label{eq:Delta}
\end{align}
The eigenvalues of this system are also very complicated. With help from the numerical analysis, this critical point is always saddle point and the existence of the saddle point exists in the range of positive $\lambda$ and $\lambda\mu>\sqrt{\frac{57}{10}}$. The plot is depicted in Fig.(\ref{fig:plotFB3}). The effective equation of state of this system is given by
\begin{equation}
    w_\text{eff}=-\frac{5}{3}.
\end{equation}
The non-vanishing energy density parameters are read
\begin{align}
\Omega _A=&\frac{\left(429-46 \lambda ^2 \mu ^2-\Delta\right)\left(78 \lambda ^2 \mu ^2+\Delta-369\right)}{8 \left(2 \lambda ^2 \mu ^2-9\right) \left(86 \lambda ^2 \mu ^2+\Delta-269\right)},\nonumber\\
\Omega _{\phi }=&\frac{\left(-56 \lambda ^2 \mu ^2+\Delta+234\right)\left(86 \lambda ^2 \mu ^2+\Delta-269\right)}{40 \left(12 \lambda ^4 \mu ^4-68 \lambda ^2 \mu ^2+63\right)},\nonumber\\
\Omega _{\eta _1}=&\frac{\left(\lambda ^2 \mu ^2+1\right) \left(56 \lambda ^2 \mu ^2-\Delta-234\right)}{12 \lambda ^4 \mu ^4-68 \lambda ^2 \mu ^2+63},\nonumber\\
\Omega _{\rho _m}=&\frac{436 \lambda ^4 \mu ^4+2 \lambda ^2 \mu ^2 \left(3 \Delta-2998\right)-57 \Delta+18693}{2 \left(2 \lambda ^2 \mu ^2-9\right) \left(86 \lambda ^2 \mu ^2+\Delta-269\right)}.
\end{align}
\newpage
\begin{figure}[H]%
    \centering
    \subfloat[\centering Stability]{{\includegraphics[height=5.5cm]{section/Python_pic/FB3_Stab.png}}}%
    \subfloat[\centering Density parameters]{{\includegraphics[height=5.5cm]{section/Python_pic/omegaFB3.png}}}%
    \caption{
    On the left panel, the graph (a) shows the stability of the $\mathcal{F}_3^{(B)}$ critical point. The saddle critical points lie in the red shaded region while non-shaded region is forbidden for critical points. On the right panel, the graph (b) illustrates the relation between density parameter of non-vanishing fields and free parameters, i.e., $\lambda$ and $\mu$. The green, black, blue and red line in graph (b) indicates $\Omega_A$, $\Omega_\phi$, $\Omega_{\eta_1}$ and $\Omega_{\rho_m}$, respectively.}

    \label{fig:plotFB3}%
\end{figure}

\subsubsection{$\mathcal{F}_i^{(C)}$ : $(\phi,\eta_0,\eta_1,\rho_m)$-system (C)}
In this system, we have assumed that $X_1=0$ and $X_6=0$. The constraint equation is written by
\begin{equation}
    1 = X_2^2 - X_3^2 - X_4^2 + X_5^2 X_7^2 + X_5 X_9^2 - \lambda^2X_7^2  - \lambda^2X_8^2.  \label{eq:contC}
\end{equation}

$\bullet$ \hypertarget{FC1}{$\mathcal{F}_1^{(C)}$}, the critical point of the system is given by
\begin{align}
    X_3&= 0,\quad X_8=0, \quad X_4=\frac{3}{\sqrt{\frac{1}{3} \left(8 \lambda ^2 \mu ^2-6\right)}}, \quad X_5 =\frac{\sqrt{4 \lambda ^2 \mu ^2-9}}{2 \mu }, \nonumber\\
    X_7 &=-\frac{2 \mu }{\sqrt{\frac{1}{3} \left(8 \lambda ^2 \mu ^2-6\right)}}, \quad X_9 =\frac{2 \sqrt{\frac{27}{8 \lambda ^2 \mu ^2-6}+\frac{1}{2}}}{\sqrt{\frac{\sqrt{4 \lambda ^2 \mu ^2-9}}{\mu }}}.
\end{align}
Using definition of stability matrix in Eq.(\ref{eq:eigenv}), we obtain the eigenvalues and they are simply reduced to
\begin{equation}
    \left(0,0,0,-\frac{i}{2}  \sqrt{4 \lambda ^2 \mu ^2-9},\frac{i}{2} \sqrt{4 \lambda ^2 \mu ^2-9},-\frac{9}{5}\right)\quad \stackrel{{\rm real~ part}}{\longrightarrow} \quad\left(0,0,0,0,0,-\frac{9}{5}\right),
\end{equation}
for positive $\lambda$ and $\lambda\mu >\frac{3}{2}$. The effective equation of state reduces to zero,
\begin{equation}
    w_\text{eff}=0.
\end{equation}
This $(\phi,\eta_0,\eta_1,\rho_m)$ system behaves as DM. The density parameters are
\begin{align}
    \Omega_\phi &=\frac{3 \left(4 \lambda ^2 \mu ^2-9\right)}{8 \lambda ^2 \mu ^2-6},\quad \Omega_{\eta_0} =0,\nonumber \\
    \Omega_{\eta_1} &=\frac{3 \left(4 \lambda ^2 \mu ^2+9\right)}{6-8 \lambda ^2 \mu ^2},\quad \Omega_{\rho_m} =\frac{27}{4 \lambda ^2 \mu ^2-3}+1.
\end{align}
\begin{figure}[H]%
    \centering
    \subfloat[\centering Stability]{{\includegraphics[height=5.5cm]{section/Python_pic/FC1_Stab.png}}
    \subfloat[\centering Density parameters]{{\includegraphics[height=5.5cm]{section/Python_pic/omegaFC1.png}}}%
    \caption{
    On the left panel, the graph (a) shows the stability of the $\mathcal{F}_1^{(C)}$ critical point. The stable critical points lie in the blue shaded region while non-shaded region is forbidden for critical points. On the right panel, the graph (b) illustrates the relation between density parameter of non-vanishing fields and free parameters, i.e., $\lambda$ and $\mu$. The black, purple, blue and red line in graph (b) indicates $\Omega_\phi$, $\Omega_{\eta_0}$, $\Omega_{\eta_1}$ and $\Omega_{\rho_m}$, respectively.}

    \label{fig:plotFC1}%
\end{figure}

\subsubsection{$\mathcal{F}_i^{(D)}$ : $(\phi,\eta_0,\eta_1,A,\rho_m)$-system (D)}
In this case, we consider that all fields in the model are not vanished. Therefore, the constraint equation reads
\begin{equation}
1 = X_2^2 - X_3^2 - X_4^2 - \frac{X_1^2}{2 X_5^2} + X_5^2 X_7^2 - X_6^2 X_7^2 + X_5 X_9^2 - \lambda^2 X7^2 - \lambda^2X_8^2. \label{eq:contD}
\end{equation}

$\bullet$ \hypertarget{FD1}{$\mathcal{F}_1^{(D)}$}, we replace the $X_6$ variable in the system by using constraint equation, Eq.(\ref{eq:contD}). The critical point is given by
\begin{align}
    X_1=&0,\quad X_2=0,\quad X_3=0,\quad X_8=0,\nonumber \\
    X_4=&\frac{3 \sqrt{3}}{\sqrt{8 \lambda ^2 \mu ^2+9}},\quad
    X_5=\frac{\sqrt{8 \lambda ^2 \mu ^2+9}}{2 \sqrt{2} \mu },\nonumber \\
    X_7=&-\frac{2 \sqrt{3} \mu }{\sqrt{8 \lambda ^2 \mu ^2+9}},\quad
    X_9=-\frac{2^{3/4} \sqrt{\mu  \left(8 \lambda ^2 \mu ^2+63\right)}}{\left(8 \lambda ^2 \mu ^2+9\right)^{3/4}}.
\end{align}
The eigenvalues of this system are lengthy and complicated. According to the numerical study, the eigenvalues can be zero, positive and negative real numbers for all positive $\lambda$ and $\mu$. This means that the critical point is always the saddle point. While the effective equation of state of this system reads
\begin{align}
    w_\text{eff}=0.
\end{align}
As a result, the critical point of this system behaves like DM phase of matter dominated epoch of the observed universe. The density parameters are given by
\begin{align}
\Omega _A=&-\frac{81}{16 \lambda ^2 \mu ^2+18},\quad \Omega _{\phi }=\frac{3}{2},\quad \Omega _{\eta _0}=0,\nonumber\\
\Omega _{\eta _1}=&-\frac{3 \left(4 \lambda ^2 \mu ^2+9\right)}{8 \lambda ^2 \mu ^2+9},\quad \Omega _{\rho _m}=\frac{8 \lambda ^2 \mu ^2+63}{8 \lambda ^2 \mu ^2+9}.
\end{align}
\begin{figure}[H]%
    \centering
    \includegraphics[height=5.5cm]{section/Python_pic/omegaFD1.png}
    \caption{The figure illustrates the relation between density parameter of non-vanishing fields and free parameters, i.e., $\lambda$ and $\mu$, for $\mathcal{F}_1^{(D)}$ critical point. Noted that the green, black, purple, blue and red line in the graph indicates $\Omega_A$, $\Omega_\phi$, $\Omega_{\eta_0}$, $\Omega_{\eta_1}$ and $\Omega_{\rho_m}$, respectively.
    }

    \label{fig:plotFD1}%
\end{figure}
\newpage
$\bullet$ \hypertarget{FD2}{$\mathcal{F}_2^{(D)}$}, we eliminate the $X_9$ from the autonomous system by using Eq.(\ref{eq:contD}). The critical point is determined and it reads
\begin{align}
    X_2=&0,\quad X_3=0,\quad X_8=0,\quad X_1=\frac{\sqrt{\Delta+46 \lambda ^2 \mu ^2-429}}{2 \sqrt{10} \mu }, \nonumber \\
    X_4=&\frac{\sqrt{\mu ^2 \left(\Delta-56 \lambda ^2 \mu ^2+234\right)}}{\mu  \sqrt{12 \lambda ^4 \mu ^4-68 \lambda ^2 \mu ^2+63}},\quad X_5=\frac{1}{2 \sqrt{10}}\sqrt{\frac{\Delta+86 \lambda ^2 \mu ^2-269}{\mu ^2}},\nonumber \\ X_6=&\frac{\sqrt{\Delta+46 \lambda ^2 \mu ^2-429}}{2 \sqrt{10} \mu },\quad X_7=\frac{\sqrt{\mu ^2 \left(\Delta-56 \lambda ^2 \mu ^2+234\right)}}{\sqrt{12 \lambda ^4 \mu ^4-68 \lambda ^2 \mu ^2+63}},
\end{align}
where $\Delta$ is defined by Eq.(\ref{eq:Delta}). Due to the complexity, the eigenvalues are studied by numerical calculation. The critical point turns out to be saddle point for the range of $\lambda > 0$ and $\lambda\mu >\sqrt{\frac{57}{10}}$. The effective equation of state of this system is given by
\begin{equation}
    w_\text{eff}=-\frac{5}{3}.
\end{equation}
In addition, the density parameters are written in the following forms,
\begin{align}
\Omega _A=&\frac{\left(\Delta+46 \lambda ^2 \mu ^2-429\right) \left(-5 \Delta \left(6 \lambda ^2 \mu ^2-7\right)-\Delta^2+4576 \lambda ^4 \mu ^4-33828 \lambda ^2 \mu ^2+61686\right)}{40 \left(12 \lambda ^4 \mu ^4-68 \lambda ^2 \mu ^2+63\right) \left(\Delta+86 \lambda ^2 \mu ^2-269\right)},\nonumber\\
\Omega _{\phi}=&\frac{\left(\Delta-56 \lambda ^2 \mu ^2+234\right) \left(\Delta+86 \lambda ^2 \mu ^2-269\right)}{40 \left(12 \lambda ^4 \mu ^4-68 \lambda ^2 \mu ^2+63\right)},\nonumber\\
\Omega _{\eta _0}=&0,\nonumber \\
\Omega _{\eta_1}=&\frac{\left(\lambda ^2 \mu ^2+1\right) \left(-\Delta+56 \lambda ^2 \mu ^2-234\right)}{12 \lambda ^4 \mu ^4-68 \lambda ^2 \mu ^2+63},\nonumber\\
\Omega _{\rho_m}=&\frac{\Delta \left(36 \lambda ^4 \mu ^4-384 \lambda ^2 \mu ^2+399\right)-6 \Delta^2+2616 \lambda ^6 \mu ^6+2468 \lambda ^4 \mu ^4-131638 \lambda ^2 \mu ^2+316755}{2 \left(12 \lambda ^4 \mu ^4-68 \lambda ^2 \mu ^2+63\right) \left(\Delta+86 \lambda ^2 \mu ^2-269\right)}.
\end{align}
\newpage
\begin{figure}[H]%
    \centering
    \subfloat[\centering Stability]{{\includegraphics[height=5.5cm]{section/Python_pic/FD2_Stab.png}}}%
    \subfloat[\centering Density parameters]{{\includegraphics[height=5.5cm]{section/Python_pic/omegaFD2.png}}}%
    \caption{On the left panel, the graph (a) shows the stability of the $\mathcal{F}_2^{(D)}$ critical point. The saddle critical points lie in the red shaded region while non-shaded region is forbidden for critical points. On the right panel, the graph (b) illustrates the relation between density parameter of non-vanishing fields and free parameters, i.e., $\lambda$ and $\mu$. The green, black, purple, blue and red line in graph (b) indicates $\Omega_A$, $\Omega_\phi$, $\Omega_{\eta_0}$, $\Omega_{\eta_1}$ and $\Omega_{\rho_m}$, respectively.}

    \label{fig:plotFD2}%
\end{figure}

\subsection{\label{sec:MB5!=0}$M_{(5)}\neq0$ case with slow-roll scalar fields}
From now on, we will assume that the slow-roll approximation is used in order to evaluate the critical point in this case. In the other word, the kinetic term of $\eta_0$ and $\eta_1$ are vanished ($X_4=0$ and $X_3=0$), if any of these fields exists in the considering system. Contrary to the previous subsection, the dilaton field can be decoupled from the autonomous.
\subsubsection{$\bar{\mathcal{F}}_i^{(B)}$ : $(\phi,\eta_1,A,\rho_m)$-system (B)}
In this case, we consider $(\phi,\eta_1,A,\rho_m)$-system and this leads to $X_3=0$, $X_4=0$ and $X_8=0$. The constraint equation is given by
\begin{equation}
    1 = X_2^2 - \frac{X_1^2}{2 X_5^2} + X_5^2 X_7^2 - X_6^2 X_7^2 + X_5 X_9^2 - \lambda^2X_7^2 \label{eq:contBs}
\end{equation}

$\bullet$ \hypertarget{FBbar1}{$\bar{\mathcal{F}}_1^{(B)}$}, the $X_1$ variable is eliminated by using Eq.(\ref{eq:contBs}). The critical point is read
\begin{align}
X_2=&0,\quad X_9=0,\nonumber \\
    X_5=& \sqrt{\frac{1}{458} \left(29 \sqrt{273}+547\right)} \lambda,\nonumber \\
    X_6 =& \sqrt{\frac{2}{229} \left(7 \sqrt{273}-18\right)} \lambda ,\nonumber \\
    X_7 =& \sqrt{\frac{3 \sqrt{273}}{8 \lambda ^2}-\frac{25}{8 \lambda ^2}}.
\end{align}
The eigenvalues have been studied with numerical calculation due to their complication. The result shows that the critical point is always saddle point for all positive $\lambda$ and $\mu$. Effective equation of state is written by
\begin{equation}
    w_\text{eff} =-1.
\end{equation}
As a result, we might interpret that the critical point can be identified as DE (de-Sitter) phase at the late time. Non-vanishing density parameters are given by
\begin{align}
    \Omega_A &=\frac{1}{2} \left(\sqrt{273}-21\right),\quad \Omega_\phi =\frac{1}{4} \left(\sqrt{273}+11\right),\nonumber \\
    \Omega_{\eta_1} &=\frac{1}{8} \left(25-3 \sqrt{273}\right), \quad \Omega_{\rho_m} =0.
\end{align}

$\bullet$ \hypertarget{FBbar2}{$\bar{\mathcal{F}}_2^{(B)}$}, we use constraint equation in Eq.(\ref{eq:contBs}) to get rid of $X_2$ from the system. Then, the critical point of this system is
\begin{align}
    X_1 =0,\quad X_6=0,\quad X_9=0,\quad
    X_5 =\frac{\sqrt{31}\lambda}{5},\quad X_7=\frac{5}{\sqrt{6}\lambda}.
\end{align}
The eigenvalues are given by
\begin{equation}
    \left(0,-\frac{24}{5},-\frac{3}{2},-1,1\right),
\end{equation}
According to the eigenvalues given above, it clearly shows that the critical point is the saddle point. The effective equation of state of the system takes the following form,
\begin{equation}
w_\text{eff}=-1.
\end{equation}
The critical point can behave as DE but the universe is not dominated by DE at the late time. The density parameters are written by
\begin{align}
    \Omega_A =0,\quad \Omega_\phi =\frac{31}{6},\quad
    \Omega_{\eta_1} =-\frac{25}{6}, \quad \Omega_{\rho_m} =0.
\end{align}

$\bullet$ \hypertarget{FBbar3}{$\bar{\mathcal{F}}_3^{(B)}$}, we exclude the $X_6$ from the system via the constraint equation in (\ref{eq:contBs}). The critical point is determined by the following form,
\begin{align}
    X_1=0,\quad X_2=0,\quad X_9=0,\quad
    X_5 =\sqrt{\frac{13}{5}} \lambda,\quad X_7 =\frac{1}{\lambda }\sqrt\frac{5}{2}.
\end{align}
The real part of eigenvalues is read
\begin{eqnarray}
    &&\left(-\frac{3}{2},1,-4.59,-1.60-3.56 i,-1.60+3.56 i\right)\nonumber\\\quad \stackrel{{\rm real~part}}{\longrightarrow}\quad &&\left(-\frac{3}{2},1,-4.59,-1.60,-1.60\right).
\end{eqnarray}
The results show that this critical point always represents saddle point for all positive $\lambda$ and $\mu$. The effective equation of state is reduced to
\begin{equation}
w_\text{eff}=-1.
\end{equation}
The density parameters can be written as
\begin{align}
    \Omega_A =-3,\quad \Omega_\phi =\frac{13}{2},\quad
    \Omega_{\eta_1} =-\frac{5}{2},\quad \Omega_{\rho_m} =0.
\end{align}

$\bullet$ \hypertarget{FBbar4}{$\bar{\mathcal{F}}_4^{(B)}$}, we find two solutions via eliminating the $X_7$ by Eq.(\ref{eq:contBs}). The first solution reads
\begin{align}
    X_1 =0,\quad X_2=0,\quad X_6=0,\quad X_9=0,\quad X_5 =\sqrt{\frac{7}{5}} \lambda.
\end{align}
The eigenvalues of stability matrix of this critical points are given by
\begin{equation}
    \left(0,-\frac{24}{5},-\frac{3}{2},-1,1\right).
\end{equation}
As a result, this critical point is always saddle point for all positive $\lambda$ and $\mu$. The effective equation of state is determined as
\begin{equation}
w_\text{eff}=-\frac{1}{3}.
\end{equation}
Non-vanishing density parameters are written by
\begin{align}
    \Omega_A =0,\quad \Omega_\phi =\frac{7}{2},\quad
    \Omega_{\eta_1} =-\frac{5}{2},\quad\Omega_{\rho_m} =0.
\end{align}

$\bullet$ \hypertarget{FBbar5}{$\bar{\mathcal{F}}_5^{(B)}$} the second solution of the system with eliminating $X_7$ by Eq.(\ref{eq:contBs}) is given by
\begin{align}
    X_2 &=0,\quad X_9=0,\nonumber\\
    X_1 &=-\sqrt{\frac{1}{373} \left(23 \sqrt{537}+81\right)} \lambda,\nonumber\\
    X_5 &=-\sqrt{\frac{1}{746} \left(47 \sqrt{537}+1171\right)} \lambda,\nonumber\\
    X_6 &=-\sqrt{\frac{1}{373} \left(23 \sqrt{537}+81\right)} \lambda.
\end{align}
The eigenvalues of the system are very complicated and we need to study the stability numerically. The result reflects that this critical point is always stable node. The effective equation of state is found to be
\begin{equation}
w_\text{eff}=-\frac{5}{3}.
\end{equation}
As a result, this critical point can be DE dominated universe at the late time
and the density parameters are found as
\begin{align}
    \Omega_A &=\frac{1}{2} \left(\sqrt{537}-33\right),\quad \Omega_\phi =\frac{1}{4} \left(\sqrt{537}+17\right),\nonumber \\
    \Omega_{\eta_1} &=\frac{1}{8} \left(43-3 \sqrt{537}\right),\quad \Omega_{\rho_m} =0.
\end{align}

\subsubsection{$\bar{\mathcal{F}}_i^{(C)}$ : $(\phi,\eta_0,\eta_1,\rho_m)$-system (C)}
We consider $(\phi,\eta_0,\eta_1,\rho_m)$-system, we assume that $X_1=0$ $,X_3=0$, $X_4=0$ and $X_6=0$. The constraint equation takes the following form,
\begin{equation}
    1 = X_2^2 + X_5^2 X_7^2 + X_5 X_9^2 - \lambda^2X_7^2- \lambda^2X_8^2.  \label{eq:contCs}
\end{equation}

$\bullet$ \hypertarget{FCbar1}{$\bar{\mathcal{F}}_1^{(C)}$}, we firstly get rid of $X_2$ from the system by the constraint equation (\ref{eq:contCs}), the critical point is given by
\begin{align}
X_5 = \frac{\sqrt{31} \lambda }{5}, \quad X_7 = \frac{5}{\sqrt{6} \lambda },\quad X_8 =0,\quad X_9=0.
\end{align}
The eigenvalues of this critical point are found as
\begin{equation}
\left(0,0,-\frac{24}{5},-\frac{3}{2}\right).
\end{equation}
The result of the eigenvalues clearly shows that the critical point is stable node. The effective equation of state is read
\begin{equation}
    w_\text{eff}=-1.
\end{equation}
This means that this critical point can be used to describe de-Sitter DE dominated at the late time. While the density parameters are given by
\begin{align}
    \Omega_\phi =\frac{31}{6},\quad
    \Omega_{\eta_0} =0,\quad
    \Omega_{\eta_1} =-\frac{25}{6},\quad
    \Omega_{\rho_m} =0.
\end{align}

$\bullet$ \hypertarget{FCbar2}{$\bar{\mathcal{F}}_2^{(C)}$}, we exclude the $X_7$ from the autonomous system in this case. The critical point is found to be
\begin{align}
X_2 =0,\quad X_8=0,\quad X_9=0,\quad X_5 = \sqrt{\frac{7}{5}}\lambda.
\end{align}
Noted that the critical point exists for all positive $\lambda$ parameter. The real part of the eigenvalues is given by
\begin{equation}
\left(-\frac{1}{2},1,-1-i \sqrt{\frac{37}{5}},-1+i \sqrt{\frac{37}{5}}\right)\quad \stackrel{{\rm real \ part}}{\longrightarrow}\quad \left(-\frac{1}{2},1,-1,-1\right).
\end{equation}
The result explicitly shows that this critical point is saddle point of the system. The effective equation of state is equal to
\begin{equation}
    w_\text{eff}=-\frac{1}{3}.
\end{equation}
The density parameters are found that
\begin{align}
    \Omega_\phi =\frac{7}{2},\quad
    \Omega_{\eta_0} =0,\quad
    \Omega_{\eta_1} =-\frac{5}{2},\quad
    \Omega_{\rho_m} =0.
\end{align}

$\bullet$ \hypertarget{FCbar3}{$\bar{\mathcal{F}}_3^{(C)}$}, we eliminate the $X_9$ from the system by Eq.(\ref{eq:contCs}). The critical point is determined as
\begin{align}
X_2 =0,\quad
X_5 = \sqrt{\frac{5}{3}}\lambda, \quad
X_7 = \frac{3}{\sqrt{2} \lambda }, \quad
X_8 = 0.
\end{align}
The eigenvalues of this critical point are very complicated. Then, the numerical calculation is required and the result shows that the critical point is stable for all positive $\lambda$ and $\mu$. The effective equation of state is
\begin{equation}
    w_\text{eff}=-1,
\end{equation}
which can be interpreted as late time dominated DE phase. The density parameters of this critical point are given by
\begin{align}
    \Omega_\phi =\frac{15}{2},\quad
    \Omega_{\eta_0} =0,\quad
    \Omega_{\eta_1} =-\frac{9}{2},\quad
    \Omega_{\rho_m} =-2.
\end{align}

\subsubsection{$\bar{\mathcal{F}}_i^{(D)}$ : $(\phi,\eta_0,\eta_1,A,\rho_m)$-system (D)}
In this case, all fields are considered except the kinetic terms of $\eta_0$ and $\eta_1$, i.e. $X_4=0$ and $X_3=0$. The constraint equation reads
\begin{equation}
    1 = X_2^2 - \frac{X_1^2}{2 X_5^2} + X_5^2 X_7^2 - X_6^2 X_7^2 + X_5 X_9^2 - \lambda^2X_7^2 - \lambda^2X_8^2. \label{eq:contDs}
\end{equation}

$\bullet$ \hypertarget{FDbar1}{$\bar{\mathcal{F}}_1^{(D)}$}, we firstly replace the $X_1$ by constraint equation (\ref{eq:contDs}). The critical point of this system is given by
\begin{align}
    X_2 =&0,\quad  X_8=0,\quad X_9=0,\nonumber\\
    X_5 =& \sqrt{\frac{1}{458} \left(29 \sqrt{273}+547\right)} \lambda ,\nonumber \\
    X_6 =& \sqrt{\frac{2}{229} \left(7 \sqrt{273}-18\right)} \lambda , \nonumber\\
    X_7 =& \sqrt{\frac{3 \sqrt{273}}{8 \lambda ^2}-\frac{25}{8 \lambda ^2}},\nonumber\\
\end{align}
In addition, the critical point exists for all positive $\lambda$ and $\mu$. According to the complicated form of the eigenvalues, the numerical analysis demonstrates that this critical point is stable. The effective equation of state is equal to
\begin{equation}
    w_\text{eff} = -1,
\end{equation}
which corresponds to the DE phase dominated at the late time. The density parameters of the system are given by
\begin{align}
    \Omega_A &=\frac{1}{2} \left(\sqrt{273}-21\right), \nonumber\\
    \Omega_\phi &=\frac{1}{2} \left(\sqrt{273}+11\right), \nonumber\\
    \Omega_{\eta_1} &=\frac{1}{8} \left(25-3 \sqrt{273}\right), \nonumber\\
        \Omega_{\eta_0} &=0,\quad\Omega_{\rho_m} =0 .
\end{align}

$\bullet$ \hypertarget{FDbar2}{$\bar{\mathcal{F}}_2^{(D)}$}, getting rid of the $X_2$ from the system via constraint equation (\ref{eq:contDs}), the critical point is found as follow
\begin{align}
    X_1 &=0, \quad X_6=0,\quad X_8=0,\quad X_9=0 ,\quad
    X_5 =\frac{\sqrt{31}\lambda}{5} ,\quad X_7=\frac{5}{\sqrt{6}\lambda}.
\end{align}
The eigenvalues of the critical point are equal to
\begin{align}
\left(0,0,-\frac{24}{5},-\frac{3}{2},-1,1\right),
\end{align}
which come with all positive $\lambda$ and $\mu$. The effective equation of state is given by
\begin{equation}
    w_\text{eff}=-1.
\end{equation}
Although this critical point behaves as DE due to its equation of state but it does not dominate the universe at the late time. While the density parameters are read
\begin{align}
    \Omega_A =0, \quad \Omega_\phi =\frac{31}{6}, \quad
    \Omega_{\eta_0} =0, \quad
    \Omega_{\eta_1} =-\frac{25}{6}, \quad
    \Omega_{\rho_m} =0 .
\end{align}

$\bullet$ \hypertarget{FDbar3}{$\bar{\mathcal{F}}_3^{(D)}$}, we eliminate the $X_7$ from the system by using Eq.(\ref{eq:contDs}). The critical point of this system equals to
\begin{align}
    X_1=0,\quad X_2=0,\quad X_6=0,\quad X_8=0,\quad X_9=0,\quad X_5 =\sqrt{\frac{7}{5}} \lambda.
\end{align}
The eigenvalues of the critical point are
\begin{align}
\left(0,-\frac{9}{5},-1,-\frac{1}{2},1,2\right).
\end{align}
This reflects that this critical point is always saddle point. The effective equation of state is given by
\begin{equation}
    w_\text{eff}=-\frac{1}{3}.
\end{equation}
The density parameters of this system are written by
\begin{align}
    \Omega_A =0, \quad \Omega_\phi=\frac{7}{2}, \quad \Omega_{\eta_0} =0,\quad \Omega_{\eta_1} =-\frac{5}{2}, \quad \Omega_{\rho_m} =0 .
\end{align}

$\bullet$ \hypertarget{FDbar4}{$\bar{\mathcal{F}}_4^{(D)}$}, as done previously for eliminating $X_7$, the critical point reads
\begin{align}
X_2 &=0, \quad X_8=0,\quad X_9=0,\nonumber \\
    X_1 &=\sqrt{\frac{1}{373} \left(23 \sqrt{537}+81\right)} ,\nonumber\\
    X_5 &=\sqrt{\frac{1}{746} \left(47 \sqrt{537}+1171\right)},\nonumber\\
    X_6 &=\sqrt{\frac{1}{373} \left(23 \sqrt{537}+81\right)} \lambda.
\end{align}
The stability analysis of this critical point is performed numerically and the result shows that this critical point is stable for all positive $\lambda$. While the effective equation of state is equal to
\begin{equation}
    w_\text{eff}=-\frac{5}{3}.
\end{equation}
The density parameters of this system are given by
\begin{align}
    \Omega_A &=\frac{1}{2} \left(\sqrt{537}-33\right), \nonumber\\
    \Omega_\phi &=\frac{1}{4} \left(\sqrt{537}+17\right), \nonumber\\
    \Omega_{\eta_1} &=\frac{1}{8} \left(43-3 \sqrt{537}\right), \nonumber\\
        \Omega_{\eta_0} &=0, \quad \Omega_{\rho_m} =0 .
\end{align}
\section{\label{sec:Dynamical System}Dynamical System}
Now, we arrive at the crucial part of this work that is to study the dynamics of the universe. The main purpose of this work is to demonstrate the existence of DM and DE in the KK inspired BD gravity with the inclusion of the barotropic fluid. Since the dynamical system is widely used and succeed to qualitatively analyse the cosmological models \cite{Bahamonde:2017ize}, we will apply this method to the model. For simplicity, we define the dimensionless parameters to construct the autonomous system of the equations,

\begin{eqnarray}
&& X_1 = \frac{\sqrt{3}a}{H}\mathcal{A}_c,\quad X_2 = \sqrt{\frac{2}{3}}\frac{\dot{\phi}}{H\phi},\quad X_3 = \frac{\dot{\eta}_0}{H},\quad X_4 = \frac{\dot{\eta}_1}{H},\quad
X_5 = \frac{1}{\phi},\nonumber \\
&& X_6 = \sqrt{3} a A,\quad X_7 = \frac{\eta_1}{H R_k},\quad X_8 = \frac{\eta_0}{H R_k}, \quad X_9=\sqrt{\frac{\rho_m}{H^2}}.
\end{eqnarray}
We can use the Friedmann equation as a constraint equation of the the dimensionless parameters, i.e.,
\begin{equation}
1=-\frac{X_1^2}{2 X_5^2}+X_2^2-X_3^2-X_4^2-X_6^2 X_7^2+X_5 \left(X_5 X_7^2+X_9^2\right)-M_{(5)}^2 R_k^2. \left(X_7^2+X_8^2\right)\label{eq:cont}
\end{equation}
Noted that the additional free parameters of the system are defined by
\begin{equation}
    \lambda \equiv M_{(5)}R_k \quad\text{and}\quad \mu \equiv \frac{1}{H R_k}.
\end{equation}
Then the conformal derivative of the defined parameters are derived as following
\begin{align}
\frac{1}{H}\frac{dX_1}{dt}=&\frac{31 X_1^3}{60 X_5^2}-\frac{3 X_1 X_2^2}{10}+\frac{3 X_1 X_3^2}{10}+\frac{3 X_1 X_4^2}{10}+\frac{1}{10} X_1 X_5^2 X_7^2-\frac{13}{30} X_1 X_6^2 X_7^2\,,\nonumber \\
&-\frac{7}{10} \lambda ^2 X_1 X_7^2-\frac{7}{10} \lambda ^2 X_1 X_8^2+\frac{17 X_1}{5}\label{eq:1}\\
\frac{1}{H}\frac{dX_2}{dt}=&\frac{31 X_1^3}{60 X_5^2}-\frac{3 X_1 X_2^2}{10}+\frac{3 X_1 X_3^2}{10}+\frac{3 X_1 X_4^2}{10}+\frac{1}{10} X_1 X_5^2 X_7^2-\frac{13}{30} X_1 X_6^2 X_7^2\nonumber \\
&-\frac{7}{10} \lambda ^2 X_1 X_7^2-\frac{7}{10} \lambda ^2 X_1 X_8^2+\frac{17 X_1}{5}\,,\label{eq:2}\\
\frac{1}{H}\frac{dX_3}{dt}=&\frac{31 X_1^2 X_3}{60 X_5^2}-\frac{3 X_2^2 X_3}{10}-\sqrt{\frac{3}{2}} X_2 X_3+\frac{3 X_3^3}{10}+\frac{3 X_3 X_4^2}{10}+\frac{1}{10} X_3 X_5^2 X_7^2\nonumber \\
&-\frac{13}{30} X_3 X_6^2 X_7^2-\frac{7}{10} \lambda ^2 X_3 X_7^2-\frac{7}{10} \lambda ^2 X_3 X_8^2-\frac{3 X_3}{5}-\lambda ^2 \mu  X_8\,,\label{eq:3}\\
\frac{1}{H}\frac{dX_4}{dt}=&\frac{31 X_1^2 X_4}{60 X_5^2}-\frac{3 X_2^2 X_4}{10}-\sqrt{\frac{3}{2}} X_2 X_4+\frac{3 X_3^2 X_4}{10}+\frac{3 X_4^3}{10}+\frac{1}{10} X_4 X_5^2 X_7^2-\lambda ^2 \mu  X_7\nonumber \\
&-\frac{13}{30} X_4 X_6^2 X_7^2-\frac{7}{10} \lambda ^2 X_4 X_7^2-\frac{7}{10} \lambda ^2 X_4 X_8^2-\frac{3 X_4}{5}+\mu  X_5^2 X_7-\mu  X_6^2 X_7\,,\label{eq:4}\\
\frac{1}{H}\frac{dX_5}{dt}=&-\sqrt{\frac{3}{2}} X_2 X_5\,,\label{eq:5}\\
\frac{1}{H}\frac{dX_6}{dt}=&X_1-X_6\,,\label{eq:6}\\
\frac{1}{H}\frac{dX_7}{dt}=&\frac{31 X_1^2 X_7}{60 X_5^2}-\frac{3 X_2^2 X_7}{10}+\frac{3 X_3^2 X_7}{10}+\frac{3 X_4^2 X_7}{10}+\mu  X_4+\frac{X_5^2 X_7^3}{10}-\frac{13 X_6^2 X_7^3}{30}\,,\nonumber \\
&+\frac{12 X_7}{5}-\frac{7 \lambda ^2 X_7^3}{10}-\frac{7}{10} \lambda ^2 X_7 X_8^2\,,\label{eq:7}\\
\frac{1}{H}\frac{dX_8}{dt}=&\frac{31 X_1^2 X_8}{60 X_5^2}-\frac{3 X_2^2 X_8}{10}+\frac{3 X_3^2 X_8}{10}+\mu  X_3+\frac{3 X_4^2 X_8}{10}+\frac{1}{10} X_5^2 X_7^2 X_8\nonumber\\
&-\frac{13}{30} X_6^2 X_7^2 X_8-\frac{7}{10} \lambda ^2 X_7^2 X_8-\frac{1}{10} 7 \lambda ^2 X_8^3+\frac{12 X_8}{5}\,,\label{eq:8}\\
\frac{1}{H}\frac{dX_9}{dt}=&\frac{31 X_1^2 X_9}{60 X_5^2}-\frac{3 X_2^2 X_9}{10}+\frac{3 X_3^2 X_9}{10}+\frac{3 X_4^2 X_9}{10}+\frac{1}{10} X_5^2 X_7^2 X_9-\frac{13}{30} X_6^2 X_7^2 X_9\nonumber\\
&-\frac{7}{10} \lambda^2 X_7^2 X_9-\frac{7}{10} \lambda^2 X_8^2 X_9+\frac{9 X_9}{10}\,,\label{eq:9} \\
\frac{1}{H}\frac{d\mu}{dt}=&\frac{12 \mu }{5}+\frac{31 \mu  X_1^2}{60 X_5^2}-\frac{3 \mu X_2^2}{10}+\frac{3 \mu X_3^2}{10}+\frac{3 \mu X_4^2}{10}+\frac{1}{10} \mu X_5^2 X_7^2-\frac{13}{30} \mu X_6^2 X_7^2\nonumber \\
&-\frac{7}{10} \lambda ^2 \mu X_7^2-\frac{7}{10} \lambda ^2 \mu X_8^2. \label{eq:mu}
\end{align}
We obtain a critical point of the dynamical system by setting all derivative equations from Eqs.(\ref{eq:1}-\ref{eq:9}) equal to zero with constraint Eq.(\ref{eq:cont}). The stability of a critical point can be analysed from its eigenvalues of the matrix
\begin{equation}
    M_{ij}=\frac{d}{dX_i}\left( \frac{1}{H} \frac{dX_j}{dt} \right) \label{eq:eigenv}.
\end{equation}
If all eigenvalues are negative, the critical point is stable and if any is positive, the critical point is unstable. The energy density of each specie can be written as

\begin{align}
    &\Omega_{A} = -\frac{X_1^2}{2X_5^2} -     X_6^2 X_7^2,\quad \Omega_{\phi} = X_2^2 + X_5^2 X_7^2,\quad
    \Omega_{\eta_0} = - \left(X_3^2 + \lambda^2 X_8^2\right),\nonumber\\
    &\Omega_{\eta_1} = - \left(X_4^2 + \lambda^2 X_7^2\right),\quad
    \Omega_{\rho_m} = X_5X_9^2.
    \label{density}
\end{align}
The effective equation of state is given by
\begin{equation}
    w_\text{eff}=-1-\frac{2\dot{H}}{3H^2}.
    \label{w-eff}
\end{equation}
We use the above equation to identify the phase of the universe, where $w_\text{eff}<-1/3$ and $w_\text{eff}>-1/3$ represent the acceleration and deceleration expansion of the universe, respectively. While $w_\text{eff}=0$ represents the matter dominated phase of the universe, $w_\text{eff}<-1$ corresponds to the phantom DE. In addition, one may find the solutions of the scale factor by integrating the Eq.(\ref{w-eff}).
\section{\label{sec:Lower Mode Cases}Lower Mode Cases}
It is well known that the KK scalar fields are composed of the infinite tower mode. In particular, a higher mode with a momentum larger than the reduced Plank scale will be physically irrelevant to the observable universe. Therefore, we can neglect such higher modes and the remaining modes of the KK excited field should satisfy the momentum condition:
\begin{equation}
n < R_k\,.
\end{equation}
Therefore, we focus to study the KK inspired BD model for, $1 < R_k < 2$. In this case, only zero mode and the first excited mode involve in the system. For simplicity and convenient to study the dynamics of the universe by using the dynamical system analysis with the dimensionless parameters, we can simply drop the radiation matter term in this work. Moreover, we found that the inclusion of radiation part would lead to an unnecessary introduction of free parameters. The Friedmann equation for the lower mode of the model reads,
\begin{eqnarray}
H^{2} = -\frac{3 \phi^2a^2}{2}\mathcal{A}_c^2 +\frac{2\dot{\phi}^2}{3\phi^2} - \dot{\eta}_0^2 - \dot{\eta}_1^2 - M_{(5)}^2 \eta_0^2 - M_{(5)}^2 \eta_1^2  +  \left(\frac{1}{\phi^2} - 3 a^2 A^2\right) \frac{\eta_1^2}{R_k^2} + \frac{\rho_m}{\phi}. \label{eq:Friedman2}
\end{eqnarray}
Noted that $\eta_0$ and $\eta_1$ scalar fields are the zero mode and the first excited mode, respectively. The equation of motion for the dilation field, $\phi$ becomes
\begin{align}
    \ddot{\phi} + 3 H \dot{\phi} =& -\frac{9}{2}\phi\left(\dot{H} + 2H^2\right) - \frac{27}{8}\phi^3a^2 \mathcal{A}_c^2 + \frac{\dot{\phi}^2}{2\phi} - \frac{3}{4}\phi\dot{\eta}_0^2 - \frac{3}{4}\phi\dot{\eta}_1^2 \nonumber \\
& + \frac{3}{4} \phi \left( \frac{1}{\phi^2} + 3 a^2 A^2 \right)\frac{\eta_1^2}{R_k^2} + \frac{3}{4}\phi M_{(5)}^2\eta_0^2 + \frac{3}{4}\phi M_{(5)}^2\eta_1^2 \label{eq:seteomphi}.
\end{align}
The pressure can be written as $p_i=w_i\rho_i$ where $w_m =0$. Solving Eq.(\ref{eq:EFEij}) and Eq.(\ref{eq:seteomphi}) gives
\begin{align}
    \dot{H} =& \frac{13  a^2 A^2 \eta_1^2}{10 R_k^2}-\frac{31}{20} a^2 \phi^2 \mathcal{A}_c^2-\frac{3}{10}  \dot{\eta}_1^2+\frac{7}{10} M_{(5)}^2  \eta_1^2-\frac{  \eta_1^2}{10 R_k^2 \phi^2}\nonumber \\
    & -\frac{12 H^2}{5}-\frac{3}{10}  \dot{\eta}_0^2+\frac{7}{10} M_{(5)}^2  \eta_0^2+\frac{\dot{\phi}^2}{5 \phi^2} ,
\end{align}
and
\begin{align}
    \ddot{\phi} =& -\frac{18  a^2 A^2 \eta_1^2 \phi}{5 R_k^2}+\frac{18}{5} a^2 \phi^3 \mathcal{A}_c^2+\frac{3}{5}   \phi \dot{\eta}_1^2-\frac{12}{5} M_{(5)}^2  \eta_1^2 \phi + \frac{6   \eta_1^2}{5 R_k^2 \phi}\nonumber \\
    & -3 H \dot{\phi} + \frac{9}{5} H^2 \phi + \frac{3}{5}   \phi \dot{\eta}_0^2-\frac{12}{5} M_{(5)}^2  \eta_0^2 \phi - \frac{2 \dot{\phi}^2}{5 \phi},
\end{align}
where we define $\mathcal{A}_c = \dot{A} + 2 H A$\,.

The specific form of equations and fields allow us to analyse the system both analytically and numerically. In the next section, we will use the results to study the dynamical system of the model.