# Correspondence learning between morphologically different robots through
task demonstrations

Hakan Aktas1,2, Yukie Nagai2, Minoru Asada3,4, Erhan Oztop3,5, Emre Ugur1
1Hakan Aktas and Emre Ugur are with the Department of Computer Engineering,
Bogazici University, Istanbul, Turkey<EMAIL_ADDRESS>Aktas and
Yukie Nagai are with IRCN, the University of Tokyo, Tokyo, Japan 3Minoru Asada
and Erhan Oztop are with Symbiotic Intelligent Systems Research Center,
Institute for Open and Transdisciplinary, Research Initiatives, Osaka
University, Japan 4Minoru Asada is also affiliated with the International
Professional University of Technology in Osaka, Osaka, Japan 5Erhan Oztop is
also affiliated with the Department of Computer Engineering, Ozyegin
University, Istanbul, Turkey.

###### Abstract

We observe a large variety of robots in terms of their bodies, sensors, and
actuators. Given the commonalities in the skill sets, teaching each skill to
each different robot independently is inefficient and not scalable when the
large variety in the robotic landscape is considered. If we can learn the
correspondences between the sensorimotor spaces of different robots, we can
expect a skill that is learned in one robot can be more directly and easily
transferred to the other robots. In this paper, we propose a method to learn
correspondences between robots that have significant differences in their
morphologies: a fixed-based manipulator robot with joint control and a
differential drive mobile robot. For this, both robots are first given
demonstrations that achieve the same tasks. A common latent representation is
formed while learning the corresponding policies. After this initial learning
stage, the observation of a new task execution by one robot becomes sufficient
to generate a latent space representation pertaining to the other robot to
achieve the same task. We verified our system in a set of experiments where
the correspondence between two simulated robots is learned (1) when the robots
need to follow the same paths to achieve the same task, (2) when the robots
need to follow different trajectories to achieve the same task, and (3) when
complexities of the required sensorimotor trajectories are different for the
robots considered. We also provide a proof-of-the-concept realization of
correspondence learning between a real manipulator robot and a simulated
mobile robot.

##  I INTRODUCTION

We observe a large variety of robots in terms of their morphology, sensors,
and actuators, such as mobile robots, fixed-based manipulators, and aerial and
underwater robots. While the robots are designed for the environments they
operate, there are commonalities in the required skills, such as grasping,
carrying, and releasing objects, grasping and using tools, and pushing and
pulling objects. In general, robots are equipped with these skills by using
the Learning from Demonstration approach with additional controller tuning
when needed. Given the commonalities in the skill sets, teaching each skill to
each different robot independently is not efficient and not scalable when the
large variety in the robotic landscape is considered. If we can learn the
correspondences between the sensorimotor spaces of different robots, we can
expect a skill that is learned in one robot can be more directly and easily
transferred to the other robots.

The learning of sensorimotor correspondences between robots has been studied
in robotics, where emphasis was given to efficiency in the transfer of skills
from one robot to the others [1]. The correspondence between the robots has
been established between all pairs of states [2] either by manually forming a
shared feature space [3] or by aligning states using unsupervised manifold
alignment [4]. State alignment from local and global perspectives has also
been used with a regularized policy update [5]. Expectation-maximization-based
dynamic time warping has been utilized to align the states and subsequently
established a shared feature space by employing non-linear embedding functions
[6]. Transferring policies using modularity have also been established in
various studies [7], with some using hierarchical decoupling between agents
[8, 9]. Methods that can find unsupervised action correspondence using
unpaired[10] and unaligned [11] have also been introduced. While the previous
work learned correspondences at the single state level, we previously proposed
a method [12] that formed a common feature space that encodes the
correspondence of the trajectories rather than single states. However, [12]
used an engineering loss function in the latent space and achieved the
transfer between the same type of robots: fixed-based manipulators between 3
and 4 degrees of freedom. The use of this additional loss function requires
extra fine-tuning of the hyper-parameters.

Figure 1:  The conceptual summary of the task-based correspondence learning
problems studied in this paper. The middle, left, and right columns represent
task space, and sensorimotor (SM) spaces of the mobile and manipulator robots,
respectively. The shapes of the markers represent tasks. Each marker in the
middle plot corresponds to the required robot execution trajectory in the task
space to achieve the task. Each marker in the left and right plots represents
the required SM trajectory of the corresponding robot to achieve the task. The
dashed line between markers represents the mapping between sensorimotor and
task spaces. In (A), the robots follow the same task-space trajectory to
achieve the same task. In (B), the robots need to follow different
trajectories to achieve the same task. In (C), the robots again need to follow
different trajectories, but in this case, the complexity of the sensorimotor
mappings is also different.

In this work, we aim to learn the correspondences between robots with
different morphologies, between a fixed-based manipulator robot and a
differential-drive mobile robot, that can learn to achieve the same tasks with
significantly different bodies, actuators, and control variables. Inspired by
[13], our system forms a common latent representation between different
robots. After learning the correspondences, and given observations of a new
task execution from one robot, this common representation is used to generate
the execution trajectory of the same task for the other robot. Each skill in
each robot is encoded by a separate Conditional Neural Movement Primitives
(CNMPs) network [14], which is an encoder-decoder network that can produce
complex target trajectories from observations on these trajectories.
Observations are given to the CNMP encoder to generate a latent representation
for the target trajectory. The core idea is to blend the latent
representations of each robot’s CNMPs in order to form a common latent
representation for both robots for the same task. This common latent
representation is then fed to separate decoders of different robots to
generate the execution trajectories of the corresponding robots. The core idea
is that the coupled CNMP system is trained such that the common latent
representation, which is able to generate execution trajectories for both
robots, can be obtained even only from the observations (and therefore the
latent encoding) of a single robot by setting the blending weights
accordingly.

A series of experiments are conducted to assess the extent of correspondence
learning between the manipulator and mobile robots in tasks with different
levels of complexity. For this, we designed these experiments to address the
question of “Can task-specific correspondences between robots be learned”

  * •

when the robots need to follow the same task-space trajectory to achieve the
same task (Fig. 1 (A))?

  * •

when the robots need to follow different trajectories to achieve the same task
(Fig. 1 (B))?

  * •

when the complexities of the required motor trajectories are different for
different robots (Fig. 1 (C))?

Figure 2: The general overview of our system that learns task-based
correspondence via forming common latent representation between different
robots. Given observations from an SM trajectory of one robot, our system can
generate the full SM trajectory for the other robot to achieve the same task
by setting the blending weight p to 0 or 1. Please refer to the text for the
details of training and generation.

##  II METHOD

The general overview of our approach is provided in Fig. 2. Two (or possibly
more) robots are provided initial proxy demonstrations for the same tasks.
Depending on the embodiment and capabilities of the robots the demonstrations
might differ in terms of the used modalities and the required trajectories.
Given these demonstrations, the encoder networks (A) of the robots are trained
together, forming a common latent space that can be used to generate the
required trajectories for both robots. After training, a new task
configuration can be demonstrated for one of the robots, either by providing
an entire execution trajectory or a number of observations sampled from this
trajectory. Using the common latent representation, our system can generate
the full execution trajectory for the other robot to achieve the same task
(B). We provide the details of the skill encoding, sharing, and generation in
the rest of this section.

In our system, the skills are encoded and learned using CNMPs, and the latent
representations of the CNMPs of different robots are blended to form a common
latent representation. More formally, let R =
[$Robot_{1},Robot_{2},Robot_{3}...,Robot_{n}$] be a set of sensorimotor
information collected by n robots during action executions. The sensorimotor
information during these executions is recorded at each time step. Let I be
the multi-robot sensorimotor interaction data set, and the $i^{th}$ task
execution is shown as $I_{i}=\\{(t,S_{t}^{R})\\}_{t=0}^{T}$ where $S_{t}^{R}$
is the collection of sensorimotor data from all robots and t is time.
$S_{t}^{R}$ is defined as
$S_{t}^{R}=[S_{t}^{Robot_{1}},S_{t}^{Robot_{2}},S_{t}^{Robot_{3}},...,S_{t}^{Robot_{n}}]$,
where each element holds the sensorimotor information collected by each robot.
In this study, we used 2 robots whose morphologies and capabilities are
significantly different from each other to test the capabilities of our
system. We aim to learn correspondences between robots even when the means of
achieving the same task is significantly different between the two robots.

The system we propose connects multiple CNMP networks to find a common latent
space between them. At the start of each training iteration, a task trajectory
$I_{d}$ is sampled randomly from $I$. From $I$, an arbitrary number of data
points of $(t,S_{t}^{R})$ are sampled randomly as observations. Here, the
number of selected data points is between 1 and $obs_{max}$. $obs_{max}$ is a
hyperparameter that is used to limit the maximum number of observations used
in each training step. The set of sampled observations is defined as
$O^{R}=\\{(t_{i},S_{t_{i}}^{R})\\}_{i}^{obs_{max}}$ where
$(t_{i},S_{t_{i}}^{R})\in I_{d}$ . In our experiments, we used two robots, a
manipulator and a mobile robot. The observation set of these will be
represented as $O^{manipulator}$ and $O^{mobile}$. Other than $O^{R}$, a
target tuple $(t_{target},S_{t_{target}}^{R})$ is also sampled from $I_{d}$
which is used to learn the distributions on $t_{target}$ for observations
collected by all robots.

The purpose of the model is to find a common latent representation between
observations of all robots. To this end, the observations of each robot
$O^{r}$ are converted to latent representations, which are calculated using
the following equation:

$L_{i}^{r}=E^{r}((t_{i},S_{t_{i}}^{r})|\theta^{r})\quad(t_{i},S_{t_{i}}^{r})\in
O^{r},r\in R$ (1)

where $E^{r}$ is a deep encoder for robot $r$ with weights $\theta^{r}$, and
$L_{i}^{r}$ is the latent representation constructed using the given
observation. After the encoding process, we obtained $L_{i}^{manipulator}$ and
$L_{i}^{mobile}$ for the manipulator and the mobile robot, respectively.
Following the construction of $L_{i}^{r}$’s, these representations are
averaged for each robot:

$L^{r}=\frac{1}{n}\sum_{i}^{n}L_{i}^{r}\quad r\in R$ (2)

where n is the number of sampled observations during this iteration. Following
this equation, we calculated $L^{manipulator}$ and $L^{mobile}$, which
represent the latent representations constructed using $O^{manipulator}$ and
$O^{mobile}$, respectively. We aim to find a common latent representation
using the constructed latent representations to generate the required
trajectories for both robots. To accomplish this, a general latent
representation $L$ is constructed by taking the convex combination of all
$L^{r}$:

$L=\sum_{r}^{R}p^{r}L^{r}\quad 0\leq p^{r}\leq 1,\sum p^{r}=1,r\in R$ (3)

This calculation is applied to blend all latent representations into one
common latent space. For this calculation to be carried out, all $L^{r}$
should have the same dimension.

After all representations are merged into one, this representation is decoded
to obtain target distributions on $t_{target}$ for all robots:

$(\mu_{t_{target}}^{r},\sigma_{t_{target}}^{r})=Q^{r}((L,t_{target})|\phi^{r})\quad
r\in R$ (4)

where $Q^{r}$ is a deep decoder with weights $\phi^{r}$ that constructs
distributions with mean $\mu_{t_{target}}^{r}$ and variance
$\sigma_{t_{target}}^{r}$ for the robot $r$. For our robots,
$(\mu_{t_{target}}^{manipulator},\sigma_{t_{target}}^{manipulator})$ and
$(\mu_{t_{target}}^{mobile},\sigma_{t_{target}}^{mobile})$ are obtained using
$Q^{manipulator}$ and $Q^{mobile}$, respectively. The learning goal of our
system is to obtain more accurate distributions for given condition point(s)
and the target point, and therefore the loss is defined similarly to [14]:

$Loss=-\sum_{r}^{R}logP(S_{t_{target}^{r}}|\mu_{t_{target}}^{r},\sigma_{t_{target}}^{r})\quad
r\in R$ (5)

where $S_{t_{target}}^{r}\in S_{t_{target}}^{R}$ is the target sensorimotor
state for agent $r$.

After sufficient training, the system can make action predictions for all time
steps of that action for all robots when a distinct $O^{R}$ is provided to the
system. More importantly, when querying, trajectory generation can be achieved
by making $p^{r}$ for the desired robot zero. In our experiments, the system
was able to generate the desired trajectories for both the manipulator and the
mobile robot when observations from only one of them were provided to the
system.

Figure 3:  Training and test paths correspond to the paths required to achieve
the obstacle avoidance tasks in correspondence learning when the task requires
the same Cartesian path experiment (Section III-A). Left and right figures
show the manipulator and mobile robots. Blue and red paths were used for
training and testing, respectively. Green lines show the action range of the
manipulator.

Figure 4: The snapshots from the mobile robot to the manipulator (on top) and
from the manipulator to the mobile robot (on bottom) for the correspondence
learning when the task requires the same Cartesian path experiment (Section
III-A). Using three observations from the given trajectory on a new task
configuration (on the left of the figure), the system is able to generate the
desired trajectories for the other robot (on the right ).

Figure 5: Generated trajectories of one of the test sets of correspondence
learning when the task requires the same Cartesian path experiment (Section
III-A) with their ground truth values. Dashed lines are expected ground truth
values and solid lines represent the generated values.

##  III EXPERIMENTS

To analyze the capabilities and limits of our method, we designed three
experiments (Fig.1. Across all experiments, we used a mobile robot named
Pioneer 3-DX and a manipulator with six degrees of freedom named UR-10. The
mobile robot is controlled by setting the desired position and orientation of
the robot body at each time step. Therefore $S^{mobile}$ includes position and
orientation. The manipulator, on the other hand, is controlled by setting the
desired joint angles of the robot arm, therefore $S^{manipulator}$ includes
the six joint angles. The simulation experiments were conducted in the
CoppeliaSim environment [15].

###  III-A Correspondence learning when the task requires the same Cartesian
path

This experiment aims to verify that our method can learn the correspondences
between two robots when they need to follow the same path in task space. For
this, an obstacle avoidance task is designed where the robots are required to
avoid the obstacle using different paths. Eight different obstacle avoidance
trajectories were collected for each robot as shown in Fig. 3 where six
trajectories were used for training and the rest were used for testing. After
paired training is completed with a number of demonstration trajectories,
where position and orientation are used for the mobile robot and joint angles
are used for the manipulator robot, the model was queried to construct the
trajectories for a robot given a number of the observation points from the
other robot at certain time steps as input. When conditioned with an arbitrary
number of points collected by any robot, the system was able to construct the
corresponding trajectories. In other words, the manipulator/mobile robot was
able to avoid the obstacles in the test and training setups following the
trajectories generated based on task input from the other robot.

Fig. 4 provides a number of snapshots for task configurations that were not
observed during training. On the top and bottom rows, transfers from mobile
robot to manipulator and from manipulator to mobile robot are shown,
respectively. On the top left, the mobile robot is provided with a new
movement trajectory that was not seen during training. Our system uses
observations (position and orientation of the mobile robot) from this
trajectory at three different time steps and is able to generate the desired
movement trajectory (joint angles) for the manipulator robot, whose
corresponding motion is shown on the right.

We analyzed the results by comparing the trajectories generated by the model
and the desired trajectories. We applied this process both to training
trajectories and test trajectories. For both cases, we obtained results with
negligible error. Generated trajectories of one of the test cases can be seen
in Fig. 5 with the desired trajectories.

###  III-B Correspondence learning when different Cartesian paths are needed
for task completion

The morphologies and the capabilities of the mobile and manipulator robots are
different, and therefore they can achieve the same task using different means
depending on the task and the involved objects. In this section, we aim to
test whether our system can learn the correspondence between robots even if
the means to achieve the same task, i.e. the paths they follow, are different
(Fig. 1 (B)). As a use case, we selected a task that required the robots to
retrieve a cup to a desired location. In this setup, while the demonstrator
used one of the fingers of the gripper of the manipulator as a hook to pull
the cup to itself, the mobile robot is driven behind the cup to push it to the
same location. The location of the cup is changed, while the performed
strategy remains the same. Note that the cup is only chosen to represent the
goal and no information about the cup is used by either of the robots or by
any part of our model. To show the generalization capability of our system,
during testing, the cup is placed in a position that was not observed during
training. Seven demonstrations were collected for each robot. In order to
provide a clearer picture, we presented only three of them in Fig. 6.

Figure 6: A number of sample demonstration paths collected for cup retrieval
tasks of correspondence learning when different Cartesian paths are needed for
task completion experiment (Section III-B). The paths the robots need to
follow are different in the task space. The blue paths are for training and
the red paths are for testing. Green lines represent the action range of the
manipulator.

Fig. 7 provides a number of snapshots for tasks that were not observed during
training. The top and bottom rows show the transfers from mobile robot to
manipulator and from manipulator to mobile robot, respectively. As in the
previous experiment, given a number of observations from the trajectory of a
robot, our system was able to generate the trajectory for the other robot for
this new task configuration.

Figure 7: Snapshots from mobile robot to manipulator (on top) and from
manipulator to mobile robot (on bottom) for the cup retrieval tasks of
correspondence learning when different Cartesian paths are needed for task
completion experiment (Section III-B). Using three observations from the given
trajectory on a new task configuration (on the left of the figure), the system
is able to generate the desired trajectories for the other robot (on the right
).

Figure 8: The results of the cup retrieval tasks in the correspondence
learning when different Cartesian paths are needed for task completion
experiment in section III-B. Graphs show how as training progresses and
training loss decreases, the distance of the cup from the desired position
changes. Orange bars represent when the model is conditioned using
observations collected by the mobile robot and the blue bars represent when
conditioned with manipulator observations. Observations from the same time
steps are used in all runs (First, thirtieth, sixtieth, and the last points of
a sample with length 128).

This experiment was designed to retrieve the cup to a goal position. Therefore
the task success can be evaluated by comparing the goal position and the final
position of the object after action execution (even though no information
related to the object is used during control). We measured the error between
the desired and actual positions for the transfer in both directions at
different time steps during training and plot the error statistics along with
the training loss in Fig. 8. First of all, the results provide the mean and
variance of the final position errors obtained with 30 runs. Note that we only
present the results for the novel task configurations. As shown, the final
position error of the object drops with training and consistently with the
decrease in the loss of our model. The left plot shows that the object was
retrieved to the desired goal position with higher accuracy with the mobile
robot at the end of training. In the figure, the error bars shown with blue
color correspond to the error in the final object position when the new task
is demonstrated by the manipulator robot, and the ones with orange color
correspond to the error in the final object position when the new task is
demonstrated by the mobile robot. It can be seen in the figure that whether
the demonstrating robot is the same as or different from the observer robot
doesn’t affect the results significantly. Furthermore, this result shows that
our system achieved generalization in transfer in interpolation cases. On the
right plot, we also observe that the training loss drops similarly to the left
one, which shows that the system learned to generate the trajectories close to
the ground truth. However, the error in the final object position did not drop
as in the left part of the figure, and there is also a higher variance in
these error values. Upon observing the behavior, we noticed that minor
deviations in the manipulator’s trajectory could lead to failure to insert the
finger into the cup, thereby making it impossible to move the cup to the
desired position. This failure is not attributed to our correspondence
learning system because the same failure happens when the same robot is used
for observation and generation.

Figure 9: (A) Training set used in cup retrieval tasks in correspondence
learning when task completion requires different levels of complexity in the
SM spaces experiment (Section III-C). (B) One of the test results for the
manipulator in the same experiment. (C) A number of the paths the system may
choose to generalize in the same experiment. (D) All trajectories used in
training and testing of one of the joints of the manipulator in the same
experiment.

Figure 10: Mobile robot test results of the correspondence learning when task
completion requires different levels of complexity in the SM spaces experiment
(Section III-C). Dashed lines are expected ground truth values and solid lines
represent the generated values.

###  III-C Correspondence learning when task completion requires different
levels of complexity in the SM spaces

In the previous experiment, the robots needed to follow different paths in
order to achieve the same task, and gradual changes in the sensorimotor spaces
of both robots were required for the gradual change in the task space. In this
section, we aim to test our system when the trajectory of one robot needs to
change significantly in order to achieve gradually changing task
configurations (Fig.1 (C)). In order to address this challenge, we updated the
cup retrieval task by keeping the cup in the same position and moving it to
different goal positions by pulling/pushing the cup in different directions
(Fig.9). This way, we managed to increase the level of complexity of the
actions of the mobile robot while decreasing it for the manipulator. This
difference in complexity can be seen in panels (A) and (D) of Fig.9. Panel (A)
shows the paths followed by the mobile robot to move the cup to the desired
position, and panel (D) shows the joint angle trajectories for different
desired positions. As shown while there is a gradual change in the SM
trajectory of the manipulator, the change in the SM of the mobile robot is
very complex. Additionally, only the last parts of the joint angle
trajectories are different for all trajectories of all joints. The identical
part represents when the gripper is moving towards the cup and putting its
finger in the cup, and the rest represents pulling the cup to different
angles. We used seven trajectories for training and four trajectories for
testing for this section. Another objective of this experiment was to present
a proof-of-the-concept realization in the real world. Therefore, we used a
real manipulator, a real UR-10 with a 3-finger gripper, instead of its
gripper.

Figure 11: Snapshots from mobile robot to manipulator (on top) and from manipulator to mobile robot (on bottom) for the cup retrieval task of correspondence learning when task completion requires different levels of complexity in the SM spaces experiment (Section III-C). Using three observations from the given trajectory on a new task (on the left of the figure), the system is able to generate the desired trajectories for the other robot (on the right ). TABLE I: Directional and positional error values of all test cases of cup retrieval tasks in correspondence learning when task completion requires different levels of complexity in the SM spaces experiment (Section III-C) for the mobile robot. Dir. error is the difference between the angle the mobile robot pushed the cup to and the desired angle for that case. Pos. error is the difference between the desired and actual final position of the cup. Test | Positional and Angular Error
---|---
Case | Desired Angle | Outcome | Directional Error | Pos. Error
1 | 30 | 26.739 | 3.261 | 0.0226
2 | -60 | -62.082 | 2.082 | 0.0587
3 | -30 | -29.646 | 0.354 | 0.0514
4 | 60 | 63.176 | 3.176 | 0.0514

We observed that the difference in complexity created uncertainty in the
trajectory level for different task configurations between robots. The reason
is that the trajectories used to approach the cup until making contact with it
are the same for the manipulator independent of where the cup is being moved,
whereas the approach trajectory depends on the final goal position for the
mobile robot. We observed this affected the way the system generalizes for the
test cases. While the system can generate the desired trajectories close to
how we envisioned them before testing (Fig. 9 (B)) for the manipulator, this
was not the case for the mobile robot (Fig. 10 ). However, while there are
significant shifts between the desired trajectories and the generated ones in
Fig. 10, it can be seen that the last parts of the trajectories are close to
the desired ones. This can be explained by the uncertainty. While there is
only one way for the manipulator to achieve the task (with the same approach
trajectory), there are numerous ways for the mobile robot (with slightly
different approach trajectories). As seen in Fig. 9 (C), to push the cup to
the desired location, the system can generalize to any of the red paths shown
and to many others. It can be seen in Table I that although the system is
unable to generate the trajectories for the mobile robot as we predicted them,
it can successfully achieve the given task with small errors.

##  IV CONCLUSION

In this work, we proposed and realized a correspondence learning framework
that enables task-specific skill transfer between robots with significantly
different morphologies and allows goal-based imitation with different means
and robots. Our system was able to learn the task-level correspondence between
the pose space of a differential drive mobile robot and the joint space of a
fixed-based manipulator. The robots have different bodies and sensorimotor
spaces and may need to follow different paths to achieve the same task, yet,
our system successfully learned this correspondence. Even though our system
was shown to enable task-level skill transfer between manipulation and mobile
robots in cup retrieving tasks that require different movement trajectories,
the geometric or visual properties of the interacted objects were not used by
our system to make generalizations across tasks. In the future, we plan to
also focus on using object features, such as object images or high-level
features, in order to transfer and generalize object affordances [16] among
different robots.

## ACKNOWLEDGMENT

This research has been funded by the JST Moonshot R&D, Japan (Grant Number:
JPMJMS2292), by the JST CREST, Japan (Grant Number: JPMJCR21P4), by the JSPS
KAKENHI, Japan (Grant Number: 21H05053), by the World Premier International
Research Center Initiative (WPI), MEXT, Japan, by Japan Society for the
Promotion of Science, Grant-in-Aid for Scientific Research – the project with
number 22H03670, the project JPNP16007 commissioned by the New Energy and
Industrial Technology Development Organization (NEDO), and the Scientific and
Technological Research Council of Turkey (TUBITAK, 118E923). The authors would
like to thank Muhammet Hatipoglu for his help in accelerating the collection
of real-world data.

## References

  * [1] M. E. Taylor and P. Stone, “Transfer learning for reinforcement learning domains: A survey.” _Journal of Machine Learning Research_ , vol. 10, no. 7, 2009.
  * [2] M. E. Taylor, N. K. Jong, and P. Stone, “Transferring instances for model-based reinforcement learning,” in _Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2008, Antwerp, Belgium, September 15-19, 2008, Proceedings, Part II 19_. Springer, 2008, pp. 488–505.
  * [3] H. B. Ammar and M. E. Taylor, “Reinforcement learning transfer via common subspaces,” in _Adaptive and Learning Agents: International Workshop, ALA 2011, Held at AAMAS 2011, Taipei, Taiwan, May 2, 2011, Revised Selected Papers_. Springer, 2012, pp. 21–36.
  * [4] H. B. Ammar, E. Eaton, P. Ruvolo, and M. Taylor, “Unsupervised cross-domain transfer in policy gradient reinforcement learning via manifold alignment,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 29, no. 1, 2015.
  * [5] F. Liu, Z. Ling, T. Mu, and H. Su, “State alignment-based imitation learning,” _arXiv preprint arXiv:1911.10947_ , 2019.
  * [6] A. Gupta, C. Devin, Y. Liu, P. Abbeel, and S. Levine, “Learning invariant feature spaces to transfer skills with reinforcement learning,” _arXiv preprint arXiv:1703.02949_ , 2017.
  * [7] C. Devin, A. Gupta, T. Darrell, P. Abbeel, and S. Levine, “Learning modular neural network policies for multi-task and multi-robot transfer,” in _2017 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2017, pp. 2169–2176.
  * [8] D. Hejna, L. Pinto, and P. Abbeel, “Hierarchically decoupled imitation for morphological transfer,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 4159–4171.
  * [9] P. Sharma, D. Pathak, and A. Gupta, “Third-person visual imitation learning via decoupled hierarchical controller,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [10] Q. Zhang, T. Xiao, A. A. Efros, L. Pinto, and X. Wang, “Learning cross-domain correspondence for control with dynamics cycle-consistency,” _arXiv preprint arXiv:2012.09811_ , 2020.
  * [11] K. Kim, Y. Gu, J. Song, S. Zhao, and S. Ermon, “Domain adaptive imitation learning,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 5286–5295.
  * [12] M. Akbulut, E. Oztop, M. Y. Seker, X. Hh, A. Tekden, and E. Ugur, “Acnmp: Skill transfer and task extrapolation through learning from demonstration and reinforcement learning via representation sharing,” in _Conference on Robot Learning_. PMLR, 2021, pp. 1896–1907.
  * [13] M. Y. Seker, A. Ahmetoglu, Y. Nagai, M. Asada, E. Oztop, and E. Ugur, “Imitation and mirror systems in robots through deep modality blending networks,” _Neural Networks_ , vol. 146, pp. 22–35, 2022.
  * [14] M. Y. Seker, M. Imre, J. H. Piater, and E. Ugur, “Conditional neural movement primitives.” in _Robotics: Science and Systems_ , vol. 10, 2019.
  * [15] E. Rohmer, S. P. Singh, and M. Freese, “V-rep: A versatile and scalable robot simulation framework,” in _2013 IEEE/RSJ international conference on intelligent robots and systems_. IEEE, 2013, pp. 1321–1326.
  * [16] M. Cakmak, M. Dogar, E. Ugur, and E. Sahin, “Affordances as a framework for robot control,” in _Proceedings of the 7th international conference on epigenetic robotics, epirob’07_ , 2007.


11institutetext: Laboratoire d’Astrophysique de Bordeaux, Univ. Bordeaux,
CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France
22institutetext: Università di Padova, Dipartimento di Astronomia, vicolo
dell’Osservatorio 3, 35122 Padova, Italy 33institutetext: Observatoire
astronomique de l’Université de Genève, chemin Pegasi 51, 1290 Versoix,
Switzerland

# Toward a multidimensional analysis of transmission spectroscopy. Part III:
Modelling 2D effects in retrievals with TauREx ††thanks: The TauREx 2D plugin
is available at https://forge.oasu.u-bordeaux.fr/falco/taurex_2d

Tiziano Zingales 11 2 2 Aurélien Falco 11 William Pluriel 11 3 3 Jérémy
Leconte 11

New-generation spectrographs dedicated to the study of exoplanetary
atmospheres require a high accuracy in the atmospheric models to better
interpret the input spectra. Thanks to space missions like James Webb Space
Telescope (JWST), ARIEL and Twinkle, the observed spectra will indeed cover a
large wavelength range from visible to mid-infrared with an higher precision
compared to the old-generation instrumentation, revealing complex features
coming from different regions of the atmosphere. For hot and ultra hot
Jupiters (HJs and UHJs), the main source of complexity in the spectra comes
from thermal and chemical differences between the day and the night sides. In
this context, one-dimensional plane parallel retrieval models of atmospheres
may not be suitable to extract the complexity of such spectra. In addition,
Bayesian frameworks are computationally intensive and prevent us from using
complete three-dimensional self-consistent models to retrieve exoplanetary
atmospheres, and they constrain us to use simplified models to converge to a
set of atmospheric parameters. We thus propose the TauREx 2D retrieval code,
which uses two-dimensional atmospheric models as a good compromise between
computational cost and model accuracy to better infer exoplanetary atmospheric
characteristics for the hottest planets. TauREx 2D uses a 2D parametrization
across the limb which computes the transmission spectrum from an exoplanetary
atmosphere assuming azitmuthal symmetry. It also includes a thermal
dissociation model of various species. We demonstrate that, given an input
observation, TauREx 2D mitigates the biases between the retrieved atmospheric
parameters and the real atmospheric parameters. We also show that having a
prior knowledge on the link between local temperature and composition is
instrumental in inferring the temperature structure of the atmosphere.
Finally, we apply such a model on a synthetic spectrum computed from a GCM
simulation of WASP-121b and show how parameter biases can be removed when
using two-dimensional forward models across the limb.

##  1 Introduction

Exoplanets have diverse and complex atmospheres, which often require very
sophisticated methods to be correctly interpreted. In particular, hot Jupiters
(HJs) and ultra hot Jupiters (UHJs) appear to have a strong dichotomy between
the inflated hot day side and the cold and shrunken night side (Showman &
Guillot, 2002; Showman et al., 2008; Menou & Rauscher, 2009; Wordsworth et
al., 2011; Heng et al., 2011; Charnay et al., 2015; Kataria et al., 2016;
Drummond et al., 2016; Tan & Komacek, 2019; Pluriel et al., 2020, 2022). For
such planets, it is decisive to use more-than-one dimensional models to, at
least, reproduce realistically the different chemistry and physical
properties.

Transit spectroscopy allows the observer to infer atmospheric features by
looking at a limited volume around the terminator line (Brown, 2001;
Kreidberg, 2018). The transmitted light carries the information of both the
day and the night side of the planet of a small region around the terminator
line Fortney et al. (2010); Caldas et al. (2019); Wardenier et al. (2022). For
cold planets, a one-dimensional plane parallel geometry can be enough to
interpret low resolution data since the atmosphere is more homogeneous
(MacDonald et al., 2020). Nevertheless, in the case of HJs and UHJs, the use
of a higher-dimensional model becomes crucial to mitigate the biases between
the retrieved model and the real three dimensional observed atmosphere. The
impact of 3D atmospheric structures cannot be neglected and will have an
impact on the spectra thus on the retrieval analysis because of chemistry
heterogeneities (Changeat et al., 2019; Baeyens et al., 2021), or due to
morning-evening asymmetry (Line & Parmentier, 2016; Parmentier et al., 2016;
MacDonald et al., 2020; Espinoza & Jones, 2021).

Atmospheric retrievals encompass a series of algorithms which try to associate
an input exoplanetary spectrum to an atmospheric model (Rodgers, 2000; Irwin
et al., 2008; Madhusudhan & Seager, 2009; Line et al., 2013; Waldmann et al.,
2015b, a; Gandhi & Madhusudhan, 2018; Al-Refaie et al., 2019; Mollière et al.,
2019; Zhang et al., 2019; Min et al., 2020). Bayesian frameworks are very
useful to find the best atmospheric models which will reproduce the input
spectra. They are often the main analysis tool to infer atmospheric properties
of exoplanets. The use of such tools poses two main challenges:

  * •

Computational cost: Bayesian analysis is typically computationally-intensive,
especially when one tries to fit several model parameters (e.g. higher than
10). In order to reach convergence in a few hours timespan, it is necessary to
use parallelized-software on a multi-core environment.

  * •

Model accuracy: Real atmospheres are complex three-dimensional structures,
whose physics can be computationally-intensive to reproduce correctly.
Moreover, such complex 3D atmospheres often lead to degeneracies when using
too broad assumptions in the models.

Inverse modeling requires the fast generation of a high number of forward
models (e.g. $10^{3}-10^{7}$). Each forward model thus needs to be computed in
a very short amount of time to allow the Bayesian Framework to converge in a
reasonable amount of time. For this reason the retrieval community is obliged
to accept very strong assumptions on the atmospheric model to find a good
compromise between computational time and the atmospheric model’s precision,
given the accuracy of an input spectrum.
The Hubble/WFC3 camera represents the state-of-the-art technology for the
study of exoplanetary atmospheres, sensitive in the 1.1$\mu$m-1.7$\mu$m
wavelength range, with an accuracy down to $\sim$50 ppm. With the information
given by the Hubble/WFC3 camera, a 1D plane parallel atmosphere can be enough
to claim detection of several molecular species present in the exoplanetary
spectra, i.e. H2O, NH3, CH4, TiO and VO (Deming, 2010; MacDonald &
Madhusudhan, 2017; Molaverdikhani et al., 2019; Changeat & Edwards, 2021). It
is however challenging to get to elementary abundances i.e., O/H, C/O etc.,
only with the information from the WFC3 wavelength range. Today, there are
several studies which demonstrate the presence of water using Hubble data with
a precision up to 0.3 dex (Wakeford et al., 2018; Tsiaras et al., 2018;
Welbanks et al., 2019; Barstow et al., 2020).

Future space missions, like ARIEL (Tinetti et al., 2016, 2018), Twinkle
(Tessenyi et al., 2016) and JWST (Gardner et al., 2006), will spend part or
most of their observation time on exoplanets. They will have a wider
wavelength range as well as a higher spectral resolution and spectral accuracy
(Beichman et al., 2014). Planets like HJs and UHJs, which will be observed by
these future space missions cannot be well interpreted by merely assuming a 1D
plane parallel atmosphere as a forward model (Pluriel et al., 2020). More
complex atmospheric models will be necessary to explain some spectral features
correctly, and a compromise between computational cost and model accuracy will
still be required. There has been several attempts to address these
challenges, taking into considerations, for example, phase curves retrievals
(Feng et al., 2016; Changeat & Al-Refaie, 2020; Feng et al., 2020; Taylor et
al., 2020; Lacy & Burrows, 2020; MacDonald et al., 2020; MacDonald & Lewis,
2022; Espinoza & Jones, 2021; Mikal-Evans et al., 2022).
All of these factors lead us to reconsider the use of one dimensional models
and to reconstruct the forward model used in the TauREx atmospheric retrieval
tool (Waldmann et al., 2015b, a; Al-Refaie et al., 2021) considering an
additional geometric dimension. More-than-1D models and retrievals have
recently been discussed by Lacy & Burrows (2020); MacDonald & Lewis (2022) and
very recently by Nixon & Madhusudhan (2022). In particular, the approach of
Lacy & Burrows (2020) relies on a similar 2D atmospheric structure, based on
Caldas et al. (2019). However, their chemical approach differs from our
parametrization. They consider chemical equilibrium everywhere and retrieve
the bulk metallicity. In TauREx 2D, we can use two modes (see Sec. 2.1); one
with more freedom: the abundances in both hemispheres are retrieved
separately. The second mode assumes a chemical equilibrium, taking into
account thermal dissociation, and we retrieve the deep abundance of each
species. Another difference between our approach and that of Lacy & Burrows
(2020) is our use of a deep temperature parameter set for pressures higher
than a threshold. These 2 parameters, based on hot Jupiter GCM simulations,
correspond to the level where the radiative time scale becomes higher than the
dynamics time scale e.g. to an homogeneous temperature layer.

In the following sections, we propose a two-dimensional atmospheric model as a
good compromise between model accuracy and computational requirement for a
Bayesian analysis. The model, discussed in its forward usage in Falco et al.
(2022), is tested here within a Bayesian framework. We show the advantages and
the weaknesses of using a more-than-one dimensional model through a benchmark
of three atmospheric configurations with increasing complexity. Finally, we
investigate the atmosphere of an ultra hot Jupiter, following the results of
Pluriel et al. (2020), in presence of molecular dissociation, to quantify the
intrinsic model biases introduced within a Bayesian framework. One of the main
characteristics of TauREx 2D is the introduction of compositional variation
between day and night side produced by thermal dissociation of atmospheric
molecules.

##  2 Two-dimensional Model

The forward model used in TauREx 2D is formally described in Falco et al.
(2022). We give here a short presentation of the model and refer to the above
paper for more information.

The model relies on the assumption that the atmosphere is symmetrical around
the star-observer axis. The computation of the transmission spectrum therefore
reduces to a two-dimensional problem.

A 1D pressure profile is first defined within the assumption of a plane-
parallel atmosphere, in a fashion similar to TauREx (Al-Refaie et al., 2021).
The temperature is then defined in a 2D coordinate system $(\alpha^{*},P)$ of
which one dimension is angular and the second is pressure-based. This
coordinate system is shown by the underlying (black) grid in Fig. 1.

Figure 1: The red-colored grid is a polar grid $(\alpha^{*},z)$ of which the
radial ticks are set as the altitude corresponding to levels of pressure in
the day side of the planet. The input pressure-based grid $(\alpha^{*},P)$ is
also represented in the same coordinate system (colored cells from yellow to
purple).

The temperature, defined by Eq. 1, follows a linear transition between
$-\beta/2$ and $\beta/2$ as proposed in Caldas et al. (2019); Wardenier et al.
(2022).

$\left\\{\begin{array}[]{lll}P>P_{iso},&T=&T_{deep},\\\
\hbox{\multirowsetup$P<P_{iso},\left\\{\begin{array}[]{l}2\alpha^{*}\geq\beta,\\\
2\alpha^{*}\leq-\beta,\\\
-\beta<2\alpha^{*}<\beta,\end{array}\right.$}&T=&T_{day},\\\ &T=&T_{night},\\\
&T=&T_{night}\\\
&&+(T_{day}-T_{night})\frac{\alpha^{*}+\beta/2}{\beta}.\end{array}\right.$ (1)

It relies on three temperature variables ($T_{day}$, $T_{night}$, $T_{deep}$),
an angle parameter ($\beta$) and a pressure level defining the upper limit of
an isothermal annulus ($P_{iso}$). This parametrization, although simplistic,
is based on global climate model (GCM) simulations (Showman et al., 2015;
Parmentier et al., 2018; Tan & Komacek, 2019) to capture the most salient
features of the thermal structure of HJs.

The temperature is shown by the colors (yellow to purple) in Fig. 1. The data
(pressure, temperature, abundances) in the input grid $(\alpha^{*},P)$ is then
interpolated in an altitude-based coordinate system $(\alpha^{*},z)$, shown by
the red grid in Fig. 1. Using the new coordinate system $(\alpha^{*},z)$, rays
of light crossing the atmosphere may be subdivided into a list of segments, as
illustrated in Falco et al. (2022), using simple geometric formula. There is
one ray per layer in the model. The length of a ray $r$ which crosses the
$i$-th atmospheric layer is identified with $\Delta\ell_{r,i}$. The optical
depth can then be computed for each ray $r$ with:

$\tau_{r}^{\lambda}=\sum_{\begin{subarray}{c}i\end{subarray}}\frac{P_{r,i}}{k_{B}T_{r,i}}\left(\sum_{m=1}^{N_{gas}}\chi_{m,r,i}\sigma_{m,\lambda}+\sum_{j=1}^{N_{con}}k_{mie,j}\right)\Delta\ell_{r,i},$
(2)

where $P_{r,i}$ and $T_{r,i}$ are the pressure and temperatures of the cell
corresponding to the segment $i$ of the ray $r$, $k_{B}$ is the Boltzmann
constant, $\chi_{m,r,i}$ is the volume-mixing ratio of the $m$-th molecule,
$\sigma_{m,\lambda}$ is the total cross-section of Rayleigh scattering and
molecular and continuum absorptions, and finally, $k_{mie,j}$ is the cross
section associated with the Mie scattering for the $j$-th aerosol. $N_{gas}$
is the number of chemical species, $N_{con}$ is the number of aerosols, We
take into account collision-induced absorption for H2-H2 and H2-He pairs from
HITRAN (Richard et al., 2012). Formally, the contribution of CIA opacities is
normalized to the number of H2 molecules and included in
$\sigma_{\mathrm{H}_{2},\lambda}$.

The volume mixing ratios for all the species can also be differentiated with a
day and a night value arrays. In this work, we have set the number of layers
to 100 and the number of angular points (or slices) to 24. These two numbers
are a good compromise between precision of the model and computational speed
(Falco et al., 2022). Two angular points correspond to the full day and night
sides, while the rest equally subdivide the space between $-\beta/2$ and
$\beta/2$, accounting for the linear transition in-between.

###  2.1 Parameters

Table 1: Input parameters for the three atmosphere configurations.

Test | $T_{\text{deep}}$ | $T_{\text{up, day side}}$ | $T_{\text{up, night side}}$ | $P_{\text{iso}}$ | log10[H2O]∗ | log10[CO] | $\beta$
---|---|---|---|---|---|---|---
1D Fully Isothermal | 2500K | 2500K | 2500K |  $10^{4}$Pa | -3.30 | -3.36 | -
Two Layers | 1800K | 2600K | 2600K |  $10^{4}$Pa | -3.30 | -3.36 | -
Day/Night | 2500K | 3300K | 500K |  $10^{4}$Pa | -3.30 | -3.36 | 20∘

∗These values change as a function of pressure and temperature when water
dissociation is assumed.

We tested TauREx 2D retrievals using three different configurations described
in Tab 1. In all three configurations the thermal structure of the atmosphere
is described using the following free parameters:

  1. 1.

$\beta$ angle: the angle across the terminator line within which we compute
the radiative transfer equation; s

  2. 2.

$P_{iso}$: the pressure of the boundary between the lower and the upper level;

  3. 3.

$R_{p}$: the 10 bar radius of the planet;

  4. 4.

$T_{day/night/deep}$: the temperature in the day and night side of the planet
and in the deeper isothermal region, respectively;

  5. 5.

$P_{clouds}$: the pressure of the upper level of the cloud deck.

The TauREx 2D retrieval code is configured to have a deep atmospheric region
where the molecular abundances and the temperature remain constant: at
pressures higher than $P_{iso}$, the temperature is equal $T_{deep}$ and the
volume mixing ratio of each molecule has a constant value. In the upper part
of the atmosphere, where $P<P_{iso}$, the temperature is different between the
day side ($T_{day}$) and the night side ($T_{night}$).

For the chemistry, we can use two different schemes:

  1. 1.

Constant (or “free”) chemistry mode, where volume mixing ratios are constant
within each column ($\alpha$) for $P<P_{iso}$. In this mode, TauREx 2D can fit
independently the abundances on each hemisphere, this implies three free
parameters per molecule ($[mol]_{side1}$, $[mol]_{side2}$ and $[mol]_{deep}$).

  2. 2.

Thermal dissociation mode, where we assume local equilibrium chemistry. In
this mode, the free parameters are the deep abundances of each molecule
($[mol]$), that are considered uniform over the planet but we recompute the
local abundances as a function of the local pressure and temperature following
Parmentier et al. (2018). In the regime we consider, these changes in
abundances are mainly driven by thermal dissociation.

The thermal dissociation module we used in this work follows the
parametrization of Parmentier et al. (2018), which uses a modified version of
the NASA CEA Gibbs minimization code (Gordon et al., 1984) in a grid of models
to investigate gas and condensate equilibrium chemistry in differrent
atmospheric conditions in substellar objects (Moses et al., 2013; Kataria et
al., 2015; Skemer et al., 2016; Wakeford et al., 2017; Burningham et al.,
2017; Marley et al., 2017).

While the constant chemistry mode could seem better as it offers more
flexibility, we will show that its higher number of parameters creates more
degeneracies than the other parametrization (see Section 4) so that the
thermal dissociation mode will be our default except when stated otherwise.

In the tests we show in the next sections, we assumed a constant 30 ppm error
bars in the considered wavelength range $0.5\mu$m-$10\mu$m. The assumed error
was kept throughout the whole spectral range with a normal distribution
(Greene et al., 2016). This value was used as a noise floor in Pluriel et al.
(2022, 2020) and it is a good representation of the average noise computed
using the estimated mean number of photoelectrons in Cowan et al. (2015) for a
WASP-121b like star at 270 pc during a single transit. This wavelength range
will be covered using MIRI and NIRSpec instruments on board of the JWST space
telescope (Gardner et al., 2006).

##  3 Test and Validation

In all the test cases, we assume the same planetary and stellar parameters,
that of the system WASP-121, whose values are specified in Tab 2.

Table 2: Input planetary and stellar parameters. Planetary parameters
---
$R_{p}$ | $1.807R_{Jupiter}$
$M_{p}$ | $1.183M_{Jupiter}$
$g_{s}$ | $9.39ms^{-2}$
$H_{e}/H_{2}$ | 0.25884514
$[\log{P_{min}},\log{P_{max}}]$ | $[-4,6](Pa)$
Stellar Parameters
$R_{*}$ | $1.458R_{\odot}$
$T_{eff}$ | $6460K$
Retrieval Bounds
$\beta$ (degrees) | $0-30$
$T$ (K) | $100-3600$
$\log{P}_{iso}$ (Pa) | -4 $-$ 6
$R_{p}$ ($R_{\textnormal{jupiter}}$) | 0.9 $-$ 2.7
$\log$ mixing ratios | -10 $-$-1
$\log{P}_{clouds}$ (Pa) | -4 $-$ 6

TauREx 2D works within the same Bayesian framework as the one dimensional
TauREx retrieval code (Al-Refaie et al., 2021). The atmospheric retrieval
which uses a two-dimensional forward model is more computationally intensive
than the one calculated with a one-dimensional forward model, due to its
additional number of parameters and geometrical dimensions.

Before focusing on a realistic test case, computed using SPARC/MITgcm (Showman
et al., 2009), it is fundamental to test the reliability of the TauREx 2D code
on some synthetic spectra generated by itself. We decided to test the code
using the three cases described in Tab 1. For these tests, the input spectra
have been generated by the TauREx 2D forward model itself. The three cases
represent three different scenarios with different PT profiles and, then H2O
mixing ratio profiles as shown in Fig. 2. All the simulated spectra have been
generated assuming a resolution $R=\lambda_{0}/\Delta\lambda=100$ where
$\lambda_{0}=5.25\mu m$ and $\Delta\lambda=[0.5\mu m-10\mu m]$.

Figure 2: Pressure Temperature (left) and $\log{[\text{H}_{2}\text{O}]}$
mixing ratio (right) profiles for the three different test cases: FI: Fully
Isothermal, TL: Two Layers and DN: Day Night cases. The DN case has been split
into three curves showing the profiles at the therminator, the day and the
night side.

The retrieval tests were performed using the Multinest algorithm (Feroz et
al., 2009) with 1000 live points and uniform prior distribution between the
bounds shown in Tab 2. In the next three subsections we discuss the results
for the three test configurations. For the sake of conciseness, the fitted
spectrum and posterior distributions will be shown only for the third, most
complex case. To summarize, the three tests show that TauREx 2D can perfectly
retrieve itself when there is enough information in the input spectrum.

###  3.1 Fully Isothermal atmosphere

A fully isothermal atmosphere represents the simplest assumption in our test
set. We assumed a constant vertical and day/night composition with H2O and CO
defined in Tab 1.

The solution found is consistent with the input model. TauREx 2D found a
solution which is consistent with the input parameters. TauREx found a volume
mixing ratio for $\log{[H_{2}O]}=-3.30^{+0.01}_{-0.01}$ for water and
$\log{[CO]}=-3.36^{+0.05}_{-0.06}$ and for CO, centered around the ground
truth for this case. Since we used a forward model with two vertical layers,
divided at pressure $P_{iso}$, TauREx found the best model where the upper
layer is suppressed at the top of the atmosphere around $10^{-3.07}$Pa and
treated the rest of the atmosphere as isothermal with $T=2502^{+10}_{-11}$,
consistently with the input model.

The 2D model of TauREx is equivalent to the 1D version of TauREx in the case
of isothermal atmospheres, as demonstrated in Falco et al. (2022). Therefore,
this retrieval test is a confirmation of the consistency of the implementation
with its 1D predecessor.

###  3.2 Two Vertical Layers atmosphere

The test case with two layers is an atmosphere where the lower layer has a
temperature $T_{deep}=1800K$ and the upper temperature has
$T_{upper,day/night}=2600K$ but we do not have any day/night difference. We
use the deep abundances defined in Tab 1 and let molecules dissociate as a
function of the pressure and temperature. This test is important to verify
whether our model can find the input vertical structure of the atmosphere,
without including day/night effects. When $T_{upper}\neq T_{deep}$, the code
should recover the input $P_{iso}$ and two different temperatures. When
$T_{upper}=T_{deep}$, it should converge to the isothermal configuration.

We chose these parameters for the two-layers configuration to demonstrate how
TauREx 2D can converge to the correct temperature structure, given that the
composition remains the same in both hemispheres and with a temperature
inversion between the lower and the upper part of the atmosphere. For such a
test, we put the pressure threshold $P_{iso}$ between the two layers at 1
mbar, which corresponds to the range of pressures probed in such hot
atmospheres, thus ensuring that the input spectrum should encode some
information about the temperature change.

TauREx 2D found a $T_{side1}=2500^{+700}_{-1400}$,
$T_{side2}=2600^{+600}_{-1400}$, $T_{deep}=2594^{+8}_{-12}$ and
$\log{P_{iso}}=-2.89^{+0.82}_{-0.51}$. TauREx favors a thin upper atmosphere
and fits the lower layer with an isothermal at $\sim 2600K$, which is the
input upper temperature. This means that we managed to probe only the upper
layer of the atmosphere. The retrieved molecular abundances found are
compatible with the input values $\log{[H_{2}O]}=-3.33^{+0.01}_{-0.01}$ for
water and $\log{[CO]}=-3.34^{+0.05}_{-0.04}$.

###  3.3 Day/night atmosphere

With this complex test, we explore the most complex atmospheric configuration
in our test cases. Indeed, the input spectrum was generated assuming an
isothermal inner annulus at $T_{deep}=2500K$ a day side temperature
$T_{day}=3300K$ and a night side of $T_{night}=500K$ with a $\beta$ angle of
20∘ (see Tab 1 for the exact setting). The day/night thermal difference is
realistic according to GCM models of this type of planet (e.g. WASP-121b) and,
within this scenario, the $\beta$ parameter can give us a hint of the 3D
structure of the planet, since within this angle we probe both, the day side
and the night side, around the terminator line. Based on this GCM equatorial
temperature maps, the atmosphere can be approximated to an isothermal day and
night sides, with a sharp transition region between the two sides (represented
by the $\beta$ parameter), and a isothermal ring under a certain pressure.
Within this scenario, the parameter can give us a hint as to the 3D structure
of the planet, since within this angle we probe both, the day side and the
night side, around the terminator line. In Fig. 3 we show the best fit model
for this last test case.

Figure 3: Best fit spectrum for the Day/Night convergence test.

In Fig. 4 we can see how TauREx finds a unique solution around the ground
truth represented by the red vertical lines.

Figure 4: Posterior distribution for the Day/Night case. The input spectrum
refers to an hot atmosphere with three different equilibrium temperatures:
lower temperatures, upper day and upper night. As in the previous models, H2O
and H2 can dissociate as a function of P (Pa) and T (K). The red lines
represent the input values to generate the input spectrum.

With this final convergence test, we can see the reliability of TauREx 2D with
handling spectra from bidimensional models. We let the $P_{iso}$ and $\beta$
as free parameters and, we give TauREx the possibility to look for the best
spectrum.

The $\beta$ parameter retrieved with TauREx 2D, in particular, would give us
valuable information to model three dimensional effects of an exoplanetary
atmospheres. We see in Fig. 4 that TauREx found the parameters of the input
spectrum.

##  4 Intrinsic degeneracies in transmission spectra: can we recover
independently the compositions of the day and night sides?

By construction, the transmission spectrum of a planet only depends on the
total optical depth of the atmosphere along the various rays of light passing
through the limbs (see Eq. 2). Therefore, it can be seen right away that the
transmission spectrum of a planet would remain unaffected if we took a mirror
image of the planet, putting the day side away from the star.

This intrinsic symmetry of the problem tells us that, technically, a retrieval
algorithm cannot determine the conditions on the day and night sides. At best,
it can infer the presence of different conditions on both hemispheres and we
have to use our physical insight to decipher what conditions correspond to
what side; for example be saying that the hottest hemisphere is probably the
dayside one. This is why we use $T_{side1}$ and $T_{side2}$ as independent
parameters in our retrievals.

But the fact that we measure only the total optical depth creates another
degeneracy: to first order, only the total column density of a species along
the ray matters. To take an extreme example, if we had a completely isothermal
planet with two hemispheres having different chemical compositions, the effect
of a trace molecule (mol) on the spectrum would be unaffected as long as

$[\mathrm{mol}]_{day}+[\mathrm{mol}]_{night}=cst.$ (3)

A homogeneous atmosphere would have the same signature as an atmosphere where
all the molecules would be in a single hemisphere.

To check the performance of TauREx 2D with a fully isothermal atmosphere and a
compositional difference between the day side and the night side, we simulated
an input synthetic spectrum as in the fully isothermal test case, with the
difference that the water abundance changes from the day
($\log[H_{2}O]_{day}=-4.3$) to the night side ($\log[H_{2}O]_{night}=-3.3$).
As shown in the posterior distributions (Fig. 5), TauREx 2D favors the
solution in which the upper layer is flattened towards the upper edge of the
planetary atmosphere, in a small annulus with a pressure lower than
$\log[P](Pa)=-2.81^{+1.17}_{-0.80}$ and fits an isothermal atmosphere using
only the lower layer. The input temperature is found within 2-sigma and the
lower water abundance $\log[H_{2}O]_{deep}=-3.53^{+0.04}_{-0.05}$ is within
1-sigma from the arithmetic average between the input abundances
$\log[H_{2}O]_{deep,input}=-3.56$.

Figure 5: Posterior distribution for a particular case of the Fully isothermal
atmosphere with different water abundances in the day
($\log[H_{2}O]_{day}=-5.3$) and the night side (($\log[H_{2}O]_{night}=-4.3$))

Of course, having different temperatures on both sides formally breaks this
degeneracy, because of the scale height effect as well as the temperature
dependence of the molecular opacities (MacDonald & Lewis, 2022). But the
question remains whether there is enough information in a transmission
spectrum to retrieve the composition of both sides independently.

To test this, we performed another retrieval on the ‘Day/night” input spectrum
used in sect. 3.3 using TauREx 2D’s ability to model the planet with
independent compositions for the day and night sides as described in Section
2.1 (”free” mode). The posterior distributions for this retrieval are shown in
Fig. 6, and the best fit spectra found by the two models (along with the
residuals) are shown in Fig. 7.

As can be seen in the bottom panel of Fig. 7, the residuals are well below the
30 ppm noise floor used as uncertainty in these retrievals. The two approaches
thus give an acceptable fit. However, it can be seen in Fig. 6 that the
posterior distributions of the ”free” retrieval show a lot of cross-
correlation between several parameters. In particular, as can be expected from
the argument above, when the CO or H2O abundances are high on one side, it can
be very small on the other. This model also seems to miss the temperature
dichotomy between the two hemispheres of the planet.

To see whether these two forward models can be statistically differentiated,
we introduce the (logarithmic) Bayes factor:

$\mathscr{B}=\log{\frac{E_{\mathrm{d}}}{E_{\mathrm{f}}}}=\log{E_{\mathrm{d}}}-\log{E_{\mathrm{f}}},$
(4)

where $E_{\mathrm{d}}$ and $E_{\mathrm{f}}$ are the Bayesian evidence for the
retrievals with dissociation (d) and free (f) models.

Given that the Bayesian evidence of the retrieval which uses thermal
dissociation and coupled day/night sides is $E_{\mathrm{d}}=2914$ while the
one related to the model which uses independent day/night sides is
$E_{\mathrm{f}}=2908$, it follows that the Bayes factor is $\sim$6 giving a
slight but significant advantage to the coupled model. In this case, it seems
that the higher number of free parameters penalizes the model with ”free”
composition.

This seems to show that there is not enough information in the input spectrum
(at least with the level of uncertainty expected for a JWST like observation)
to retrieve simultaneously and independently the composition of both sides of
the planet. Going further, the results even imply that the information on the
temperature dichotomy in a transmission spectrum mainly comes from the link
between temperature and composition in the observed atmosphere. Assuming
chemical equilibrium brings a considerable amount of prior knowledge on the
atmosphere. If this link seems robust in Ultra hot Jupiters where molecular
dissociation is fast enough to happen close to chemical equilibrium, it should
be revisited for cooler atmospheres. For this reason, the ability to retrieve
the temperature dichotomy on cooler planets might be greatly reduced compared
to what is shown in Lacy & Burrows (2020), where they assumed chemical
equilibrium to hold down to dayside temperatures near 1200 K. This issue
should be thoroughly investigated.

For all these reasons, in the following, we will keep using the model assuming
thermal dissociation at chemical equilibrium as it possesses less free
parameters and results in less correlated posterior distributions between the
parameters. We should however keep in mind that this assumption should be
revisited when modelling cooler atmospheres.

Figure 6: Posterior distribution for a retrieval test made using the same
”Day/Night” input spectrum as in Fig. 4 but retrieved using a forward
atmospheric model with independent compositions for the day and night sides.
One can see the strong cross-correlations between the abundances of a given
molecule on the day and night sides.  Figure 7: Spectral comparison between
the two Day/Night retrieval cases. We can see that both the coupled model (the
model that considers thermal dissociation in the atmosphere) and the “free”
model (which is the one that considers a constant chemistry, but independent
on both sides of the planets) lead to a consistent spectrum.

##  5 A GCM retrieval

After validating the TauREx 2D retrieval code on some self-generated spectra,
we now test its abilities on a more realistic atmosphere, generated by the
SPARC/MIT GCM and described in detail in Parmentier et al. (2018). We
simulated the transmission spectra using 30ppm photon noise as described in
section 2.1. We used the temperature map obtained by the GCM simulation and
compute the radiative transfer equation through a hot Jupiter atmosphere with
Pytmosph3R 2.0 (Falco et al., 2022). The representative planet for this test
is a WASP-121 b model, where molecular species can thermally dissociate as a
function of the pressure and temperature. In this model, we include TiO and VO
in the atmosphere as they have been detected on WASP-121 b by Evans et al.
(2018). These molecules are optical absorbers that heat the day side of the
planet and create a strong thermal inversion (Fortney et al., 2008; Parmentier
et al., 2015, 2018). Other molecules, ions and metals (such as Iron detected
in the UHJ WASP-76 b by Ehrenreich et al. (2020)) also contribute
significantly to thermal inversion in UHJs, as pointed out by Lothringer et
al. (2018). Based on this GCM model of WASP-121 b, we simulate two spectra
from two different atmospheric configurations to test its impact on the
retrieval analysis:

  * •

A solar metallicity atmosphere with H2O and CO as trace gases;

  * •

A solar metallicity atmosphere with H2O, CO, TiO and VO as trace gases.

The presence of TiO and VO covers indeed some fundamental wavelengths where
Rayleigh scattering otherwise dominates. In the absence of clouds, this
wavelength range allows the model to calibrate the pressure in the atmosphere
and thus to have better absolute abundances. We thus want to test whether
hiding these regions leads to additional biases on the retrieved parameters.

A crucial parameter related to the three-dimensionality of the planet is the
$\beta$ parameter. The $\beta$ angle is a 2D parameter and therefore its
equivalent for a 3D GCM simulation is hard to define. For the following input
spectra, the best approximation we can make is $5_{-3}^{+3}$ degrees based on
the GCM equatorial temperature map, indicating a very sharp transition between
the day and the night sides.

###  5.1 Atmosphere without optical absorbers

We start analyzing an atmosphere without the presence of optical absorbers
(e.g. molecules like TiO and VO). The retrieval performed by TauREx 2D
converged to the best-fit solution shown in Fig. 8. In this test we can also
see a degenerate solution converging to two different $\beta$ angles and
temperatures for the two sides of the planet.

Figure 8: Best fit model for an atmosphere without optical absorbers. TauREx
2D retrieval code found two possible solutions for this case. The comparison
with the retrieval made with a one-dimensional isothermal forward model shows
that TauREx 2D leads to a better interpretation of the atmosphere.

The two solutions found by TauREx have different physical properties and the
best fit parameters are shown in Fig. 9.

Figure 9: Posterior distribution for the case where we used as input a
WASP-121b GCM’s synthetic atmosphere using only H2O and CO as trace gases. The
red lines correspond to the input deep abundances from the GCM model
atmosphere.

TauREx 2D managed to converge around a realistic solution which is compatible
with the ground truth, i.e. the parameters used to generate the input
spectrum. Indeed, as shown by the red line in Fig. 9, the deep atmosphere
trace gas parameter volume mixing ratios are compatible with the ground truth
of the input model within $1\sigma$ error bars for both solutions. In
addition, the temperatures obtained in both solutions reflects realistically
the ones of WASP-121b. This two-dimensional model is more realistic and
accurate than a simpler one-dimensional model, since it is able to catch day-
night differences that, as shown in Pluriel et al. (2020), are crucial in Hot
Jupiters like WASP-121b.

We find $\beta_{0}=1.2_{-0.9}^{+1.2}$ and $\beta_{1}=11.4_{-2.8}^{+2.5}$,
respectively, the $\beta$ angle for first and second solution, whose measures
are correlated with the temperature values found in the two sides,
respectively $T_{night}=190_{-60}^{+150}$ and $T_{day}=2000_{-100}^{+100}$ for
the first solution, $T_{night}=980_{-80}^{+80}$ and $T_{day}=2320_{-30}^{+30}$
for the second solution. Interestingly, both solutions give the same pressure
level and the same temperature for the deep layer ($P_{iso}\sim 10$ Pa and
$T_{deep}\sim 1780K$). It means that the retrieval can probe approximately up
to that atmospheric region, after which we approximate the atmosphere with an
isothermal annulus of $T_{deep}\sim 1780K$. The $\beta$ parameter gives us a
hint of the three-dimensional structure of the planet, since it spans above
the day side and the night side, and a dramatic change on this parameter can
lead to a dramatic change in the retrieved temperature on both sides of the
planet. In both solutions, the transition between the day and the night side
retrieved is sharp since it is less than $12^{\circ}$.

The mixing ratio retrieved by the 2D model and the PT profile approximates the
input profiles better than the 1D model as shown in figures 10 and 13. Both
solutions from the two-dimensional retrieval are closer to the input mixing
ratio profiles than the solution of the one-dimensional isothermal model. In
particular, TauREx 2D found the correct abundance of water in the night side
and the correct value of CO within 2$\sigma$ error.

Figure 10: Mixing ratio for H2O and CO for the GCM case without optical
absorbers. The dashed lines are the retrieved values while the continuous
lines represent the input value, In yellow the results from the 1D isothermal
retrieval. For this retrieval, the terminator and the night mixing ratios
converged to the same values.

To evaluate the contribution of a 2D retrieval model, we compared the solution
given by a two-dimensional model with the one provided with a 1D simple model,
where we also take into account molecular dissociation along the vertical
dimension of the atmosphere.

From Pluriel et al. (2022), we saw how one-dimensional models cannot interpret
the input spectrum of two-dimensional models. To quantify how good the 2D
model is compared to the one dimensional one, we can calculate the logarithmic
Bayes factor expressed in Eq 4.

The resulting logarithmic Bayes factor for the 2D model compared to the 1D one
is $\mathscr{B}=+89.2$, which favours the solution of the two dimensional
retrieval over the one dimensional one with a $\sigma$-significance higher
than $5\sigma$ (Trotta, 2008).

For this test case, without molecules hiding the spectral part with a
significant contribution from the Rayleigh scattering, the two dimensional
retrieval converges to the input spectrum compatibly with the input parameters
of the synthetic spectrum. An important goal of this result is to have managed
to unravel the biases observed in Pluriel et al. (2020) in the molecular
abundances of the atmospheric gases. In addition, the parametrisation of our
2D retrieval model seize a thermal structure of the atmosphere which is
compatible with 3D GCM models.

###  5.2 Atmosphere with optical absorbers

After studying the case of an atmosphere without TiO and VO as optical
absorbers, we want to measure how the solution of a spectral retrieval in
presence of optical absorbers is impacted. Given the degeneracy of the TiO and
VO and the Rayleigh scattering at short wavelengths and, depending on the
spectral resolution and wavelength range of the data set, it could be common
not be able to disentangle correctly the contribution of the Rayleigh
scattering from the one of the optical absorbers. The suggested Bayesian model
has 11 free parameters and run in $\sim 20$ days on 20 CPU cores using 1000
live points. The best fit models found by TauREx 2D is shown in Fig. 11.

Figure 11: Best fit model for a GCM’s synthetic atmosphere with TiO and VO as
optical absorbers. TauREx 2D retrieval code found two possible solutions for
this case (the red and the blue for the solution with, respectively,
$\beta=4.3$ and $\beta=10.0$ degrees). In yellow and green is shown the
solutions from the one-dimensional retrieval using, respectively, an
isothermal profile and a 4 PT point profile.

The results are similar to the previous retrieval: (i) TauREx 2D finds two
different solutions that provide the same atmospheric spectrum (ii) converging
around two beta angles, $\beta_{0}=4.3_{-0.9}^{+0.9}$ and
$\beta_{1}=10.0_{-0.9}^{+0.8}$, (iii) getting the expected deep abundance of
each molecules considered. As we increase the degrees of freedom of the model
(in our case, we have an additional spatial dimension compared to the 1D
atmospheric model), the degeneracy of the solution could be a normal
consequence that we need to handle. The $\beta_{0}$ angle found in this
simulation is $\sim 3$ degrees higher than the $\beta_{0}$ solution found in
the retrieval without optical absorbers. In presence of TiO and VO, we can
probe a more extensive part of the day side in the high atmosphere and, then,
TauREx finds a greater retrieved $\beta$ angle for this latter case.

From Fig. 11 it is difficult to understand the contribution of the Rayleigh
scattering under the TiO and VO absorption bands. In that spectral region,
also the H2-H2 CIA provides a continuum contribution, with a peak at around
$2.2\mu$m. Campared to the case without optical absorption, we have a
precision decrease from 0.02 dex to 0.10 dex, which it is related to the loss
of a clear Rayleigh slope signal.

Figure 12: Top and middle: Temperature model for the two different solutions
found in the GCM retrieval with optical absorbers. For visual reasons, the
size of the atmosphere has been doubled. The temperature structure on the top
represents the one of the first solution, the green one in Fig. 9. On the
middle there is the thermal structure corresponding to the second, violet
solution of Fig. 15. Bottom: Temperature map (equatorial cut) of WASP-121 b
GCM simulation.

In Fig. 12 we show the retrieved thermal structure for this simulation,
together with the equatorial cut view of the model. In Fig. 15, we can see
that the two-dimensional retrieval managed to converge around the input mixing
ratios for the four molecules. The temperatures retrieved are respectively
$T_{night}=780_{-80}^{+80}$K and $T_{day}=2730_{-40}^{+400}$K for the first
solution, $T_{night}=690_{-40}^{+40}$K and $T_{day}=3470_{-120}^{+80}$K for
the second solution and $T_{deep}\sim 1800$K for both solutions.

Also in this test, the volume mixing ratio retrieved by the 2D model and the
PT profile approximates better the input profiles than the 1D retrieval, as
shown in figures 13 and 14. TauREx 2D found two solutions, approximating the
input day and night temperatures and estimating the deep atmosphere around its
average temperature (Fig. 13). Moreover, TauREx 2D found a chemical
distribution of elements closer to the input distribution than the one found
with the one dimensional retrievals (Fig. 14).

Figure 13: PT profile retrieved compared to the input PT profile. In both
scenarios, with and without TiO/VO, has been used PT profile with temperature
inversion. In blue, purple and red are shown, respectively, the night,
terminator and day side PT profile of both the input (continuous line) and the
retrieved models (dashed lines). In the upper and lower part we show,
respectively, the first and the second solution from the two-dimensional
retrievals. In green and yellow are shown, respectively, the solution from a
one-dimensional retrieval using the 4 PT point profile (only for the GCM case
with optical absorbers, and the isothermal PT profile.

.

Figure 14: Mixing ratio for H2O, CO, TiO and VO in the GCM test retrieval. The
dashed lines are the values retrieved by TauREx while the continuous lines
represent the input value. In yellow the results from the 1D isothermal
retrieval and, in green, we show the solution found using a 4-point PT profile
forward model. For this retrieval, the terminator and the night mixing ratios
converged to the same values.

The crucial result of this test is that TauREx 2D manage to retrieve the
expected deep abundance of each molecule with a very good agreement to fit the
input data. TauREx 2D correctly separate the contribution of the Rayleigh
scattering from the one of the optical absorbers, unraveling also here the
biases explained in Pluriel et al. (2020).

Also in this case, a comparison with an equivalent 1D retrieval leads to a
Bayes factor $\mathscr{B}=+47.9$ in favor of the 2D retrieval.

Figure 15: Best fit model for an atmosphere with TiO and VO as optical
absorbers. TauREx 2D retrieval code found two possible solutions for this
case. The red lines correspond to the known deep input abundances in the GCM
model.

##  6 Comparison with a more complex 1D model retrieval

The 2D retrieval, in this case, helps us digging more on the three-dimensional
structure of the planet. We see that the two different solutions differ mainly
from two possible $\beta$ angles which lead to two possible values of
$T_{day}$. We remark that a higher $\beta$ angle results in a higher day-side
temperature.

To quantify the goodness of the fit, and see the advantages of a two-
dimensional retrieval over the ”classic” one-dimensional, we fitted the same
input spectrum with a one dimensional model and calculated the Bayes factor
between the two best-fit models.

We tested the one-dimensional model twice, assuming two different pressure-
temperature profiles: one isothermal atmosphere and one four-point PT profile
atmosphere.

Table 3: Comparison of models Retrieval type | Number of parameters | Log Bayesian Evidence | $\tilde{\chi}^{2}$ |  $\mathscr{B}$ with 1D iso |  $\mathscr{B}$ with 1D 4-point PT
---|---|---|---|---|---
2D sol0 | 11 | 2795.0 | 0.07 | 47.9 | 14.5
2D sol1 | 11 | 2795.0 | 0.07 | 47.9 | 14.5
1D 4 PT | 12 | 2780.5 | 0.18 | 33.4 | -
1D iso | 7 | 2747.1 | 0.45 | - | -33.4

In Tab 3 we can find a resume of all the comparison made between the 2D
retrieval including TiO and VO and the 1D retrievals. Together with the
Bayesian evidence, we also provide the reduced $\tilde{\chi}^{2}$ of each best
model which is defined as follows:

$\tilde{\chi}^{2}=\frac{1}{N-p}\sum_{i=1}^{N}\left(\frac{O_{i}-C_{i}}{\sigma_{i}}\right)^{2},$
(5)

Where $O$ and $C$ are, respectively the observed and the computed spectral
points, $\sigma$ is the uncertainty, $p$ is the number of parameters in the
model and $N$ is the total number of measurements. As described in Tab 3, the
2D models approximate better the input spectrum, returning the lowest
$\tilde{\chi}^{2}$ with respect with the 1D retrievals.

In Fig. 11 we can see that the four PT points retrieval leads to a better
retrieval than the one done using an isothermal atmosphere. The logarithmic
Bayes factor $\mathscr{B}=\log{E_{4pt}}-\log{E_{iso}}=33$, which confirms that
the four-points retrieval is statistically better than the isothermal
retrieval with a $\sigma$-significance higher than $5\sigma$.

A comparison between the two dimensional and the four-points retrieval leads
to a logarithmic Bayes factor of $\mathscr{B}=\log{E_{2D}}-\log{E_{4pt}}=14$
and, then the two-dimensional retrieval is better also than the four points
retrieval with a $\sigma$-significance higher than $5\sigma$.

The reduced $\tilde{\chi}^{2}$ of the four points retrieval tells us that the
atmospheric model found is compatible with the input spectrum.

All the molecular abundances appear to be overestimated, and then, biased
compared to the input parameters used to generate the input spectrum.

Figure 16: Posterior distributions related to the four-points retrieval.
Confirming the parameter biases found by Pluriel et al. (2022) .

The calculation of the Bayes factor allows us to quantify the goodness of the
2D retrieval compared to the one-dimensional one. Indeed, if we consider only
the reduced $\tilde{\chi}^{2}$ of the fit, we see that the single models are
compatible with the input spectra with a high significance. However, it is
impossible to get rid of the parameters biases when using the one-dimensional
model. The key-point of the calculation of the logarithmic Bayes factor is to
highlight the fundamental differences of the two model and the limits given by
their geometry. With 30ppm error bars it is, indeed, possible to reveal and
explain some spectral features by introducing an additional geometric
dimension on the atmospheric model.

Indeed, even though the one-dimensional retrieval converges toward a
statistically significant atmospheric model, the parameters found can differ
from the ground truth by several orders of magnitudes. In Fig. 16, as an
example, we show the posterior distributions for a retrieval of the GCM input
spectrum with optical absorbers, using a one-dimensional atmospheric model and
a four temperature-pressure points. The best model fits very well the input
spectrum, but the retrieved amount of atmospheric trace gases differ from the
ground truth up to 2 orders of magnitudes.

We see that, the logarithmic Bayes factor defined in Eq 4, that expresses the
difference of the 2D model compared with the 1D one is $\mathscr{B}=+47.9$,
which privileges the solution of the bi-dimensional model with a
$\sigma$-significance higher than $5\sigma$. Also in this case the 2D model
leads to a statistically better fit than the one-dimensional model.

The statistical significance of the 2D model is higher than that of the 1D
model, and it is less affected by the parameters biases, since it approximates
better the three-dimensional geometry of the input model.

The comparison with the one-dimensional retrieval emphasizes the importance of
running more-than-one dimensions retrieval models when using 30ppm noise input
spectra.

In the era of the JWST, ARIEL and the Twinkle space missions it will be
fundamental to use better models which rely on more complex geometrical
assumptions than the plane parallel atmospheres.

##  7 Conclusion

More-than-one dimensional retrievals will be crucial to better interpret the
input spectral observation of new generation instruments. Even though the
approximation of a plane-parallel atmosphere can lead to a statistically
acceptable fit, it inevitably leads to parameters biases that can be very far
away from the reality (Pluriel et al., 2020, 2022). This is especially true
for HJs and UHJs, where the great dichotomy between the day and the night side
of the planet may lead to a misinterpretation of the input observations,
without an accurate model.

We developed a 2D version of TauREx, which demonstrates that a two dimensional
forward model is a good compromise between model accuracy and computational
requirements. This model can solve the issues of the parameters biases induced
by a one dimensional model which is essential in the context of future mission
such like JWST, Twinkle or Ariel. Thanks to TauREx 2D we are able to get with
much higher confidence than with 1D retrieval models the abundances of the
species.

We showed how the model choice could be crucial to properly retrieve the
thermal structure of an exoplanetary atmosphere (in Section 4), due to the
intrinsic degeneracies of a high dimensional transmission model.

Finally, in Section 6 we demonstrated how a two-dimensional retrieval with an
isothermal atmosphere returns a statistically more significant solution than a
one-dimensional retrieval using a more complex 4-points pressure temperature
profile.

###### Acknowledgements.

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement n∘679030/WHIPLASH). TZi acknowledges the funding support from
Italian Space Agency (ASI) regulated by ‘Accordo ASI-INAF n. 2013-016-R.0 del
9 luglio 2013 e integrazione del 9 luglio 2015 CHEOPS Fasi A/B/C’.

## References

  * Al-Refaie et al. (2019) Al-Refaie, A. F., Changeat, Q., Waldmann, I. P., & Tinetti, G. 2019
  * Al-Refaie et al. (2021) Al-Refaie, A. F., Changeat, Q., Waldmann, I. P., & Tinetti, G. 2021, ApJ, 917, 37
  * Baeyens et al. (2021) Baeyens, R., Decin, L., Carone, L., et al. 2021, MNRAS, 505, 5603
  * Barstow et al. (2020) Barstow, J. K., Changeat, Q., Garland, R., et al. 2020, Monthly Notices of the Royal Astronomical Society, 493, 4884
  * Beichman et al. (2014) Beichman, C., Benneke, B., Knutson, H., et al. 2014, PASP, 126, 1134
  * Brown (2001) Brown, T. M. 2001, ApJ, 553, 1006
  * Burningham et al. (2017) Burningham, B., Marley, M. S., Line, M. R., et al. 2017, MNRAS, 470, 1177
  * Caldas et al. (2019) Caldas, A., Leconte, J., Selsis, F., et al. 2019, Astronomy & Astrophysics, 623, A161
  * Changeat & Al-Refaie (2020) Changeat, Q. & Al-Refaie, A. 2020, ApJ, 898, 155
  * Changeat & Edwards (2021) Changeat, Q. & Edwards, B. 2021, ApJ, 907, L22
  * Changeat et al. (2019) Changeat, Q., Edwards, B., Waldmann, I. P., & Tinetti, G. 2019, The Astrophysical Journal, 886, 39
  * Charnay et al. (2015) Charnay, B., Meadows, V., Misra, A., Leconte, J., & Arney, G. 2015, ApJ, 813, L1
  * Cowan et al. (2015) Cowan, N. B., Greene, T., Angerhausen, D., et al. 2015, PASP, 127, 311
  * Deming (2010) Deming, D. 2010, Nature, 468, 636
  * Drummond et al. (2016) Drummond, B., Tremblin, P., Baraffe, I., et al. 2016, A&A, 594, A69
  * Ehrenreich et al. (2020) Ehrenreich, D., Lovis, C., Allart, R., et al. 2020, arXiv e-prints, arXiv:2003.05528
  * Espinoza & Jones (2021) Espinoza, N. & Jones, K. 2021, AJ, 162, 165
  * Evans et al. (2018) Evans, T. M., Sing, D. K., Goyal, J. M., et al. 2018, AJ, 156, 283
  * Falco et al. (2022) Falco, A., Zingales, T., Pluriel, W., & Leconte, J. 2022, in , A41
  * Feng et al. (2020) Feng, Y. K., Line, M. R., & Fortney, J. J. 2020, AJ, 160, 137
  * Feng et al. (2016) Feng, Y. K., Line, M. R., Fortney, J. J., et al. 2016, ApJ, 829, 52
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
  * Fortney et al. (2008) Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419
  * Fortney et al. (2010) Fortney, J. J., Shabram, M., Showman, A. P., et al. 2010, ApJ, 709, 1396
  * Gandhi & Madhusudhan (2018) Gandhi, S. & Madhusudhan, N. 2018, MNRAS, 474, 271
  * Gardner et al. (2006) Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Science Reviews, 123, 485
  * Gordon et al. (1984) Gordon, S., McBride, B., & Zeleznik, F. J. 1984, Computer program for calculation of complex chemical equilibrium compositions and applications. Supplement 1: Transport properties
  * Greene et al. (2016) Greene, T. P., Line, M. R., Montero, C., et al. 2016, ApJ, 817, 17
  * Heng et al. (2011) Heng, K., Frierson, D. M. W., & Phillipps, P. J. 2011, MNRAS, 418, 2669
  * Irwin et al. (2008) Irwin, P. G. J., Teanby, N. A., de Kok, R., et al. 2008, J. Quant. Spec. Radiat. Transf., 109, 1136
  * Kataria et al. (2015) Kataria, T., Showman, A. P., Fortney, J. J., et al. 2015, ApJ, 801, 86
  * Kataria et al. (2016) Kataria, T., Sing, D. K., Lewis, N. K., et al. 2016, ApJ, 821, 9
  * Kreidberg (2018) Kreidberg, L. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte, 100
  * Lacy & Burrows (2020) Lacy, B. I. & Burrows, A. 2020, ApJ, 905, 131
  * Line & Parmentier (2016) Line, M. R. & Parmentier, V. 2016, ApJ, 820, 78
  * Line et al. (2013) Line, M. R., Wolf, A. S., Zhang, X., et al. 2013, ApJ, 775, 137
  * Lothringer et al. (2018) Lothringer, J. D., Barman, T., & Koskinen, T. 2018, ApJ, 866, 27
  * MacDonald et al. (2020) MacDonald, R. J., Goyal, J. M., & Lewis, N. K. 2020, ApJ, 893, L43
  * MacDonald & Lewis (2022) MacDonald, R. J. & Lewis, N. K. 2022, ApJ, 929, 20
  * MacDonald & Madhusudhan (2017) MacDonald, R. J. & Madhusudhan, N. 2017, ApJ, 850, L15
  * Madhusudhan & Seager (2009) Madhusudhan, N. & Seager, S. 2009, ApJ, 707, 24
  * Marley et al. (2017) Marley, M. S., Saumon, D., Fortney, J. J., et al. 2017, in American Astronomical Society Meeting Abstracts, Vol. 230, American Astronomical Society Meeting Abstracts #230, 315.07
  * Menou & Rauscher (2009) Menou, K. & Rauscher, E. 2009, ApJ, 700, 887
  * Mikal-Evans et al. (2022) Mikal-Evans, T., Sing, D. K., Barstow, J. K., et al. 2022, Nature Astronomy, 6, 471
  * Min et al. (2020) Min, M., Ormel, C. W., Chubb, K., Helling, C., & Kawashima, Y. 2020, A&A, 642, A28
  * Molaverdikhani et al. (2019) Molaverdikhani, K., Henning, T., & Mollière, P. 2019, ApJ, 883, 194
  * Mollière et al. (2019) Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67
  * Moses et al. (2013) Moses, J. I., Line, M. R., Visscher, C., et al. 2013, ApJ, 777, 34
  * Nixon & Madhusudhan (2022) Nixon, M. C. & Madhusudhan, N. 2022, arXiv e-prints, arXiv:2201.03532
  * Parmentier et al. (2016) Parmentier, V., Fortney, J. J., Showman, A. P., Morley, C., & Marley, M. S. 2016, ApJ, 828, 22
  * Parmentier et al. (2015) Parmentier, V., Guillot, T., Fortney, J. J., & Marley, M. S. 2015, A&A, 574, A35
  * Parmentier et al. (2018) Parmentier, V., Line, M. R., Bean, J. L., et al. 2018, Astronomy & Astrophysics, 617, A110
  * Pluriel et al. (2022) Pluriel, W., Leconte, J., Parmentier, V., et al. 2022, A&A, 658, A42
  * Pluriel et al. (2020) Pluriel, W., Zingales, T., Leconte, J., & Parmentier, V. 2020, A&A, 636, A66
  * Richard et al. (2012) Richard, C., Gordon, I. E., Rothman, L. S., et al. 2012, Journal of Quantitative Spectroscopy & Radiative Transfer, 113, 1276
  * Rodgers (2000) Rodgers, C. D. 2000, Inverse Methods for Atmospheric Sounding: Theory and Practice
  * Showman et al. (2008) Showman, A. P., Cooper, C. S., Fortney, J. J., & Marley, M. S. 2008, ApJ, 682, 559
  * Showman et al. (2009) Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564
  * Showman & Guillot (2002) Showman, A. P. & Guillot, T. 2002, A&A, 385, 166
  * Showman et al. (2015) Showman, A. P., Lewis, N. K., & Fortney, J. J. 2015, ApJ, 801, 95
  * Skemer et al. (2016) Skemer, A. J., Morley, C. V., Zimmerman, N. T., et al. 2016, ApJ, 817, 166
  * Tan & Komacek (2019) Tan, X. & Komacek, T. D. 2019, ApJ, 886, 26
  * Taylor et al. (2020) Taylor, J., Parmentier, V., Irwin, P. G. J., et al. 2020, MNRAS, 493, 4342
  * Tessenyi et al. (2016) Tessenyi, M., Savini, G., Tinetti, G., et al. 2016, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 48, AAS/Division for Planetary Sciences Meeting Abstracts #48, 123.33
  * Tinetti et al. (2018) Tinetti, G., Drossart, P., Eccleston, P., et al. 2018, Experimental Astronomy, 46, 135
  * Tinetti et al. (2016) Tinetti, G., Drossart, P., Eccleston, P., et al. 2016, in Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave, ed. H. A. MacEwen, G. G. Fazio, M. Lystrup, N. Batalha, N. Siegler, & E. C. Tong, Vol. 9904, International Society for Optics and Photonics (SPIE), 658 – 667
  * Trotta (2008) Trotta, R. 2008, Contemporary Physics, 49, 71
  * Tsiaras et al. (2018) Tsiaras, A., Waldmann, I. P., Zingales, T., et al. 2018, AJ, 155, 156
  * Wakeford et al. (2018) Wakeford, H. R., Sing, D. K., Deming, D., et al. 2018, AJ, 155, 29
  * Wakeford et al. (2017) Wakeford, H. R., Visscher, C., Lewis, N. K., et al. 2017, MNRAS, 464, 4247
  * Waldmann et al. (2015a) Waldmann, I. P., Rocchetto, M., Tinetti, G., et al. 2015a, ApJ, 813, 13
  * Waldmann et al. (2015b) Waldmann, I. P., Tinetti, G., Rocchetto, M., et al. 2015b, ApJ, 802, 107
  * Wardenier et al. (2022) Wardenier, J. P., Parmentier, V., & Lee, E. K. H. 2022, MNRAS, 510, 620
  * Welbanks et al. (2019) Welbanks, L., Madhusudhan, N., Allard, N. F., et al. 2019, ApJ, 887, L20
  * Wordsworth et al. (2011) Wordsworth, R. D., Forget, F., Selsis, F., et al. 2011, ApJ, 733, L48
  * Zhang et al. (2019) Zhang, M., Chachan, Y., Kempton, E. M. R., & Knutson, H. A. 2019, PASP, 131, 034501