# Learned Tree Search for Long-Horizon Social Robot Navigation in Shared
Airspace

Ingrid Navarro1∗, Jay Patrikar1∗, Joao P. A. Dantas2, Rohan Baijal3, Ian
Higgins1,
Sebastian Scherer1 and Jean Oh1 1Authors are with the Robotics Institute,
Carnegie Mellon University, Pittsburgh, PA, USA. {ingridn, jaypat, ihiggins,
basti<EMAIL_ADDRESS>is with the Institute for Advanced Studies,
Sao Jose dos Campos, SP, Brazil. Work done as a visiting researcher at
Carnegie Mellon University<EMAIL_ADDRESS>is with the Indian
Institute of Technology Kanpur, Kanpur, Uttar Pradesh, India. Work done as a
visiting researcher at Carnegie Mellon University<EMAIL_ADDRESS>contribution.

###### Abstract

The fast-growing demand for fully autonomous aerial operations in shared
spaces necessitates developing trustworthy agents that can safely and
seamlessly navigate in crowded, dynamic spaces. In this work, we propose
Social Robot Tree Search (SoRTS), an algorithm for the safe navigation of
mobile robots in social domains. SoRTS aims to augment existing socially-aware
trajectory prediction policies with a Monte Carlo Tree Search planner for
improved downstream navigation of mobile robots. To evaluate the performance
of our method, we choose the use case of social navigation for general
aviation. To aid this evaluation, within this work, we also introduce
X-PlaneROS, a high-fidelity aerial simulator, to enable more research in full-
scale aerial autonomy. By conducting a user study based on the assessments of
26 FAA-certified pilots, we show that SoRTS performs comparably to a competent
human pilot, significantly outperforming our baseline algorithm. We further
complement these results with self-play experiments, showcasing our
algorithm’s behavior in scenarios with increasing complexity. [Code]111
Codebase: https://github.com/cmubig/sorts$\mid$ [Video]222 Video:
https://youtu.be/PBE3O4cW2rI

##  I Introduction

A social robot strives to synthesize decision policies that enable it to
seamlessly interact with humans, ensuring social compliance while attaining
its desired goal. While marked progress has been made in social navigation and
motion prediction [1, 2, 3], achieving seamless navigation among humans while
balancing social and self-interested objectives remains challenging. Social
navigation can be formulated as a Partially Observable Stochastic Game (POSG)
[4, 5]. Deep Reinforcement Learning (DRL) methods [6] explicitly formulate the
POSG reward to derive a policy from simulated self-play experiments. While
such techniques are promising in sparse-data domains where the robot is easily
distinguishable from human, tuning reward parameters for homogeneous
navigation among humans is not trivial [7]. DRL policies are also a function
of the underlying simulator, which often translates to undesirable behavior
with the sim-to-real transfer owing to the lack of compatibility with real-
world scenarios.

On the other hand, data-driven approaches are common in social trajectory
prediction. They aim to directly characterize human interactions observed in
the data [2], eliminating the need for reward shaping and accurate
simulations. Recent sequence-to-sequence models, for instance, have achieved
promising results in intent prediction [8, 9, 10]. However, using these models
for downstream navigation is difficult as they often suffer from prediction
failures [11] which hurt their generalization capabilities. This potential for
unsafe behavior prompts the need for robustifying models deployed in the real
world.

Figure 1: The flight simulator setup used in our user study where, in randomly
selected order, each User interacted with our Human pilot, SoRTS and the
Baseline. The figure also shows a resulting interaction between a User and
SoRTS.

In social navigation, the actions of one agent influence those of another and
vice-versa [4, 12]. This temporally recursive decision-making intuition has
been used for modeling human-like gameplay [13, 14]. Leveraging this insight,
we propose using a recursive search-based policy to robustify offline-trained
models with downstream social navigation objectives. Specifically, we use
Monte Carlo Tree Search (MCTS) [15] as our search policy which provides long-
horizon recursive simulations, collision checking and goal conditioning. We
combine it with a socially-aware intent prediction model to provide short-
horizon agent-to-agent context cues and motion naturalness. We use MCTS to
fuse these short-horizon cues with long-horizon planning by including global
reference paths to guide the tree expansion. We refer to our framework as
Social Robot Tree Search (SoRTS).

The growing operations of Unmanned Aerial Vehicles are leading to a demand for
using airspace concurrently with human pilots [16, 17]. We, therefore, select
the domain of general aviation (GA) to showcase our approach. GA was recently
framed within the paradigm of social navigation [9, 10], where pilots are
expected to follow flying guidelines to coordinate with each other and respect
other’s personal space to ensure safe operations. This is analogous to
following etiquette in human crowds and vehicular settings.

Safety-critical domains, like GA, demand high competence to guarantee seamless
and safe operations. This entails developing trustworthy robots which are able
to understand and follow navigation norms but also understand long-term
dynamic interactions to avoid causing danger or discomfort to others. In this
paper, we separate these aspects into two axes, navigation performance and
safety. We center our framework design and evaluations around these two axes.
Through a user study conducted on a custom simulator framework, X-PlaneROS,
with 26 experienced pilots, we investigate how pilots interact with our model
in a realistic flight setting. We further analyze how they gauge our model
along these axes as compared to a competent human pilot. We also complement
our evaluations with self-play experiments in scenarios with increasing
complexity.

Statement of contributions:

  1. 1.

We introduce and open-source SoRTS, an MCTS-based algorithm for long-horizon
navigation that robustifies offline learned socially-aware intent prediction
policies for downstream navigation.

  2. 2.

We introduce X-PlaneROS, a high-fidelity simulation environment for navigation
in shared aerial space, and;

  3. 3.

Through a user study with 26 FAA-certified pilots and through self-play
simulations in more complex scenarios, we show that SoRTS is perceived
comparably to a competent human pilot in terms of navigation performance and
safety and significantly outperforms its baseline algorithm.

##  II Related Work

###  II-A Social Navigation Algorithms

Social navigation has a rich body of work focused on pedestrian and autonomous
driving domains [1]. In pedestrian settings, classical model-based approaches
[18, 19] have been proposed and remain prominent baselines. Yet, their
extension to other domains is non-trivial. Recent DRL methods [6, 20, 21] that
use handcrafted safety-focused reward functions [22] have produced promising
results in these domains. However, shortcomings in simulator design [23], and
domain-specific reward functions limit real-world performance [22]. Achieving
scalability and robustness is challenging, often requiring expensive
retraining. Similar to our approach, [5] introduces a DRL method that uses
MCTS to train and deploy policies. While their method relies on pre-defined
reward functions and simulator training, our work extends these methods to use
offline expert-based policies, providing domain-specific treatment for social
navigation.

Data-driven approaches focus on learning policies from datasets that record
interactions between agents [22, 1]. As these models do not need explicit
reward construction, they can capture the rich, joint social dynamics.
However, these methods are challenging to deploy directly owing to noisy and
missing demonstrations [24, 25]. To alleviate this, [4] used the gradients of
a Q-value function for Model Predictive Control, and [12] proposed a
generalization to this method using dual control for belief state propagation.
These methods rely on Inverse Reinforcement Learning as an additional step to
generate the Q-value functions. Using gradients from sequence models directly
in optimizations has also been proposed [26], but the convergence properties
were not examined. Our method is more direct in its use of sequential models
and calculating gradients or Q-values is not required. Instead, we transform
the model’s outputs into action distributions for the downstream planning
task.

###  II-B Social Navigation Evaluation

Different metrics have been considered for the evaluation of social robot
navigation [1, 27, 2]. Some of the main axes of analysis for evaluation
include behavior naturalness based on a reference trajectory or irregularity
of the executed path [28, 29], performance and efficiency [30, 31, 32], and
notions of physical personal space or discomfort [33, 34]. User studies are
often conducted to evaluate more subjective aspects such as the perceived
discomfort and trust that a robot induces during an interaction [35].
Following prior works, we focus on navigation performance to measure our
agent’s smoothness and ability to follow navigation guidelines, and safety to
judge its ability to respect others’ personal space. We also conduct a user
study where we ask experienced pilots to interact with our algorithm in a
realistic flight setting and rate the robot’s performance, perceived safety
and trust.

Figure 2: An overview of SoRTS, a Monte Carlo Tree Search (MCTS)-based planner
for social robot navigation whose tree search is guided by three components; a
Social Module, a Reference Module and a Cost Map. The Social Module leverages
an offline-trained intent prediction model to characterize agent-to-agent
interactions and predict a set of possible future states for the ego-agent.
The Reference Module provides the agent with the closest reference state from
a global plan which embodies navigation guidelines the agent must follow. The
Cost Map encodes a global visitation map to further guide navigation. MCTS
uses all these modules to provide online collision checking and long-horizon
simulations by searching for choosing between different decision modalities.

##  III Problem Formulation

We formulate social navigation as an approximate POSG, a framework for
decentralized finite-horizon planning. For more details about POSGs and their
use within social navigation, we refer the reader to [36, 5, 4].Following [10,
9], we assume M agents with ${\mathbf{s}}^{i}_{t}\in R^{3}$ representing state
of agent $i$ at time-step $t$ and ${\mathbf{a}}^{i}_{t}\in\mathcal{A}$ is the
discrete set of actions or motion primitives that follow
$\dot{{\mathbf{s}}}=f({\mathbf{s}},{\mathbf{a}})$. Let ${\mathbf{s}}^{i}_{0}$
and ${\mathbf{s}}^{i}_{g}$ represent the start and goal states respectively.
The system also has access to a set of offline expert demonstrations
$\mathcal{D}$ and set of global reference trajectories $\tau_{R}$. We omit
superscripts to refer to the joint state for all agents.

Thus, given the set of start ${\mathbf{s}}_{0}$ and goal ${\mathbf{s}}_{g}$
states of M agents, the objective is to find a sequence of control inputs
$\pi=\\{{\mathbf{a}}_{0},\dots,{\mathbf{a}}_{g}\\}$ such that the agents
follow collision-free trajectories
$\tau=\\{{\mathbf{s}},\cdots,{\mathbf{s}}_{g}\\}$. The generated trajectories
need to ensure $||{\mathbf{s}}^{i}_{t}-{\mathbf{s}}^{j}_{t}||\geq d~{}\forall
i,j\in\\{1,\dots,M\\}~{}\forall t$ where d is the minimum separation distance
to satisfy the safety objective. Furthermore, the trajectories also need to
stay close to reference trajectories $\min||\tau-\tau_{R}||$ to satisfy the
navigation objective and follow $\tau\sim\mathcal{D}$ to satisfy the social
objective. Without loss of generality, we assume first agent (i=0) to be the
robot ego-agent and at each time step execute the optimal action.

##  IV Approach

We designed SoRTS along the axes of navigation performance and safety. The
core insight of SoRTS is to use social modules and global reference paths to
bias the MCTS to search for and choose between various decision modalities.
Fig. 2 shows an example of two aircraft merging on a single path. Each
aircraft can safely execute the merger by choosing from three actions:
forward, speeding up, and slowing down. MCTS, in its forward simulations, not
only prunes branches that lead to future collisions but also uses the social
module to choose between cutting-in-front (socially undesirable) and yielding
(socially desirable), thereby producing socially compliant and safe behavior.

###  IV-A Modules

SoRTS is a search-based planner which uses Monte Carlo Tree Search (MCTS)
whose tree expansion is guided by three modules. First, a Social Module
handles the short-horizon dynamics in the scene, characterizing social
interactions and cues. Second, a Reference Module provides the agent with a
global navigation guideline, e.g., an airport traffic pattern. Finally, a Cost
Map encodes global value map. MCTS uses these components together and provides
collision checking and long-horizon socially-compatible simulations. Its
corresponding pseudo-code is shown in Algorithm 1 and 2.

####  IV-A1 Social Module

We leverage an offline-trained intent prediction algorithm parameterized in
$\theta$, $P_{\theta}(\cdot)$, to account for the short-term agent-to-agent
dynamics following the expert trajectories in $\mathcal{D}$.

${p}_{\theta}({\mathbf{s}}_{t}^{i},{\mathbf{a}}_{t}^{i})\sim
P_{\theta}({{{\mathbf{a}}}}_{t}^{i}\mid{{{\mathbf{s}}}}_{t-t_{obs}:t},{\mathbf{s}}_{g}^{i})$
(1)

where $p_{\theta}(\cdot)$ provides a distribution of future actions
${{{\mathbf{a}}}_{t}^{i}}$ for agent i conditioned on the past trajectories
${{{\mathbf{s}}}}_{t-t_{obs}:t}$ of all the agents and the goal
${{\mathbf{s}}_{g}^{i}}$ where $t_{obs}$ is the observation time horizon.

Here, we use Social-PatteRNN [10], an algorithm which predicts multi-future
trajectories from learned interactions that exploit motion pattern information
in the data.

####  IV-A2 Reference Path Module

Given the start and goal state of agent i the algorithm samples a suitable
reference trajectory $\tau_{r}^{i}\in\tau_{R}$, 1 in Algorithm 1. Similar to
section IV-A1, this trajectory is used to compute a reference action
distribution $p_{r}(\cdot)$,

${p}_{r}({\mathbf{s}}_{t}^{i},{\mathbf{a}}_{t}^{i})\sim
P_{r}({\mathbf{a}}_{t}^{i}|{\mathbf{s}}_{t}^{i},\tau_{r}^{i})$ (2)

proportional to the L2 norm between ${\mathbf{s}}$ and $\tau_{r}$ at time $t$.
In Algorithm 2, the reference action is obtained in 5. $\tau_{R}$ can be drawn
from expert distributions $\mathcal{D}$, global path planning algorithms like
A-Star, logic specifications [37] or can also be handcrafted.

####  IV-A3 Cost Map

The algorithm also uses a cost map of the environment representing the value
function, $v({\mathbf{s}})$. For our case, we use the cost map to represent
state visitation frequency to bias the search towards more desirable areas.
Algorithm 2 uses the cost map in 3. This value function can be either learned,
e.g., via self-play [5] or pre-computed from a prior distribution and captures
the value of the joint state distribution.

###  IV-B Social Monte Carlo Tree Search

MCTS is a search-based algorithm that expands its tree search toward high
rewarding trajectories. In principle, this is done by selecting nodes along
the search that maximize upper confidence bound [15], which balances the
degree of exploitation and exploration. SoRTS, leverages MCTS and uses the
components presented in the previous section to guide its trajectory roll-
outs. Formally, it follows the UCT shown below,

$\displaystyle U({\mathbf{s}},{\mathbf{a}})$
$\displaystyle=Q({\mathbf{s}},{\mathbf{a}})+c_{1}\cdot
S({\mathbf{s}},{\mathbf{a}})+c_{2}\cdot R({\mathbf{s}},{\mathbf{a}})$ (3)

where $Q({\mathbf{s}},{\mathbf{a}})$ represents the expected value for taking
action ${\mathbf{a}}$ at state ${\mathbf{s}}$; $S({\mathbf{s}},{\mathbf{a}})$
is visitation normalized component according to the socially-aware module;
$R({\mathbf{s}},{\mathbf{a}})$ is the visitation normalized component
according to the reference path; $c_{1}$ and $c_{2}$ are hyperparameters. We
drop the time $t$ subscript for ease of notation. In 14 of algorithm 2, these
values follow the update:

$\displaystyle Q({\mathbf{s}},{\mathbf{a}})$
$\displaystyle=\frac{N({\mathbf{s}},{\mathbf{a}})\cdot
Q({\mathbf{s}},{\mathbf{a}})+v({\mathbf{s}})}{N({\mathbf{s}},{\mathbf{a}})+1}$
(4) $\displaystyle R({\mathbf{s}},{\mathbf{a}})$
$\displaystyle=\frac{N({\mathbf{s}},{\mathbf{a}})\cdot
R({\mathbf{s}},{\mathbf{a}})+p_{r}({\mathbf{s}},{\mathbf{a}})}{N({\mathbf{s}},{\mathbf{a}})+1}$
(5) $\displaystyle S({\mathbf{s}},{\mathbf{a}})$
$\displaystyle=\frac{\sqrt{N({\mathbf{s}})}}{N({\mathbf{s}},{\mathbf{a}})+1}\cdot
p_{\theta}({\mathbf{s}},{\mathbf{a}})$ (6)

where $N({\mathbf{s}})$ is the state visitation count, and
$N({\mathbf{s}},{\mathbf{a}})$ the visitation given action ${\mathbf{a}}$.
These updates are done iteratively by following the states within a time-
budget, or until a new state is found. At each time-step, a new forward
simulation tree is iteratively constructed by alternately expanding the
agents’ future states in a round-robin fashion. Branches that lead to
collision states are pruned.333Note: In practice, for $M>2$, we only use the
ego-agent and the closest agent to the ego agent for tree expansion. While the
tree is explicitly constructed only for two agents, $p_{\theta}$ provides the
high-level social context for all the agents. This approximation preserves the
real-time nature of the algorithm and is shown to perform well in practice. At
the end of $planHorizon$, the ego agent takes the first action that maximizes
the visitation count in 4 of Algorithm 1. The tree is reset and the process
continues till goal is reached.

Algorithm 1 SoRTS(${\mathbf{s}}_{0},{\mathbf{s}}_{g},\theta,\tau_{R},v$)

1:$\tau_{r}\leftarrow$
GetReferencePaths$({\mathbf{s}}_{0},{\mathbf{s}}_{g},\tau_{R})$

2:while ${\mathbf{s}}_{t}^{0}\neq{\mathbf{s}}^{0}_{g}$ or $timeElapsed\leq
planHorizon$ do

3: $N(\cdot)\leftarrow$
SocialMCTS$({\mathbf{s}}_{t},{\mathbf{s}}_{G},\tau_{r},\theta,v,0)$

4:
${\mathbf{a}}^{0}\leftarrow\arg\max_{\mathbf{a^{\prime}}\in\mathcal{A}}N({\mathbf{s}}_{t}^{0},\mathbf{a^{\prime}})$

5: ${\mathbf{s}}_{t+1}^{0}\leftarrow f({\mathbf{s}}_{t}^{0},{\mathbf{a}}^{0})$
$\triangleright$ Robot’s transition model.

6:end while

Algorithm 2
SocialMCTS(${\mathbf{s}}_{t},{\mathbf{s}}_{G},\tau_{r},\theta,v,p$)

1:while p $\leq$ M do

2: if ${\mathbf{s}}_{t}^{p}\not\in S(\cdot)$ then $\triangleright$ New tree
node

3: $v_{s}\leftarrow$ GetValue$(v,{\mathbf{s}}_{t})$

4: $S(\cdot)\leftarrow$
GetSocialActionProbabilities$({\mathbf{s}}_{t},\theta)$

5: $p_{r}(\cdot)\leftarrow$
GetReferenceActionProbabilities$(\tau^{p}_{r},{\mathbf{s}}_{t}^{p})$

6: $N({\mathbf{s}}_{t}^{p})\leftarrow 1$

7: return $v_{s}$

8: end if

9: $\mathcal{A^{\prime}}\leftarrow$
CollisionCheck$(\mathcal{A},d,{\mathbf{s}}_{t})$

10:
${\mathbf{a}}^{p}\leftarrow\arg\max_{\mathbf{a^{\prime}}\in\mathcal{A^{\prime}}}U({\mathbf{s}}_{t}^{p},\mathbf{a^{\prime}})$
$\triangleright$ See eq. 3

11: ${\mathbf{s}}_{t+1}^{p}\leftarrow
f({\mathbf{s}}_{t}^{p},{\mathbf{a}}^{p})$

12: p $\leftarrow$ p+1 $\triangleright$ Choose next agent to expand

13: SocialMCTS$({\mathbf{s}},{\mathbf{s}}_{g},\tau_{r},\theta$, v, p)

14: Update$(Q,R,S,N)$ $\triangleright$ See eq. 4 to 6

15:end while

16:return $N(\cdot)$

##  V Experimental Setup

Our experiments focus on assessing SoRTS along our axes of interest, i.e.,
navigation performance and safety. As such, this section describes the main
aspects of our evaluation setup and implementation details for our case study.

###  V-A User Study

We recruited 26 FAA-certified pilots (14 private, 8 commercial, 3 student
pilots and 1 airline transport pilot), who on average have 986 flight hours.
Using the flight deck setup and simulator shown in fig. 1, each pilot was
asked to land an aircraft on a specified runway. Simultaneously, a second
pilot was solving the same task, thus, requiring the pilot’s coordination to
land safely. Here, the second pilot was either a human, SoRTS, or the baseline
algorithm444We use second pilot and algorithm interchangeably..

Our experimental setup follows the design delineated next. We followed a
within-subject design in which each user tests against each algorithm. We
allowed the user to get familiarized with the simulator and controls prior to
the tests. The pilots are spawned at a 10 km radius from the airport, where
their incoming direction is either north (N), south (S), or west (W), defining
six possible scenarios, i.e., $\\{(N,S),(S,N),(N,W),(W,N),(S,W),(W,S)\\}$. The
algorithm order and scenario were selected randomly at the beginning of the
experiment. The scenario, the initial states and final goals remained fixed
throughout the three tests.

After each test, the user completed a 5-scale Likert questionnaire evaluating
the second pilot along various factors we deemed relevant for high navigation
performance and safety. Since operation in safety-critical domains demands
high competence, we also take an interest in understanding which aspects the
users prioritize for deeming a pilot as trustworthy and competent. We, thus,
ask users to also rate how trustworthy and human-like they perceived the
second pilot, based on their interaction with them and their responses along
the axes of navigation and safety. We summarize the components of our user
study questionnaire in table I. Finally, we also collected the trajectory data
from each experiment for further analysis using relevant metrics discussed in
section V-C.

###  V-B Self-Play

We complement our user study with self-play simulations to assess SoRTS and
our baseline on a wider variety of scenarios. The simulations follow a similar
design as the user study, but now allow each agent’s location to be anywhere
around the 10km radius to increase the number of possible scenarios. We also
consider multi-agent scenarios varying from 2 to 4 agents. We randomly
generate 100 episodes for each setting, where in each episode an agent is
deemed unsuccessful if it breaches a minimum separation distance, or if it
reaches a maximum number of allowed steps.

###  V-C Metrics

To quantify the trajectory data collected in our user study and self-play
experiments along our axes of analysis, we consider the metrics listed below;

  1. 1.

Reference Error (RE): The Euclidean distance between the agent’s reference
trajectory and its executed path.

  2. 2.

Loss of Separation (LS): The duration that two agents break a minimum distance
of each other. This metric is relevant for the aviation domain [38], but is
akin to commonly used metrics for social robots in pedestrian settings, e.g.,
personal space [1, 27].

###  V-D Baseline

To showcase the benefits of SoRTS, we introduce a baseline which naively
chooses the next action by balancing the reference and social values over the
state-action space following the equation below,

$\displaystyle{\mathbf{a}}=\arg\max_{{\mathbf{a}}^{{}^{\prime}}\in\mathcal{A}}\big{[}\lambda\cdot
p_{r}({\mathbf{s}},{\mathbf{a}})+(1-\lambda)\cdot
p_{\theta}({\mathbf{s}},{\mathbf{a}})\big{]}$

where $\lambda\in{\mathbb{R}}:[0,1]$ balances the importance we give to
$p_{\theta}(\cdot)$ and $p_{r}(\cdot)$. This baseline translates to replacing
3 and 4 for the above equation in algorithm 1.

###  V-E Simulator

To evaluate the proposed algorithm and enable future research in the domain of
full-scale aerial autonomy, we introduce X-PlaneROS. The system combines two
main components via a ROS bridge, X-Plane-11 and ROS-Plane autopilot [39],
enabling the use of high and low-level commands to control a general aviation
aircraft in realistic world scenarios. X-Plane-11 is a high-fidelity simulator
which provides realistic aircraft models and visuals. ROS-Plane is a widely
accepted tool which provides reliable autonomous flight control loops. Our
simulator further adds support for following a select set of motion
primitives, as well as visualization utilities that aid in tuning the control
loops. The documentation and source code is open-sourced555Code for
X-PlaneROS: https://github.com/castacks/xplane_ros.

###  V-F Implementation Details

The modules in section IV-A leverage TrajAir [9], a dataset consisting of
aircraft trajectory data collected in non-towered terminal airspace; the
Social Module uses TrajAir to train its intent prediction module offline,
where we followed the training details in [10]. We also use TrajAir to build a
library of FAA-abiding paths used by the Reference Module. Finally, to build
our cost map we discretized TrajAir’s flight frequency based on aircraft poses
and wind direction.

The action space, $\mathcal{A}$, of our planner is comprised of the set of 252
motion primitives for aerial navigation used in [37]. Empirically, we set
$planHorizon$=10s in algorithm 2 and the exploration parameters in the UCT
equation to $c_{1}=2$ and $c_{2}=5$, and $\lambda=0.3$ for the baseline
planner.

TABLE I: We asked users to rate the second pilot along these aspects of
navigation performance and safety, and to rate it along trustworthiness and
humanness.

Navigation Performance | Safety | Was the second pilot…
---|---|---
Follow FAA Guidelines | Collision Risk | trustworthy?
Overall Flying Skill | Comfort | a human?
Flight Smoothness | Cooperative |
| Abrupt |
| Predictable |

##  VI Results

Figure 3: User study results. Left: Examples of two experiments. Each row
shows the resulting trajectories of a User interacting with our Human pilot,
SoRTS, and the Baseline. The top row shows a user that successfully followed
the expected guideline; we can see that the baseline did not fly as smoothly
as SoRTS, and abruptly cut short close to the goal. The bottom row shows a
user that did not follow the expected guideline; our algorithm still managed
to navigate properly, whereas the baseline behaved erratically, unsafely
crossing over the runway. Right: Box-plots showing the per-algorithm
distribution of the results from the user study: the top ones show the average
scores given by the user for navigation performance (a) and safety (b). The
bottom ones show the average reference error (a) and loss-of-separation (b)
metrics obtained from the trajectory data.

We now present our results and insights. For our analysis below, a competent
algorithm refers to one which performs highly along the axes of navigation
performance and safety. We first examine the relevance of the factors in table
I in characterizing the trustworthiness of a competent algorithm. Leveraging
the results from our trustworthiness analysis, we then perform a comparative
analysis between SoRTS, the baseline and our human pilot. Finally, we briefly
explore how pilots perceive human competence along the factors in our
questionnaire.

###  VI-A On competence and trustworthiness

As a means to determine how to better assess the competence of a robot, e.g.,
which metrics to prioritize, we compare the factors in table I with the user’s
perceived trust. We do so via a Pearson’s correlation analysis with repeated
measures, whose results are summarized in table II.

We find strong correlations between trust and the components within navigation
performance, with flight smoothness being the strongest one. Within safety, we
see that comfort, followed by predictability showed the strongest
correlations. Although one would expect for abruptness and collision risk to
be as crucial for rating trust. We find this result is, in part, due to our
experiments not being explicitly designed to exhibit adversarial behavior. We
also believe these results showcase the strength of our model in unprecedented
scenarios. As said earlier, because data-driven models can heavily misbehave
when exposed to out-of-distribution states [11], SoRTS is a means to augment
these models online. Since, the user study serves as a mechanism for testing
our model in settings that weren’t observed in the training data, the lower
correlation between the aforesaid components may indicate robustness toward
unsafe behavior.

TABLE II: Using repeated measures, we show the Pearson correlations (R,
p-value=0.05) for each factor listed in table I vs. perceived trustworthiness
and humanness.

Axis | Factor | Trustworthy? | Human?
---|---|---|---
R | p-value | R | p-value
Nav. Performance | Flight Smoothness | 0.81 | 4.80e-19 | 0.08 | 0.47
Follow FAA Guidelines | 0.76 | 2.99e-15 | 0.10 | 0.37
Overall Flying Skill | 0.71 | 4.56e-13 | 0.15 | 0.20
Safety | Comfortable | 0.92 | 3.55e-31 | 0.17 | 0.14
Predictable Behavior | 0.77 | 5.19e-16 | 0.21 | 0.07
Cooperative | 0.65 | 1.61e-10 | 0.06 | 0.62
Collision Risk | -0.56 | 1.10e-07 | -0.11 | 0.36
Abrupt | -0.56 | 1.13e-07 | -0.09 | 0.42

###  VI-B On each algorithm’s competence

The previous section studied how each factor in table I tied to the perception
of a trustworthy pilot. This section will now focus on providing a pairwise
comparison between the algorithms for both the user study and self-play,
leveraging the results obtained in the previous analysis.

####  VI-B1 User study

Leveraging the trustworthiness results in table II, we compute each user’s
average score over the questions with highest correlations along each axis,
i.e., for navigation performance, we compute the mean score between following
FAA and flight smoothness, and for safety we use predictability and comfort.
The scores are shown in fig. 3, and the pairwise statistical differences using
ANOVA with repeated measures in table III. Here, we did not find statistical
evidence that the scores for the human pilot and SoRTS were different,
suggesting that the users rated their competence similarly. Our results
further show evidence that our baseline was generally rated lower on both axes
of competence, while also exhibits more variance.

We then examine if the RE and LS metrics in section V-C—which here tie to
navigation performance and safety, respectively—reflect the users’ assessments
for each algorithm. Thus, we compute them on the collected trajectories and
show their respective mean scores and statistical comparisons in fig. 3 and
table III. As before, we show the mean scores for each algorithm in fig. 3 and
the statistical results in table III. For the RE metric, we find significantly
different results between the algorithms; wherein SoRTS yields higher error
compared the human pilot, but markedly lower than the baseline and less
variance than the other two. Similarly, the LS metric computed at 0.5km, shows
that generally neither the human pilot nor SoRTS breach this distance. In
contrast, the baseline more frequently invades others’ personal space,
creating more situations for potential collisions.

TABLE III: Statistical significance between algorithmic pairs for results in
fig. 3 with $t^{\ast}\geq 2.060$ and $p\leq 0.05$.

Algorithmic | NP | Safety | RE | LS
---|---|---|---|---
Pair | t-val | p-val | t-val | p-val | t-val | p-val | t-val | p-val
Baseline-Human | 3.121 | 0.009 | 3.062 | 0.016 | 5.782 | 0.000 | 1.321 |  0.199
Baseline-SoRTS | 3.018 | 0.009 | 2.626 | 0.022 | 2.105 | 0.011 | 0.397 |  0.694
Human-SoRTS | 0.415 | 0.682 | 0.322 | 0.322 | 2.944 | 0.022 | 1.211 |  0.237

  * •

NP: Navigation Performance, LS: Loss of Separation, RE: reference error.

Figure 3 shows examples of trajectories from our experiments. Each row
represents one user vs. the three algorithms. We show a reference trajectory
along with the actual executed trajectory666Human pilots do not see the
reference trajectory, but they see the map in the simulator, and are
instructed to follow the appropriate traffic patterns.. The top row shows an
example where the user performed successfully. In this experiment, we see that
the ablation unexpectedly cut short while approaching the runway instead of
following the traffic pattern, while SoRTS did so smoothly and safely. The
bottom row shows a user that did not follow the standard traffic pattern.
Despite this, SoRTS manages to successfully complete the task, while the
baseline behaves erratically, not following the traffic pattern and unsafely
traversing the runway twice.

As per the user study, we conclude that SoRTS performs comparable to a
competent human pilot and significantly better than the ablation, both as
perceived by the users and through our metrics. Our results further highlight
the benefits of using the long-horizon socially-aware simulations via MCTS, as
opposed to purely using a data-driven model on the wild with a simple
weighting over the objectives of interest, leading to an erratic and unsafe
behavior.

####  VI-B2 Self-play

TABLE IV: Summary of the task performance of the baseline and SoRTS’ agents in
self-play.

# Agents | Algorithm | Success (%) | Failure (%) | RE
---|---|---|---|---
LS | Timeout
2 | Baseline | 69.0 | 30.0 | 1.0 | 1.59
SoRTS | 96.5 | 2.0 | 1.5 | 2.01
3 | Baseline | 48.3 | 49.3 | 2.4 | 1.55
SoRTS | 89.7 | 8.3 | 2.0 | 2.01
4 | Baseline | 43.5 | 54.0 | 2.5 | 1.56
SoRTS | 71.0 | 21.0 | 8.0 | 1.96

  * •

LS: Loss of Separation at 0.5km, RE: reference error for successful agents.

We now analyze the performance of our algorithm in more complex scenarios. The
results summarized in Table IV show the percentage of agents that were able to
land on the runway, i.e., successful agents, the agents that were
unsuccessful, due to either a loss of separation, or because they exceeded the
maximum allowed time. Finally, we show the average reference error for
successful agents.

Analyzing the table, we observe a decrease in task success as the number of
agents increases. Nonetheless, we can see that SoRTS performs significantly
better than the ablation algorithm improving the task success by 28%, 46% and
38% for the 2, 3 and 4 agent scenarios, respectively. We also find that, the
average reference error for SoRTS is higher than the ablation. Though this
error was computed on successful agents only, we hypothesize that it being
higher is due to conflict resolution to avoid collisions with other agents.

###  VI-C On competence and human performance

Often arising debate, is the question of whether an algorithm can pass as a
human. However, in domains such as aviation, where high competence is central,
one would expect for humans to be perceived as such. If so was the case, one
would strive for an algorithm to exhibit a performance similar to that of a
human. We study the aforesaid by asking users to gauge whether the second
pilot was a human based on their assessments along table I.

Surprisingly, we find a marked disagreement within the users’ responses to
this question, with almost a 50-50% split between the responses for the Human
Pilot (14: No, 12: Yes) and the Baseline (14: No, 12: Yes), whereas for SoRTS
(8: No, 18: Yes) more users perceived its performance as human-like. Further,
table II shows that the user’s assessments for navigational performance and
safety correlate weakly with the humanness prediction. To explain this
contrasting result, we isolated the responses given to SoRTS from the other
two, and found that predictability had the highest correlation value among all
factors, with R=0.53 and p=0.01. In contrast, the values for the Human pilot
and the Baseline were R=0.22, p=0.03 and R=-0.06, p=0.06, respectively. We
believe this result may suggest that because users generally perceived the
behavior of SoRTS as more predictable, they rated its performance as human-
like more frequently.

We believe this contrasting result does not affect our analyses from the
previous to sections, as we find that the user’s responses for trustworthiness
also correlate weakly with their assessment of humanness, with R=0.26, p=0.02.
As such, we conclude that users value more their perception of trust than the
nature of type of agent they interact with.

##  VII Conclusion

We present SoRTS, an algorithm for long-horizon social robot navigation. SoRTS
is a MCTS-based planner which expands its search tree guided by an offline-
trained intent prediction model and a global path which embodies navigation
guidelines. We introduce X-PlaneROS, a high fidelity simulator for research in
full-scale aerial autonomy. We use it to conduct a user study with experienced
pilots to study our algorithm’s performance in realistic flight settings. We
find that users perceive SoRTS comparable to a competent human pilot and
significantly better than our baseline. In self-play, we show that SoRTS
outperforms the baseline by 28-46% as the complexity of scenarios increases.
To the best of our knowledge, this is the first work in social navigation for
general aviation and attempts to bring unique problems in general aviation
within the purview of the larger robotics community.

We identify various avenues for future work. First, we assumed task
homogeneity, i.e., agents landing on the same runway, whereas in a real
scenario pilots with different objectives may require to interact. Thus,
future work includes studying interactions with agents with heterogeneous
tasks. We also assumed perfect intent and state estimation. Accordingly,
robustifying prediction models with uncertainty- and adversarial-awareness is
another promising direction. Finally, improvements on the scalability of the
model, as well as exploring other domains are also potential avenues.

## ACKNOWLEDGMENT

This work was supported by the Army Futures Command Artificial Intelligence
Integration Center (AI2C), the Ministry of Trade, Industry and Energy (MOTIE),
the Korea Institute of Advancement of Technology (KIAT), and the Brazilian Air
Force. The views expressed in this article do not necessarily represent those
of the aforementioned entities. We thank Jasmine Jerry Aloor for her support.
We also thank the pilots from Condor Aero Club (KPJC) and ABC Flying Club
(KAGC) for their participation in our user study.

## References

  * [1] C. Mavrogiannis, F. Baldini, A. Wang, D. Zhao, P. Trautman, A. Steinfeld, and J. Oh, “Core challenges of social robot navigation: A survey,” _arXiv preprint arXiv:2103.05668_ , 2021.
  * [2] A. Rudenko, L. Palmieri, M. Herman, K. M. Kitani, D. M. Gavrila, and K. O. Arras, “Human motion trajectory prediction: A survey,” _The International Journal of Robotics Research_ , vol. 39, no. 8, pp. 895–935, 2020\.
  * [3] R. Tian, L. Sun, A. Bajcsy, M. Tomizuka, and A. D. Dragan, “Safety assurances for human-robot interaction via confidence-aware game-theoretic human models,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 11 229–11 235.
  * [4] D. Sadigh, N. Landolfi, S. S. Sastry, S. A. Seshia, and A. D. Dragan, “Planning for cars that coordinate with people: leveraging effects on human actions for planning and active information gathering over human internal state,” _Autonomous Robots_ , vol. 42, no. 7, pp. 1405–1426, 2018.
  * [5] B. Riviere, W. Hönig, M. Anderson, and S.-J. Chung, “Neural tree expansion for multi-robot planning in non-cooperative environments,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 4, pp. 6868–6875, 2021.
  * [6] S. Matsuzaki and Y. Hasegawa, “Learning crowd-aware robot navigation from challenging environments via distributed deep reinforcement learning,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 4730–4736.
  * [7] C. I. Mavrogiannis, V. Blukis, and R. A. Knepper, “Socially competent navigation planning by deep learning of multi-agent path topologies,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2017, pp. 6817–6824.
  * [8] T. Salzmann, B. Ivanovic, P. Chakravarty, and M. Pavone, “Trajectron++: Dynamically-feasible trajectory forecasting with heterogeneous data,” in _European Conference on Computer Vision_. Springer, 2020, pp. 683–700.
  * [9] J. Patrikar, B. Moon, J. Oh, and S. Scherer, “Predicting like a pilot: Dataset and method to predict socially-aware aircraft trajectories in non-towered terminal airspace,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 2525–2531.
  * [10] I. Navarro and J. Oh, “Social-patternn: Socially-aware trajectory prediction guided by motion patterns,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2022, pp. 9859–9864.
  * [11] A. Farid, S. Veer, B. Ivanovic, K. Leung, and M. Pavone, “Task-relevant failure detection for trajectory predictors in autonomous vehicles,” in _Conference on Robot Learning_. PMLR, 2023, pp. 1959–1969.
  * [12] H. Hu and J. F. Fisac, “Active uncertainty reduction for human-robot interaction: An implicit dual control approach,” _arXiv preprint arXiv:2202.07720_ , 2022.
  * [13] D. Silver, T. Hubert, J. Schrittwieser, I. Antonoglou, M. Lai, A. Guez, M. Lanctot, L. Sifre, D. Kumaran, T. Graepel, _et al._ , “A general reinforcement learning algorithm that masters chess, shogi, and go through self-play,” _Science_ , vol. 362, no. 6419, pp. 1140–1144, 2018.
  * [14] A. P. Jacob, D. J. Wu, G. Farina, A. Lerer, H. Hu, A. Bakhtin, J. Andreas, and N. Brown, “Modeling strong and human-like gameplay with kl-regularized search,” 2021.
  * [15] L. Kocsis and C. Szepesvári, “Bandit based monte-carlo planning,” in _Machine Learning: ECML 17th European Conference on Machine Learning, Berlin, Germany, September 18-22, 2006, Proceedings_ , ser. Lecture Notes in Computer Science, J. Fürnkranz, T. Scheffer, and M. Spiliopoulou, Eds., vol. 4212. Springer, 2006, pp. 282–293.
  * [16] J.-P. Aurambout, K. Gkoumas, and B. Ciuffo, “Last mile delivery by drones: An estimation of viable market potential and access to citizens across european cities,” _European Transport Research Review_ , vol. 11, no. 1, pp. 1–21, 2019.
  * [17] M. Grote, A. Pilko, J. Scanlan, T. Cherrett, J. Dickinson, A. Smith, A. Oakey, and G. Marsden, “Sharing airspace with uncrewed aerial vehicles (uavs): Views of the general aviation (ga) community,” _Journal of Air Transport Management_ , vol. 102, p. 102218, 2022.
  * [18] J. v. d. Berg, S. J. Guy, M. Lin, and D. Manocha, “Reciprocal n-body collision avoidance,” in _Robotics research_. Springer, 2011, pp. 3–19.
  * [19] C. Mavrogiannis, P. Alves-Oliveira, W. Thomason, and R. A. Knepper, “Social momentum: Design and evaluation of a framework for socially competent robot navigation,” _ACM Transactions on Human-Robot Interaction (THRI)_ , vol. 11, no. 2, pp. 1–37, 2022.
  * [20] Y. F. Chen, M. Liu, M. Everett, and J. P. How, “Decentralized non-communicating multiagent collision avoidance with deep reinforcement learning,” in _2017 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2017, pp. 285–292.
  * [21] C. Chen, Y. Liu, S. Kreiss, and A. Alahi, “Crowd-robot interaction: Crowd-aware robot navigation with attention-based deep reinforcement learning,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 6015–6022.
  * [22] C.-E. Tsai and J. Oh, “A generative approach for socially compliant navigation,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 2160–2166.
  * [23] A. Biswas, A. Wang, G. Silvera, A. Steinfeld, and H. Admoni, “Socnavbench: A grounded simulation testing framework for evaluating social navigation,” _ACM Transactions on Human-Robot Interaction (THRI)_ , vol. 11, no. 3, pp. 1–24, 2022.
  * [24] M. A. Bashiri, B. Ziebart, and X. Zhang, “Distributionally robust imitation learning,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 24 404–24 417, 2021.
  * [25] F. Codevilla, E. Santana, A. M. López, and A. Gaidon, “Exploring the limitations of behavior cloning for autonomous driving,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2019, pp. 9329–9338.
  * [26] S. Schaefer, K. Leung, B. Ivanovic, and M. Pavone, “Leveraging neural network gradients within trajectory optimization for proactive human-robot interactions,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 9673–9679.
  * [27] Y. Gao and C.-M. Huang, “Evaluation of socially-aware robot navigation,” _Frontiers in Robotics and AI_ , p. 420, 2021.
  * [28] C. Mavrogiannis, A. M. Hutchinson, J. Macdonald, P. Alves-Oliveira, and R. A. Knepper, “Effects of distinct robot navigation strategies on human behavior in a crowded environment,” in _2019 14th ACM/IEEE International Conference on Human-Robot Interaction (HRI)_. IEEE, 2019, pp. 421–430.
  * [29] A. J. Sathyamoorthy, J. Liang, U. Patel, T. Guan, R. Chandra, and D. Manocha, “Densecavoid: Real-time navigation in dense crowds using anticipatory behaviors,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 11 345–11 352.
  * [30] S. Liu, P. Chang, W. Liang, N. Chakraborty, and K. Driggs-Campbell, “Decentralized structural-rnn for robot crowd navigation with deep reinforcement learning,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 3517–3524.
  * [31] M. Everett, Y. F. Chen, and J. P. How, “Motion planning among dynamic, decision-making agents with deep reinforcement learning,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 3052–3059.
  * [32] J. Liang, U. Patel, A. J. Sathyamoorthy, and D. Manocha, “Realtime collision avoidance for mobile robots in dense crowds using implicit multi-sensor fusion and deep reinforcement learning,” _arXiv preprint arXiv:2004.03089_ , 2020.
  * [33] E. Torta, R. H. Cuijpers, and J. F. Juola, “Design of a parametric model of personal space for robotic social navigation,” _International Journal of Social Robotics_ , vol. 5, no. 3, pp. 357–365, 2013.
  * [34] Y. F. Chen, M. Everett, M. Liu, and J. P. How, “Socially aware motion planning with deep reinforcement learning,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2017, pp. 1343–1350.
  * [35] J. Butler and A. Agah, “Psychological effects of behavior patterns of a mobile personal robot,” _Autonomous Robots_ , vol. 10, pp. 185–202, 03 2001.
  * [36] R. Emery-Montemerlo, G. J. Gordon, J. G. Schneider, and S. Thrun, “Approximate solutions for partially observable stochastic games with common payoffs,” in _3rd International Joint Conference on Autonomous Agents and Multiagent Systems (AAMAS 2004), 19-23 August 2004, New York, NY, USA_. IEEE Computer Society, 2004, pp. 136–143.
  * [37] J. J. Aloor, J. Patrikar, P. Kapoor, J. Oh, and S. Scherer, “Follow the rules: Online signal temporal logic tree search for guided imitation learning in stochastic domains,” _arXiv preprint arXiv:2209.13737_ , 2022.
  * [38] T. Glozman, A. Narkawicz, I. Kamon, F. Callari, and A. Navot, “A vision-based solution to estimating time to closest point of approach for sense and avoid,” in _AIAA Scitech 2021 Forum_ , 2021, p. 0450.
  * [39] G. Ellingson and T. McLain, “Rosplane: Fixed-wing autopilot for education and research,” in _2017 International Conference on Unmanned Aircraft Systems (ICUAS)_. IEEE, 2017, pp. 1503–1507.


capbtabboxtable[][] aainstitutetext: Department of Physics, Istanbul Technical
University,
34469 Maslak, Istanbul, Turkey bbinstitutetext: Department of Mathematics,
Istanbul Technical University,
34469 Maslak, Istanbul, Turkey

# Stability of the Cosmological Dynamics of $O(D,D)$-complete Stringy Gravity

A. Savaş Arapoğlu a, 1 Sermet Çağan111Corresponding author. b and Aybike
Çatal-Özer<EMAIL_ADDRESS>
###### Abstract

The massless fields in the universal NS-NS sector of string theory form
$O(D,D)$ multiplets of Double Field Theory, which is a theory that provides a
T-duality covariant formulation of supergravity, leading to a stringy
modification of General Relativity. In this framework, it is possible to write
down the extensions of the Einstein field equations and the Friedmann
equations in such a way that the coupling of gravitational and matter sectors
is dictated by the $O(D,D)$ symmetry universally. In this paper, we obtain the
autonomous form of the $O(D,D)$-complete Friedmann equations, find the
critical points and perform their stability analysis. We also include the
phase portraits of the system. Cosmologically interesting cases of scalar
field, radiation, and matter are separately considered and compared with the
Chameleon models in a similar setting. Accelerating phases and the conditions
for their existence are also given for such cases.

##  1 Introduction

Alternative theories of gravity to general relativity have become an active
field of research, especially in the last 30 years, for both observational and
conceptual reasons. When observations at cosmological scales are tried to be
understood within the framework of general relativity, the fact that unknown
matter-energy components such as dark matter and dark energy have to be
included in the models of the universe is one of the most important reasons
for this field to become active. Moreover, some tensions/incompatibilities of
the $\Lambda$CDM model built within the framework of general relativity, which
has recently emerged as a result of observations made at various scales, has
made it a kind of necessity to address the theoretical framework of the
cosmological standard model. Considering all these about the early universe
processes suggests that alternative theories of gravity should be examined
together with consistent quantum gravity theory candidates such as string
theory. The low energy limits of string theories provide an important
alternative set of gravitational theories for understanding how the early
universe processes and the present large-scale structure of the universe
emerged with such processes. Cosmology, thus, offers a very promising window
for examining the observational consequences of string theory and is also
easier in the sense that studying only time-dependent backgrounds is simpler
than considering space and time coordinates together.

‘Stringy gravity’ emerging from the low energy limit of string theory has more
fundamental gravitational fields than general relativity: in addition to the
metric $g_{\mu\nu}$, there are also other massless modes of the NS-NS sector,
which are the dilaton and the two-form field $B_{\mu\nu}$. The dilaton and
B-field are, in general, treated as moduli to be dynamically stabilized by
some means in the quest for superstring vacua damour-poly94-1 ; damour-
poly94-2 , and the dilaton is included in the framework of SUGRA cosmology
gasper2003 ; vafa1992 ; gasper2008 ; veneziano91 ; gasper92 . A crucial
difficulty with the inclusion of the dilaton and/or B-field in the search for
cosmological solutions is how to couple these stringy gravitational fields
with matter. The usual strategy is to minimally couple the metric to the
matter sector; but the couplings of dilaton and B-field to matter are omitted
since these couplings lead to some difficulties concerning observational
restrictions, like solar system constraints (damour-poly94-2, ). Actually,
coupling dilaton (or any light scalar) to the matter sector is highly
restricted by similar constraints coming from different scales. A way of
evading these problems caused by scalar fields is to introduce some
‘screening’ mechanism to build cosmological models consistent with
observations (brax2013, ). The coupling of the scalar field to matter,
together with the form of the scalar field potential may lead to a rich
phenomenology, including non-perturbative phenomena such as spontaneous
scalarization also damour1993 ; fethi2024  and references therein.

From this point of view, it would be crucial to determine how the scalar field
couples to matter, not just on a phenomenological ground, but in a way
dictated by basic symmetry principles of the underlying (quantum) theory of
gravity. Observational implications resulting from this coupling, imposed by
fundamental symmetry principles of the theory, would also be critical to test
the consistency of the theory. In string theory, the dilaton’s interaction
with matter is usually interpreted as an ad hoc addition, necessitated by the
conformal anomaly and specific field content (GSW). Double Field Theory (DFT)
changes this perspective by introducing an $O(D,D)$ T-duality covariant action
that reduces to the low energy action for the universal sector of string
theory, when a constraint, again $O(D,D)$ covariant, is imposed (hull-
zwiebach2009, ). This is achieved by introducing $D$ coordinates conjugate to
string winding modes in addition to the $D$ space-time coordinates conjugate
to momentum modes, in such a way that they form an $O(D,D)$ mutliplet. The
dynamics of the theory is governed by field equations, which involves a
generalized Ricci tensor and a generalized Ricci scalar. The generalized Ricci
tensor $\mathcal{R}(\mathcal{H},d)$ involves the generalized dilaton $d$ and
the generalized metric $\mathcal{H}$, which comprises of the Riemannian metric
and the antisymmetric B-field hohmetal2010-1 ; hohmetal2010-2 . DFT has also
been extended so as to incorporate matter fields and in this framework it is
possible to write the generalization of Einstein Field Equations by
generalizing not just the Einstein tensor but also the energy-momentum tensor.
This makes it possible to fix the coupling of matter to metric, the B-field
and, most importantly for the purposes of this paper, the dilaton field by
requiring the interactions to be consistent with the $O(D,D)$ T-duality
symmetry. The generalized equations describing the interaction between the
gravitational and the matter sectors are called the Einstein Double Field
Equations (EDFEs) (angus2018, ), and they may be considered as
$O(D,D)$-completion of equations of GR (For the earlier discussions on the
cosmological implications of the $O(D,D)$ symmetry include wu2014 ; ma2014 ;
branden2018-1 ; branden2018-2 ; brandenetal2019 ; branden-bernardo2019 ). The
generalized Friedmann equations following from this comprehensive
$O(D,D)$-symmetric framework, called the $O(D,D)$-complete Friedmann Equations
(OFEs), is derived in (Angus2020, ). This $O(D,D)$-symmetric formulation
modifies SUGRA cosmology in such a way that whenever the dilaton is
stabilized, the standard Friedmann equations are obtained not only for a
radiation-dominated universe as in the traditional string cosmology but also
in the presence of any matter sources.

Our primary motivation in this work is to examine the viability of the
coupling of dilaton to the matter sector dictated by the $O(D,D)$ symmetry of
the Double Field Theory (DFT) in the cosmological setting. In this study, we
compare the $O(D,D)$-extended system with the so-called “Chameleon Models”
where the form of the potential of the scalar field and the coupling of the
dilaton to matter sector may lead to a special type of a screening effect for
the scalar field. In the analysis of $O(D,D)$ system and comparison with
Chameleon model, we use the dynamical system approach (bohmer2018, ). This
allows to set up a framework that takes into account plausible cosmic
situations and explores the viability of the models (bernardo2020, ).

The organization of the paper is as follows. In section 2, $O(D,D)$-complete
Friedmann Equations are put into the form of a system of autonomous equations.
In section 3, the critical points are determined and their stability analysis
is considered, together with phase portraits. In section 4, some
cosmologically interesting cases for phenomenological applications are
considered and possible implications are mentioned. In section 5,
$O(D,D)$-complete cosmological models are compared with another string-
inspired phenomenological model, called Chameleon gravity chameleon ;
chameleon2 ; chameleon3 . We conclude with some comments and possible future
work in section 6.

##  2 Autonomous Equation System for $O(D,D)$-complete Cosmology

The $O(D,D)$ complete Friedmann equations was constructed in Ref. Angus2020
by imposing the most general ansatz for a homogeneous and isotropic
cosmological background in four dimensions. The ansatz is dictated by
isometries of DFT and substituting them in Einstein Double Field Equations
yields $O(D,D)$-complete Friedmann Equations as given in Angus2020 . Here, our
aim is to do linear stability analysis, so we rearrange these equations using
the cosmic gauge by setting the lapse function, $N(t)$, to unity. As a result,
we arrive at the following set of $O\left(D,D\right)$-complete Friedmann
equations,

$\displaystyle\frac{8\pi G}{3}\rho e^{2\phi}+\frac{h^{2}}{12a^{6}}$
$\displaystyle=H^{2}-2\dot{\phi}H+\frac{2}{3}\dot{\phi}^{2}$ (1a)
$\displaystyle\frac{4\pi
G}{3}\left(\rho+3p\right)e^{2\phi}+\frac{h^{2}}{6a^{6}}$
$\displaystyle=-H^{2}-\dot{H}+\dot{\phi}H-\frac{2}{3}\dot{\phi}^{2}+\ddot{\phi}$
(1b) $\displaystyle\frac{8\pi G}{3}\left(\rho
e^{2\phi}-\frac{1}{2}T_{(0)}\right)$
$\displaystyle=-H^{2}-\dot{H}+\frac{2}{3}\ddot{\phi}.$ (1c)

To have a complete set of equations, we can write the continuity equation as
below;

$\displaystyle\dot{\rho}+3H\left(\rho+p\right)+\dot{\phi}T_{(0)}e^{-2\phi}=0.$
(2)

The analytic solution of the continuity equation (2) can be found by
considering non-interacting perfect fluids with constant equation of state
parameters, $\omega$, and the parameter $\lambda$ describing the coupling of
dilaton to matter:

$\displaystyle\rho=\rho_{0}a^{-3\left(1+\omega\right)}e^{-\lambda\phi}$ (3)

where,

$\displaystyle\omega\equiv p/\rho,\qquad\lambda\equiv T_{(0)}\left(\rho
e^{2\phi}\right)^{-1}$ (4)

An autonomous equation system can be constructed starting with defining the
boundaries of the phase space by dividing both sides of the eq. (1a) by
$H^{2}$ as,

$\displaystyle\frac{8\pi G}{3H^{2}}\rho
e^{2\phi}+\frac{h^{2}}{12a^{6}H^{2}}+2\frac{\dot{\phi}}{H}-\frac{2}{3}\left(\frac{\dot{\phi}}{H}\right)^{2}=1$
(5)

for which the eq. (5) now gives the constraint of our phase space. We use the
constraint equation as our reference to define dimensionless phase space
dynamical variables, that ultimately will describe the behavior of the kinetic
term of the scalar field, the flux term, and the energy density term in our
equations.

$\displaystyle x\equiv\frac{\dot{\phi}}{H}\quad\quad
y\equiv\frac{h}{2\sqrt{3}a^{3}H}\quad\quad\Omega\equiv\frac{\sqrt{8\pi
G}}{\sqrt{3}H}\sqrt{\rho}e^{\phi}$ (6)

Substituting the dynamical variables into the constraint equation (5) yields
the following form that defines the phase space of the hyperboloid of two
sheets,

$\displaystyle\Omega^{2}+y^{2}+2x-\frac{2}{3}x^{2}=1$ (7)

The autonomous equation system can be constructed by calculating the evolution
of the dynamical variables in the phase space by taking their derivative with
respect to $dN\equiv Hdt$, denoted by the prime (′) symbol.

$\displaystyle x^{\prime}=\frac{\ddot{\phi}}{H^{2}}-x\frac{\dot{H}}{H^{2}}$
(8a) $\displaystyle y^{\prime}=-y\left(3+\frac{\dot{H}}{H^{2}}\right)$ (8b)
$\displaystyle\Omega^{\prime}=\Omega\left(\frac{1}{2}\frac{\dot{\rho}}{H\rho}+x-\frac{\dot{H}}{H^{2}}\right)$
(8c)

Explicit forms of the terms in the above equations can be written in terms of
the dynamical variables as follows.

$\displaystyle\frac{\dot{H}}{H^{2}}=\Omega^{2}\left(3\omega+\frac{3\lambda}{2}-2\right)-2x+\frac{4x^{2}}{3}+4y^{2}-1$
(9a)
$\displaystyle\frac{\ddot{\phi}}{H^{2}}=\frac{3}{2}\left(\frac{\dot{H}}{H^{2}}+1\right)+\frac{3\Omega^{2}}{2}\left(1-\frac{\lambda}{2}\right)$
(9b) $\displaystyle\frac{\dot{\rho}}{H\rho}=-3(1+\omega)-x\lambda$ (9c)

The region of acceleration can be realized where $\ddot{a}$ has positive
values. In terms of the dynamical variables and equation of state parameters,
this can be written as,

$\displaystyle\dfrac{\dot{H}}{H^{2}}+1=\dfrac{\ddot{a}}{a}>0$ (10)
$\displaystyle\Omega^{2}\left(3\omega+\frac{3\lambda}{2}-2\right)-2x+\frac{4x^{2}}{3}+4y^{2}>0$
(11)

As it has been stated in Angus2020 , since we take the equation of state and
the dilaton’s coupling parameter as constants, the flux can be taken as zero
for which its contribution can then be incorporated into the energy density by
setting specific values of $\omega$ and $\lambda$, i.e.
$\left(\omega,\lambda\right)_{\text{h}}=(1,2)$. Therefore, from now on, we
will set the flux contribution to zero, $h=0$, and we will investigate the
contribution of the flux term by setting parameters accordingly. Moreover, by
setting $h=0$ we reduced our phase space dimension to 2 with phase space only
constructed with $x$ and $\Omega$.

The initial phase space of the system has the hyperboloid of two sheets form
as can be seen in 1(a). One can make a variable redefinition as shown below to
reconstruct the phase space as a unit sphere,

$\displaystyle\tilde{x}\equiv\frac{\sqrt{3}}{2\left(x-3/2\right)}\quad\quad\tilde{y}\equiv\frac{\sqrt{3}y}{\sqrt{2}\left(x-3/2\right)}\quad\quad\tilde{\Omega}\equiv\frac{\sqrt{3}\Omega}{\sqrt{2}\left(x-3/2\right)}$
(12)

for which the phase space has now taken the form as in 1(b) and the constraint
equation now reads,

$\displaystyle\tilde{x}^{2}+\tilde{y}^{2}+\tilde{\Omega}^{2}=1.$ (13)

(a) Hyperboloid phase space in eq. (7)

(b) Spherical phase space in eq. (13)

Figure 1: Visual representation of the phase spaces.

The evolution of the phase space variables can be found as before while taking
into account the evolution of the old set in terms of the new variable set. It
takes a similar form as before with the addition of new terms as,

$\displaystyle\tilde{x}^{\prime}=-\frac{2\tilde{x}^{2}}{\sqrt{3}}\left(\frac{\ddot{\phi}}{H^{2}}-\left(\frac{3}{2}+\frac{\sqrt{3}}{2\tilde{x}}\right)\frac{\dot{H}}{H^{2}}\right),$
(14a)
$\displaystyle\tilde{y}^{\prime}=-\tilde{y}\left(3+\frac{\dot{H}}{H^{2}}\right)-\frac{2\tilde{x}\tilde{y}}{\sqrt{3}}\left(\frac{\ddot{\phi}}{H^{2}}-\left(\frac{3}{2}+\frac{\sqrt{3}}{2\tilde{x}}\right)\frac{\dot{H}}{H^{2}}\right),$
(14b)
$\displaystyle\tilde{\Omega}^{\prime}=\tilde{\Omega}\left(\frac{1}{2}\frac{\dot{\rho}}{H\rho}+\frac{3}{2}+\frac{\sqrt{3}}{2\tilde{x}}-\frac{\dot{H}}{H^{2}}-\frac{2\sqrt{3}}{3}\tilde{x}\left(\frac{\ddot{\phi}}{H^{2}}-\left(\frac{3}{2}+\frac{\sqrt{3}}{2\tilde{x}}\right)\frac{\dot{H}}{H^{2}}\right)\right).$
(14c)

Again we kept some of the terms closed to improve the readability of the
autonomous equation set and provide their explicit form below.

$\displaystyle\frac{\dot{H}}{H^{2}}=\frac{1}{2}\frac{\tilde{\Omega}^{2}}{\tilde{x}^{2}}\left(3\omega+\frac{3\lambda}{2}-2\right)+2\frac{\tilde{y}^{2}}{\tilde{x}^{2}}+\frac{1}{\tilde{x}^{2}}+\frac{\sqrt{3}}{\tilde{x}}-1$
(15a)
$\displaystyle\frac{\ddot{\phi}}{H^{2}}=\frac{3}{2}\left(\frac{\dot{H}}{H^{2}}+1\right)+\frac{3\tilde{\Omega}^{2}}{4\tilde{x}^{2}}\left(1-\frac{\lambda}{2}\right)$
(15b)
$\displaystyle\frac{\dot{\rho}}{H\rho}=-3(1+\omega)-\left(\frac{3}{2}+\frac{\sqrt{3}}{2\tilde{x}}\right)\lambda$
(15c)

The acceleration region of the phase space can be written as before with
$\frac{\dot{H}}{H}+1>0$ condition,

$\displaystyle\frac{1}{2}\frac{\tilde{\Omega}^{2}}{\tilde{x}^{2}}\left(3\omega+\frac{3\lambda}{2}-2\right)+2\frac{\tilde{y}^{2}}{\tilde{x}^{2}}+\frac{1}{\tilde{x}^{2}}+\frac{\sqrt{3}}{\tilde{x}}>0$
(16)

Ultimately, both phase spaces have the same physics, for which one critical
point found in one phase space can easily be converted to the critical point
found in the other phase space by utilizing the definitions of the dynamical
variables. Moreover, the discontinuity in the hyperboloid of two-sheet space
is also present in the spherical phase space. Even though the spherical phase
space looks closed, the conversion between them has a discontinuity, thus the
original behavior of the system remains unchanged. Therefore, in the rest of
the paper, we will use the original equation set in eq. (8) to determine the
critical points of the phase space and their characteristics.

##  3 Critical Points and Stability Analysis

Autonomous equation system in eq. (8) has a total of $7$ critical points
listed in Table 1 along with the critical point’s existence requirement, the
condition for the characteristics of the critical point and the point’s
behavior. We also give the corresponding bifurcation diagrams in Figure 2 to
visualize the characteristics of the critical points of the autonomous
equation system based on the values of $\omega$ and $\lambda$.

# | $x$ | $\Omega$ | Eigenvalues | Existence | Condition | Behavior
---|---|---|---|---|---|---
A | $0$ | $0$ | $\small{\begin{array}[]{ll}&\lambda_{1}=-2\\\ &\lambda_{2}=\frac{-(3\omega+1)}{2}\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\forall(\omega,\lambda)\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\omega<-1/3\\\ &\omega>-1/3\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\text{saddle}\\\ &\text{attractor}\end{array}}$
B | $\frac{3+\sqrt{3}}{2}$ | $0$ | $\small{\begin{array}[]{ll}&\lambda_{1}=-2-2\sqrt{3}\\\ &\lambda_{2}=\frac{-\lambda(3+\sqrt{3})-2(3w+\sqrt{3})}{4}\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\forall\left(\omega,\,\lambda\right)\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\lambda_{2}>0\\\ &\lambda_{2}<0\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\text{saddle}\\\ &\text{attractor}\end{array}}$
C | $\frac{3-\sqrt{3}}{2}$ | $0$ | $\small{\begin{array}[]{ll}&\lambda_{1}=-2+2\sqrt{3}\\\ &\lambda_{2}=\frac{-\lambda(3-\sqrt{3})-2(3w-\sqrt{3})}{4}\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\forall\left(\omega,\,\lambda\right)\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\lambda_{2}>0\\\ &\lambda_{2}<0\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\text{source}\\\ &\text{saddle}\end{array}}$
D± | $\frac{g(\omega,\lambda)}{\lambda+2\omega}$ | $\frac{\pm\sqrt{h(\omega,\lambda)}}{\lambda+2\omega}$ | $\small{\begin{array}[]{ll}&\lambda_{1}=\frac{g(\omega,\lambda)^{2}-3(\omega-1)^{2}}{2(\lambda+2\omega)}\\\\[5.0pt] &\lambda_{2}=2\lambda_{1}-\frac{4(\omega-1)}{\lambda+2\omega}\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\lambda\neq-2\omega\\\ &h(\omega,\lambda)\geq 0\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\lambda_{1,2}>0\\\ &\lambda_{1}\cdot\lambda_{2}<0\\\ &\lambda_{1,2}<0\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\text{source}\\\ &\text{saddle}\\\ &\text{attractor}\end{array}}$
E± | $\frac{-(3\omega+1)}{\lambda-2}$ | $\pm\sqrt{\frac{-4(3\omega+1)}{3(\lambda-2)^{2}}}$ | $\small{\begin{array}[]{ll}&\lambda_{1}=\frac{g(\omega,\lambda)+f(\omega,\lambda)}{2-\lambda}\\\ &\lambda_{2}=\frac{g(\omega,\lambda)-f(\omega,\lambda)}{2-\lambda}\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\omega\leq-1/3\\\ &\lambda\neq 2\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\lambda_{1}\cdot\lambda_{2}<0\\\ &\lambda_{1,2}<0\\\ &g<0,f^{2}<0\end{array}}$ | $\footnotesize{\begin{array}[]{ll}&\text{saddle}\\\ &\text{attractor}\\\ &\text{stable spiral}\end{array}}$
$\small{\begin{array}[]{ll}&f(\omega,\lambda)\equiv\sqrt{6\lambda^{2}\omega+3\lambda^{2}+36\lambda\omega^{2}+6\lambda\omega-6\lambda+36\omega^{3}-3\omega^{2}-2\omega+5}\\\
&g(\omega,\lambda)\equiv\lambda+3\omega-1\qquad
h(\omega,\lambda)=(\omega-1)^{2}-\frac{1}{3}g(\omega,\lambda)^{2}\end{array}}$
Table 1: Critical points of autonomous system in eq. (8)

(a) Point A

(b) Point B

(c) Point C

(d) Point D±

(e) Point E±

Figure 2: Critical point behavior analysis of Table 1 with respect to equation
of state parameters $\omega$ and $\lambda$. White curves between regions
represent bifurcation curves, changes in the thickness of the lines are purely
for the enhancement of the visual representation. Green: Source node, Blue:
Saddle node, Red: Attractor node, Grey: Stable spiral, Black line
($\lambda=-2\omega$, $\lambda=-2$) and the shaded regions
($h(\omega,\lambda)<0$, $\omega>-1/3$) represents the forbidden regions.

  * •

Point A always exists without any restrictions coming from the parameters
$\omega$ and $\lambda$. It lies at the origin of the phase space as
$(x,\Omega)=(0,0)$ with the behavior depending solely on the parameter
$\omega$. In the case of $\omega<-1/3$ the point behaves as a saddle point,
whereas when the $\omega$ passes the value $-1/3$ the point behaves as an
attractor node. The behavior at the intersection of those two regions where
the $\omega$ takes the value exactly $-1/3$, which characterizes the curvature
contribution, can not be determined by the linear stability analysis since the
point then becomes non-hyperbolic as it will have a zero eigenvalue. As can
also be seen in the 2(a), for the majority of the allowed values of
$(\omega,\lambda)$ pair, the critical point appears as an attractor node, in
which the point corresponds to the complete flux contribution universe, i.e.
$y=1$.

  * •

Point B is independent of the other critical points, meaning that the field
lines that live in the sheet where critical point B exists, never cross to the
other sheet that houses the remaining critical points because of the
hyperboloid of two sheet characteristics which can be seen in 1(a). Moreover,
the point B always exists without any dependence on the parameters $\omega$ or
$\lambda$. One needs to choose phase space sheets from the existing two to
continue the analysis, for which we considered the sheet that constrains the
$x$ as $x<3/2$. Justification of this particular choice can be realized when
we think about the constraint equation in Equation 7. The $x$ contribution to
the constraint equation is $2x-2x^{2}/3$ and each squared contribution that is
in the constraint equation is constrained by $0$ from the below and by $1$
from the above, i.e. for
$\epsilon\in\\{2x-2x^{2}/3,\,\Omega^{2},\,y^{2}\\},0\leq\epsilon\leq 1$
demanded by the constraint. Therefore, the value of $x\geq 3/2$ is not in the
scope. Thus, we will not make any further analysis on critical point B, since
it is not in the scope of the phase space we have chosen.

  * •

Point C also exists without any dependence on the values of $\omega$ and
$\lambda$. It can have the characteristics of both a source node and a saddle
node for specific values of $\omega$ and $\lambda$. For the majority of the
values of $\omega$ and $\lambda$, the point acts as a source node. The
dominant contribution to the evolution at the critical point C comes from the
$x=\frac{\dot{\phi}}{H}$ with $\Omega=0$. The only possible source
contribution in our phase space can be realized in the critical point C and,
for the case of point A being the attractor node and point C being the source
node the restriction to the equation of state parameters can be written as
follows,

$\displaystyle-\dfrac{1}{3}<\omega<\dfrac{1}{\sqrt{3}}-\dfrac{\lambda}{2}\left(1-\dfrac{1}{\sqrt{3}}\right)$
(17)

For which, the upper bound on $\omega$ given above can be visually understood
as in the 2(c). Moreover, since there is no restriction on the range of
allowed values of the parameter $\lambda$ from $-2$ to $2$, the above
inequality has the upper bound stated by the parameter $\lambda$ as,

$\displaystyle-\dfrac{1}{3}<\omega<\begin{cases}\quad 1\ \
\,\quad\rightarrow\lambda=-2\\\
\frac{2\sqrt{3}}{3}-1\rightarrow\lambda=2\end{cases}$ (18)

In the case of $\lambda=2$ the upper bound becomes $\omega\lesssim 0.15$, for
which the choice of $\omega=1$ for the flux contribution integrated into the
$\rho$, $(\omega,\lambda)_{\rm{h}}=(1,2)$, results in with a phase space that
does not have a source node.

  * •

Point D± points can not exist simultaneously, because of the $\Omega>0$
condition. In the positive case of the critical point, the condition of
existence is constrained by the $\lambda+2\omega>0$ statement, and in the
negative case, the constraining condition becomes the opposite as
$\lambda+2\omega<0$. Without taking the conditions for existence the point’s
stability characteristics can be categorized as; for the negative case, the
condition reduces the possible values that can be taken for the
$\omega,\,\lambda$ to the below of the black curve in 2(d). For this, the only
possible characteristic of the critical point becomes a saddle-node or an
attractor node. For the positive case, the point can take all three
characteristics as a saddle, an attractor, or a source node without the change
in location of the point in phase space. Moreover, the shaded areas in the
figure represent the values of $\omega$ and $\lambda$ where the
$h(\omega,\lambda)$ becomes negative, which renders $\Omega$ parameter
imaginary, therefore can not be taken as a possible solution and completely
removes the source characteristics from the critical point.

  * •

Point E±’s condition of existence is restricted by both parameters $\omega$
and $\lambda$ where the restricted region that is dependent on the parameter
$\omega$ is given as $\omega>-1/3$, for which makes the $\Omega$ variable
complex-valued as can be seen in 2(e). On the other hand, the restriction
coming from the parameter $\lambda$ is for its specific value $\lambda=2$
which in that case both dynamical variables $x$ and $\Omega$ go to infinity.
In the allowed region where $\omega\geq-1/3$ and $\lambda\neq 2$, the critical
point E can be a saddle-node, an attractor node, or a stable spiral node based
on the values of $\omega$ and $\lambda$ that is given in the Table 1.

##  4 Cosmologically Interesting Cases

We investigate the phase space dynamics with different parameter sets of
cosmologically interesting cases. The initial condition values for the
dynamical variables $x$ and $\Omega$ are chosen such that the evolution of the
effective equation of state parameter $\omega_{\rm{eff}}$ starts with the
value close to that of the radiation, i.e. $\omega_{\rm{eff}}=1/3$, where the
$\omega_{\rm{eff}}$ is defined as,

$\displaystyle\omega_{\rm{eff}}\equiv-1-\dfrac{2}{3}\dfrac{\dot{H}}{H^{2}}.$
(19)

We also provide in the figures below, the evolution of the deceleration
parameter $q$ defined as,

$\displaystyle q\equiv-1-\dfrac{\dot{H}}{H^{2}}$ (20)

The deceleration parameter $q$ is useful to determine the acceleration of the
universe: the universe is accelerating when $q<0$ and decelerating when $q>0$.
We reintroduced the dynamical variable $y$ as a derived parameter from the
constraint equation eq. (7). We had previously set it to zero without loss of
generality as the flux contribution can be integrated into the energy density
by employing specific values for the equation of state parameters and provided
its evolution along with other two dynamical variables $x$ and $\Omega$.
Consequently, the flux domination, $y=1$ is the attractor of the system for
the values of $\omega>-1/3$ and no restrictions on the value of $\lambda$, as
it can be seen in the bifurcation figures Figure 2. Apart from the
cosmological constant which has the $\omega=-1$, the critical point $A$ is the
attractor of the phase space. Below, we present the phase space trajectories
of configurations, accompanied by the evolution of the dynamical variables
$x$, $y$ and $\Omega$, the effective equation of state $\omega_{\rm{eff}}$ and
the deceleration parameter $q$ for a given $w$, $\lambda$ and the initial
condition for the $x$ and $\Omega$.

###  4.1 Scalar field, $(\lambda=-2\omega),\,(\omega,\lambda)=(-1/3,2/3)$

In the case of scalar-field-like contribution to the energy density $\rho$ we
take the equation of state parameters in such a way that $\lambda=-2\omega$.
By choosing the value of the $\omega$ as $-1/3$ we arrive at the following
phase space diagram and the parameter evolution in Figure 3. There are only
two critical points in the phase space $A$ and $C$ corresponding to the points
with the same name in Table 1. The acceleration region of the phase space can
be seen with the blue shaded area on the left-hand side of the Figure 3. The
phase space trajectories start from the source node $C$ for which the
dominating contribution comes from the kinetic energy of the scalar field and
evolves to the attractor node $A$ for which both $x$ and $\Omega$ at that
point vanishes with the only contribution coming from the flux term $y$. The
evolution of the parameters of the model can be seen on the right-hand side of
the above figure. One can see that the universe begins its acceleration phase
when the flux contribution starts to dominate the universe. Furthermore, as
the flux domination settles to unity, the effective equation of state reaches
the value of $\omega_{\rm{eff}}=-1/3$.

Figure 3: On the left-hand side, we give the phase space trajectory for the
scalar field with initial conditions $(x,\Omega)=(0.62,0.05)$ that corresponds
to the source node C. The blue shaded area corresponds to the accelerating
region of the phase space with $\frac{\dot{H}}{H^{2}}+1>1$. On the right-hand
side, with the same initial conditions, we give the evolution of the dynamical
system variables ($\Omega$, $x$, and $y$), effective equation of state
parameter ($\omega_{\rm{eff}}$) and the deceleration parameter ($q$).
Horizontal dashed lines correspond to the values $\pm 1/3$.

###  4.2 Radiation, $(\omega,\lambda)=(1/3,0)$

In the configuration for a radiation contribution, we have a total of three
critical points in the phase space as can be seen on the left-hand-side of the
Figure 4. Taking the initial conditions for the variables $x$ and $\Omega$ as
$(0.28,0.69)$, the field line that shows the evolution of the phase space in
this configuration starts from the source point $C$ and follows a path towards
the saddle-node $D$ and eventually the trajectory resolves to the attractor-
node $A$. The dominant contributions for each node in the phase space can be
understood by looking at the right-hand side of the same figure, for which the
universe now starts in the radiation-dominated phase with the major
contribution coming from the variable $\Omega$. In this starting
configuration, the initial effective equation of state has the value of
$\omega_{\rm{eff}}\simeq 1/3$, while the deceleration parameter starts close
to unity, $q\simeq 1$. As the trajectory approaches the saddle node, the flux
contribution starts to dominate and eventually reaches unity, i.e. $y=1$ when
the trajectory arrives at the attractor-node $A$. The effective equation of
state settles on the value $w_{\rm{eff}}=-1/3$ and the deceleration parameter
reaches $0$ but does not go beyond that point for which the universe
accelerates, $q\nless 0$. In other words, the configuration does not evolve to
a universe where it accelerates. During the evolution, the kinetic energy term
of the scalar field relaxes to the value of $0$.

Figure 4: The phase space trajectory and parameter evolution for the radiation
contribution with the initial conditions $(x,\Omega)=(0.28,0.69)$ where again
the blue shaded area corresponds to the accelerating region of the phase space
with $\frac{\dot{H}}{H^{2}}+1>1$. Horizontal dashed lines correspond to the
values $\pm 1/3$.

###  4.3 Dust, $(\omega,\lambda)=(0,0)$

In the case of a dust contribution, the equation of state parameters takes the
value $(\omega,\lambda)=(0,0)$ for which it also coincides with the scalar
field contribution since in general it is given as $\lambda=-2\omega$. In
Figure 5, we present the phase space diagram for the dust configuration and
parameter evolution for the initial condition $(x,\Omega)=(0.55,0.30)$.
Similar to the scalar field case, phase space has two critical points, $A$ and
$C$ being attractor and source nodes respectively. The initial value for the
effective equation of state, $\omega_{\rm{eff}}$, has a value close to $1/3$,
and after evolution, it settles on the value $-1/3$. The evolution graph and
the trajectory show that the acceleration again starts as the flux
contribution starts to dominate the universe. The kinetic energy of the scalar
field contribution, $x$, rapidly decreases as the universe evolves until the
flux domination, after which it starts to grow and eventually resolves on the
zero.

Figure 5: The phase space trajectory and parameter evolution for dust
contribution with the initial values are taken as $(x,\Omega)=(0.55,0.30)$.
The black dashed line corresponds to the effective equation of state
evolution, whereas the red dashed line corresponds to the deceleration
parameter evolution.

##  5 Chameleon Cosmology

The chameleon mechanism is a scenario that, in accordance with string theory
predictions, permits scalar fields to evolve on current cosmological time
scales with couplings of order unity to matter. This mechanism is considered
as one of the so-called ‘screening’ mechanisms222A closely related one is the
Damour-Polyakov mechanism and the other is Vainstein mechanism. in which the
mass of the scalar field depends on the density of local matter: the field can
be sufficiently massive to fulfill limits on equivalence principle violations
and fifth force in high-density areas like the solar system but in
cosmological scales evolve with a mass of the order of the Hubble parameter.
The screened dynamics of the scalar field in the chameleon mechanism is the
result of the interplay between a non-linear potential and a linear coupling
to matter. In a symmetry-dictated, and thus more physically motivated setting
like $O(D,D)$-symmetric formulation, the coupling of the dilaton to matter
sector is determined on firm ground and anticipated to include the chameleon-
like effects. In this section, we will elaborate on this issue.

We have the action for the Chameleon cosmology as shown below with the
function $f\left(\phi\right)$ being the coupling to the matter sector,

$\displaystyle S=\int
d^{4}x\sqrt{-g}\left\\{\dfrac{R}{2}-\dfrac{1}{2}g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi-V\left(\phi\right)+f(\phi)L_{\text{m}}\right\\}\quad
L_{\text{m}}\sim\rho_{\text{m}}$ (21)

Varying the action in eq. (21) with respect to the scalar field $\phi$ and
inverse metric we arrive at the following equations of motion,

$\displaystyle\ddot{\phi}+3H\dot{\phi}+V^{\prime}(\phi)-f^{\prime}(\phi)\rho_{\text{m}}=0$
(22) $\displaystyle 3H^{2}=\dfrac{\dot{\phi}^{2}}{2}+V+f(\phi)\rho_{\text{m}}$
(23)

Furthermore, by taking the derivative of the eq. (23) and using the equation
of motion for the scalar field in eq. (22), we have the following expression,

$\displaystyle
2\dot{H}+3H^{2}=-\left(\dfrac{\dot{\phi}^{2}}{2}-V(\phi)+f(\phi)p_{\text{m}}\right)$
(24)

The final expression we will utilize in the linear stability analysis of the
Chameleon cosmology is the conservation equation which can be written as
follows,

$\displaystyle\dot{\rho}_{\text{m}}+3H(\rho_{\text{m}}+p_{\text{m}})+2\dfrac{f^{\prime}}{f}\dot{\phi}\rho_{\text{m}}=0$
(25)

As previously done, we will construct the constraint equation for the phase
space by dividing both sides of the eq. (23) by $3H^{2}$ and define the phase
space dynamical variables accordingly.

$\displaystyle
1=\dfrac{\dot{\phi}^{2}}{6H^{2}}+\dfrac{V}{3H^{2}}+\dfrac{f\rho_{\text{m}}}{3H^{2}}$
(26)

Thus, we have the following dynamical variables for the phase space of
Chameleon cosmology,

$\displaystyle x\equiv\dfrac{\dot{\phi}}{\sqrt{6}H}\qquad
y\equiv\dfrac{\sqrt{V}}{\sqrt{3}H}\qquad\Omega\equiv\dfrac{f\rho_{\text{m}}}{3H^{2}},\quad
f(\phi)\geq 0,\ \Omega\geq 0$ (27)

Therefore, the constraint equation in eq. (26) can be rewritten using the
dynamical variables of the phase space of the Chameleon setting as,

$\displaystyle 1=x^{2}+y^{2}+\Omega$ (28)

One can immediately see that the dynamical variable $x$ again gives the scalar
field’s kinetic energy contribution to the evolution with an overall constant
factor. Moreover, the variable $\Omega$ in the Chameleon setting in eq. (27)
resembles the one in $O(D,D)$ configuration defined with the same name in eq.
(6). Lastly, the variable $y$ in the Chameleon setting provides the
contribution of the scalar field’s potential, yet in the $O(D,D)$ setting it
shows the flux contribution, $h$.

The autonomous equation system of the Chameleon cosmology can now be found by
again taking the time derivative of the dynamical variables and dividing it by
$H$ which we denote it as a prime, ′, followed by rewriting them in terms of
the dynamical variables of the phase space,

$\displaystyle
x^{\prime}=\dfrac{1}{\sqrt{6}}\dfrac{\ddot{\phi}}{H^{2}}-x\dfrac{\dot{H}}{H^{2}}$
(29) $\displaystyle
y^{\prime}=y\left(\sqrt{\dfrac{3}{2}}\dfrac{V^{\prime}}{V}x-\dfrac{\dot{H}}{H^{2}}\right)$
(30)
$\displaystyle\Omega^{\prime}=-\Omega\left(\sqrt{6}\dfrac{f^{\prime}}{f}x+3(1+w)+2\dfrac{\dot{H}}{H^{2}}\right).$
(31)

Explicit forms of the acceleration term of the scalar field and the derivative
of the Hubble constant $H$ are given as follows,

$\displaystyle\dfrac{\dot{H}}{H^{2}}=-3\left(x^{2}+\dfrac{1}{2}(1+w)\Omega\right)$
(32)
$\displaystyle\dfrac{\ddot{\phi}}{H^{2}}=-3\left(\sqrt{6}x+\dfrac{V^{\prime}}{V}y^{2}-\dfrac{f^{\prime}}{f}\Omega\right)$
(33)

We now choose the potential $V(\phi)$ and the coupling $f(\phi)$ in
exponential form as follows:

$\displaystyle V(\phi)\propto e^{A\phi},\qquad f(\phi)\propto e^{B\phi},$ (34)

where $A$ and $B$ are constants. With this choice, the systems of equations
for the $O(D,D)$ extended system and the Chameleon cases are manifestly
comparable. One can see that the coupling term $f\left(\phi\right)\propto
e^{B\phi}$ is closely related to the parameter $\lambda$ in the $O(D,D)$
complete formulation since both measure the scalar field’s coupling to the
matter sector. Furthermore, the scalar potential term, $y$, in the dynamical
system formulation of the Chameleon cosmology eq. (27) corresponds to the flux
term, $\propto h$ in eq. (6). Thus, to compare the Chameleon cosmology and the
$O(D,D)$ complete cosmology in the setting of linear stability analysis one
can set $A$ and $B$ constants and the equation of state parameter $\omega$
accordingly, so that the contribution of each term in Chameleon cosmological
setting corresponds to the ones in the previously examined cosmological cases
in the $O(D,D)$ cosmology.

In our $O(D,D)$-complete cosmology linear stability analysis, we set the flux
term to zero without loss of generality. Therefore, to match the contribution
in the Chameleon formulation we set $A$ to be zero to have constant scalar
field potential. Furthermore, by direct comparison of the $\Omega$ parameters
in eq. (6) and eq. (27) we set $B=2$.

We performed the same linear stability analysis to the corresponding Chameleon
cosmology counterparts for all three previously studied cases of $O(D,D)$
cosmologies. The dynamical system phase space in eq. (29) has three distinct
critical points with source characteristics at the node
$\left(x,\,\Omega\right)=(-1,\,0)$, with saddle characteristics at
$\left(x,\,\Omega\right)=(1,\,0)$, and with the attractor node at the
$\left(x,\,\Omega\right)=(0,\,0)$ which is named $A$, $B$, and $C$
respectively.

###  5.1 Scalar field, $\omega=-1/3$

To study the scalar field in the context of Chameleon cosmology, we take the
$\omega$ as $-1/3$, and on the left-hand side of Figure 6 we give the
resulting phase space trajectories along with on the right-hand side the phase
space parameter evolution can be seen.

Figure 6: Phase space evolution for scalar field contribution in the Chameleon
cosmology. Initial values are taken as $(x,\Omega)=(-0.97,0.05)$.

Phase space trajectories starts from the source node $A(x,\,\Omega)=(-1,0)$
goes near the saddle point $B(x,\,\Omega)=(1,0)$ and finally reaches the
attractor node $C(x,\,\Omega)=(0,0)$. The accelerating region again is a blue-
shaded area where the brief acceleration occurs with $\Omega$ domination after
which the second acceleration happens when the scalar field potential term $y$
comes to domination and eventually reaches the attractor node $C$. Comparison
with the $O(D,D)$ symmetric setting shows that the effective equation of state
$w_{\textrm{eff}}$ in the Chameleon configuration relaxes to $-1$, whereas in
the $O(D,D)$ analysis the effective equation of state goes to $-1/3$ value. By
looking at the deceleration parameter $q$, one can see that in the scalar
field case of the Chameleon setting, there are two accelerating regions for
which the first acceleration happens during the $\Omega$ domination and the
second occurs when the scalar field potential term $y$ comes to dominate in
the neighbor of the attractor node $C$ at the end of evolution which is a
similar case in the $O(D,D)$ setting where the attractor point is the one
where the flux term dominates.

###  5.2 Radiation, $\omega=1/3$

In Figure 7 we give the phase space diagram and parameter evolution in the
radiation setting. The phase space has a source node, a saddle node, and an
attractor node for which we give the initial condition for the start of the
trajectory to coincide with the source node as $(x,\,\Omega)=(-0.97,0.05)$.
With the given initial condition, the acceleration begins only when the scalar
field potential term, $y$, becomes the dominating factor as the deceleration
parameter goes below $0$ and eventually reaches $-1$, where the trajectory
reaches the attractor node at $(x,\,\Omega)=(0,0)$. The attractor region
window significantly reduces in the radiation case for which the two
accelerating regions in the previous case are now not apparent. Acceleration
only happens when the potential of the scalar field terms becomes the
dominating factor of the system. Similar to the scalar field case in the
Chameleon setting, in the radiation case, the effective equation of state goes
to $-1$ as opposed to $-1/3$ which was the case in the $O(D,D)$ complete
cosmology analysis. The attractor point of the Chameleon system again matches
to the one of a flux contribution in the $O(D,D)$ setting.

Figure 7: Phase space evolution for radiation contribution in the Chameleon
cosmology. On the left-hand side, phase space trajectory with the initial
condition $(x,\Omega)=(-0.97,0.05)$ is plotted with the red line The blue
shaded region represents the accelerated region and the black solid line shows
the constraint. On the right-hand side, we give the parameter evolution as
same as in the $O(D,D)$ symmetric setting.

###  5.3 Dust, $\omega=0$

By taking the equation of state parameter as zero, we can study the
contribution of dust, i.e. pressureless matter. In Figure 8, we give the phase
space portrait of the dynamical system along with the parameter evolution. The
acceleration window is narrower than the scalar field contribution yet it is
wider than the radiation one. The starting point of the trajectory given in
the left-hand side of the figure is taken as before,
$(x,\,\Omega)=(-0.97,0.05)$ such that it originates from the source point. The
acceleration begins again when the potential term dominates which is the case
for all analyzed systems for both $O(D,D)$ complete cosmology setting and the
Chameleon cosmology setting. The equation of state parameter
$w_{\textrm{eff}}$ again goes to the one of a scalar field or a cosmological
constant $-1$, while the potential term slowly takes over the domination.

Figure 8: Phase space evolution for dust contribution in the Chameleon
cosmology. Initial values are taken as $(x,\Omega)=(-0.97,0.05)$.

In all three cases that we have analyzed, similar aspects of $O(D,D)$ complete
cosmology and the Chameleon cosmology show themselves in the eventual
dominating factors of the system, which are the flux contribution and the
potential of the scalar field, respectively. The striking difference is that
the final value for the effective equation of state value reaches to the one
of a scalar field or a cosmological constant $-1$ in the Chameleon setting,
whereas in the $O(D,D)$ complete cosmology in all cases the effective equation
of state goes to $-1/3$. The area of acceleration regions in the Chameleon
cosmology significantly depends on the contribution of the contents of the
system whereas in the $O(D,D)$ setting the difference in the accelerating
regions is not that striking almost constant yet in the $O(D,D)$ symmetric
setting, the acceleration window of the system very much depends on the
content for which it can be seen in the case of radiation contribution of the
$O(D,D)$ cosmology where the acceleration barely happens.

##  6 Results and Discussion

One of the most elegant approaches to studying the cosmological equations that
describe the dynamics of a homogeneous and isotropic Universe is to transform
them into dynamical systems. The equations of such a maximally symmetric
cosmological setting form systems of ordinary differential equations that can
be brought into the form of autonomous systems. This enables a quantitative
comprehension of the cosmic dynamics provided by the models under study to be
obtained through the application of strong analytical and numerical
approaches. In this manuscript, we incorporate these methods in the context of
$O(D,D)$-complete Friedmann cosmology. In doing this, one of our primary aims
is to study the dynamical system approach for comparing the $O(D,D)$-extended
Friedmann cosmology with the Chameleon cosmology to see if the
$O(D,D)$-symmetric system has possible phenomenological applications to both
early and late universe problems. It is expected at the onset of analysis that
the $O(D,D)$-extended cosmological system has richer possibilities, and this
is realized in the form of a hyperboloid of two sheets in the phase space
construction. On the other hand, the two-sheeted hyperboloid can be put into
the form of a sphere, and we prefer to consider such a compact region in phase
space for our analysis. In the cosmologically interesting cases of scalar
fields, radiation, and non-relativistic matter there is at least an attractor
point and a source point. We also considered the physical meaning and
importance of these critical points for whether they lead to acceleration (or
deceleration) of the universe. We compared with the similar analysis of the
Chameleon cosmology and saw that accelerating regions in both models coincide
with the flux and scalar field potential domination in $O(D,D)$ and Chameleon
models, respectively. Furthermore, in $O(D,D)$ setting, when the flux
contribution comes to domination and eventually reaches the attractor node,
the effective equation of state parameter relaxes to $-1/3$; but in the
Chameleon setting, when the scalar field potential starts to dominate, the
effective equation of state parameter goes to $-1$.

The reason for the comparison with the chameleon cosmic model was that we
expected that the $O(D,D)$-extended system would give rise to similar, or even
richer, dynamics when analyzed and compared in terms of dynamical systems, and
the aforementioned facts confirm this expectation. This would allow
$O(D,D)$-symmetric systems to provide a physically sound alternative to
general relativity in addressing cosmic problems, instead of modifications
suggested on pure phenomenological ground. There are some other important
cosmological problems in this setting, such as checking the generation of
scale-free cosmological perturbations, properties of a possible stochastic
gravitational wave background, and nucleosynthesis effects. Furthermore, the
coupling of the dilaton to matter part in a certain symmetry-dictated way can
lead to observable effects on the structure of compact objects like neutron
stars. These are possible further research problems that we are working on.

##  7 Acknowledgement

Authors are supported by The Scientific and Technological Research Council of
Türkiye (TÜBİTAK) through the ARDEB 1001 project with grant number 121F123.

## References

  * (1) T. Damour and A.M. Polyakov, _String theory and gravity_ , _General Relativity and Gravitation_ 26 (1994) 1171 [gr-qc/9411069].
  * (2) T. Damour and A.M. Polyakov, _The string dilation and a least coupling principle_ , _Nuclear Physics B_ 423 (1994) 532 [hep-th/9401069].
  * (3) M. Gasperini and G. Veneziano, _The pre-big bang scenario in string cosmology_ , __ 373 (2003) 1 [hep-th/0207130].
  * (4) A.A. Tseytlin and C. Vafa, _Elements of string cosmology_ , _Nuclear Physics B_ 372 (1992) 443 [hep-th/9109048].
  * (5) M. Gasperini, _Dilaton Cosmology and Phenomenology_ , in _String Theory and Fundamental Interactions_ , M. Gasperini and J. Maharana, eds., vol. 737, p. 787 (2008), DOI.
  * (6) K.A. Meissner and G. Veneziano, _Manifestly O(d,d) Invariant Approach to Space-Time Dependent String Vacua_ , _Modern Physics Letters A_ 6 (1991) 3397 [hep-th/9110004].
  * (7) M. Gasperini and G. Veneziano, _O(d, d)-covariant string cosmology_ , _Physics Letters B_ 277 (1992) 256 [hep-th/9112044].
  * (8) P. Brax, _Screening mechanisms in modified gravity_ , _Classical and Quantum Gravity_ 30 (2013) 214005.
  * (9) T. Damour and G. Esposito-Farese, _Nonperturbative strong-field effects in tensor-scalar theories of gravitation_ , __ 70 (1993) 2220.
  * (10) D.D. Doneva, F.M. Ramazanoǧlu, H.O. Silva, T.P. Sotiriou and S.S. Yazadjiev, _Spontaneous scalarization_ , _Reviews of Modern Physics_ 96 (2024) 015004 [2211.01766].
  * (11) C. Hull and B. Zwiebach, _Double field theory_ , _Journal of High Energy Physics_ 2009 (2009) 099 [0904.4664].
  * (12) O. Hohm, C. Hull and B. Zwiebach, _Generalized metric formulation of double field theory_ , _Journal of High Energy Physics_ 2010 (2010) 8 [1006.4823].
  * (13) O. Hohm, C. Hull and B. Zwiebach, _Background independent action for double field theory_ , _Journal of High Energy Physics_ 2010 (2010) 16 [1003.5027].
  * (14) S. Angus, K. Cho and J.-H. Park, _Einstein double field equations_ , _European Physical Journal C_ 78 (2018) 500 [1804.00964].
  * (15) H. Wu and H. Yang, _Double field theory inspired cosmology_ , __ 2014 (2014) 024 [1307.0159].
  * (16) C.T. Ma and C.M. Chen, _Cosmological implications from O(D, D)_ , _Fortschritte der Physik_ 62 (2014) 921 [1405.4073].
  * (17) R. Brandenberger, R. Costa, G. Franzmann and A. Weltman, _Point particle motion in double field theory and a singularity-free cosmological solution_ , __ 97 (2018) 063530 [1710.02412].
  * (18) R. Brandenberger, R. Costa, G. Franzmann and A. Weltman, _Dual spacetime and nonsingular string cosmology_ , __ 98 (2018) 063521 [1805.06321].
  * (19) R. Brandenberger, R. Costa, G. Franzmann and A. Weltman, _T -dual cosmological solutions in double field theory_ , __ 99 (2019) 023531 [1809.03482].
  * (20) H. Bernardo, R. Brandenberger and G. Franzmann, _T -dual cosmological solutions in double field theory. II._ , __ 99 (2019) 063521.
  * (21) S. Angus, K. Cho, G. Franzmann, S. Mukohyama and J.-H. Park, _O (D ,D ) completion of the Friedmann equations_ , _European Physical Journal C_ 80 (2020) 830 [1905.03620].
  * (22) S. Bahamonde, C.G. Böhmer, S. Carloni, E.J. Copeland, W. Fang and N. Tamanini, _Dynamical systems applied to cosmology: Dark energy and modified gravity_ , __ 775 (2018) 1 [1712.03107].
  * (23) H. Bernardo and G. Franzmann, _$\alpha$ ’-cosmology: solutions and stability analysis_, _Journal of High Energy Physics_ 2020 (2020) 73 [2002.09856].
  * (24) J. Khoury and A. Weltman, _Chameleon Fields: Awaiting Surprises for Tests of Gravity in Space_ , __ 93 (2004) 171104 [astro-ph/0309300].
  * (25) J. Khoury and A. Weltman, _Chameleon cosmology_ , __ 69 (2004) 044026 [astro-ph/0309411].
  * (26) P. Brax, C. van de Bruck, A.-C. Davis, J. Khoury and A. Weltman, _Detecting dark energy in orbit: The cosmological chameleon_ , __ 70 (2004) 123518 [astro-ph/0408415].