# Efficient Search and Learning for Agile Locomotion on Stepping Stones

Adithya Kumar Chinnakkonda Ravi1∗, Victor Dhédin1∗, Armand Jordana2, Huaijiang
Zhu2, Avadesh Meduri2,
Ludovic Righetti2, Bernhard Schölkopf3, Majid Khadiv1 ∗The authors contributed
equally.1Munich Institute of Robotics and Machine Intelligence, Technical
University of Munich, Germany<EMAIL_ADDRESS>{victor.dhedin,
<EMAIL_ADDRESS>School of Engineering, New York University, USA
{hzhu, am9789<EMAIL_ADDRESS>Institute for Intelligent
Systems, Tübingen, Germany<EMAIL_ADDRESS>work was supported by the
Max-Planck Institute for Intelligent Systems, New York University and the
National Science Foundation (CMMI-1825993).

###### Abstract

Legged robots have become capable of performing highly dynamic maneuvers in
the past few years. However, agile locomotion in highly constrained
environments such as stepping stones is still a challenge. In this paper, we
propose a combination of model-based control, search, and learning to design
efficient control policies for agile locomotion on stepping stones. In our
framework, we use nonlinear model predictive control (NMPC) to generate whole-
body motions for a given contact plan. To efficiently search for an optimal
contact plan, we propose to use Monte Carlo tree search (MCTS). While the
combination of MCTS and NMPC can quickly find a feasible plan for a given
environment (a few seconds), it is not yet suitable to be used as a reactive
policy. Hence, we generate a dataset for optimal goal-conditioned policy for a
given scene and learn it through supervised learning. In particular, we
leverage the power of diffusion models in handling multi-modality in the
dataset. We test our proposed framework on a scenario where our quadruped
robot Solo12 successfully jumps to different goals in a highly constrained
environment (video).

##  I INTRODUCTION

Controlling legged robots is a challenging problem, in particular, due to the
need to decide over both continuous (e.g., contact forces) and discrete (e.g.,
which surface patch to step onto next) decision variables. While there has
been much progress in designing controllers based on model predictive control
(MPC) to tackle the continuous part efficiently [1, 2, 3, 4], solving the
mixed problem is still computationally intractable. The main goal of this
paper is to propose an efficient framework based on a combination of nonlinear
MPC (NMPC), Monte Carlo tree search (MCTS), and supervised learning to find
feedback policies that decide over discrete variables of the gait.

Several works have tried to solve the generation of whole-body motion through
a holistic approach, e.g., using differential dynamic programming (DDP) [5],
contact-invariant optimization (CIO) [6], contact-implicit trajectory
optimization [7], and phase-based gait parameterization [8]. While these
approaches have shown impressive behaviors in simulation for legged robots,
they are not suitable for online motion re-planning, mainly due to their large
computation load and being sensitive to initialization.

To enable real-time motion re-planning, most state-of-the-art frameworks
decompose the problem into contact planning and trajectory generation. In this
setting, the contact planner decides which end-effector goes to which contact
patch, i.e., the discrete part of the decisions to be made. The search over
these discrete variables is traditionally done through search-based algorithms
[9], or mixed-integer quadratic programming (MIQP) [10]. While it has been
shown that a relaxed version of the MIQP with L1 norm minimization [11] can be
run on a real quadruped robot in real-time [12], these approaches do not scale
well when the number of discrete variables grows. Furthermore, these
approaches do not take into account the robot dynamics. This makes them
impractical for automatically generating dynamic locomotion behaviours with
flight phases which is the main focus of this paper.

In [13, 14], the authors have included a simplified convex version of the
centroidal momentum dynamics in the MIQP formulation. While an efficient
implementation of the approach can become real-time capable for a small set of
discrete decision variables, this approach can quickly explode as the number
of available contact patches increases. To speed up the search problem, [15]
proposes to learn the outcome of contact transitions using centroidal momentum
dynamics. While this is a valid approach for fast contact planning, it was
never implemented in a physical simulation or real robot.

As an alternative to model-based planning and control, deep reinforcement
learning (DRL) has shown impressive results for agile quadrupedal locomotion
[16, 17, 18]. However, none of these works have shown locomotion in highly
constrained environments like stepping stones. Recently, [19] demonstrated
locomotion on risky terrains like stepping stones. Nevertheless, as mentioned
by the authors, even with the design of several non-trivial reward terms they
failed to learn these motions from scratch, mainly due to the sparsity of the
environment. Hence, they needed to first train a generalist policy that is
able to walk on a simple stepping stone and then fine-tune the policy for each
new environment. In contrast, as we demonstrate in this paper, by exploiting
the structure of the problem and the environment, our approach requires
neither heavy reward shaping nor multi-stage training.

In this paper, we aim to tackle the problem of contact planning for dynamic
maneuvers of legged robots on stepping stones. To do so, we propose to use
Monte Carlo Tree Search (MCTS) which has recently been shown to scale well and
outperform MIQP for dexterous manipulation [20] and gait discovery in
locomotion [21]. In particular, and unlike previous approaches, we use
nonlinear model predictive control (NMPC) inside the MCTS formulation to
ensure the feasibility of the motion in a closed-loop sense. Compared to MIQP
approaches that resort to a simplified model of the system dynamics [13, 22],
here we do not have such simplification and roll out the whole-body NMPC to
evaluate a contact plan. This way, we can make sure that the contact plan is
consistent with the limitations and constraints of the whole control pipeline.
This is in contrast with [23] which learns the contact planner in isolation of
the low-level policy that realizes it. We show that through learning the
feasible solutions of the MCTS, we can have reactive policies that can select
the next contact patches as a function of the current state of the robot and
the environment. Through extensive simulation experiments of quadruped jumping
on stepping stones, we show the effectiveness of our proposed framework.

The main contributions of the paper are as follows:

  * •

We present a novel framework that takes the structure of the dynamics of the
robot and the environment to efficiently control legged robots in highly
limited environments (stepping stones),

  * •

Through a systematic comparison, we show that diffusion models are an ideal
candidate to learn from expert data when multi-modality is present in the
dataset,

  * •

We demonstrate automatic online surface selection for dynamic quadrupedal
locomotion through a learned feedback policy.

Figure 1: A schematic block diagram of the framework. We use the NMPC
formulation in [3].

##  II Contact Planning using MCTS

We formalize our contact planning problem as a Markov-decision process (MDP),
where the state $\bm{s}$ represents the location of end-effectors in contact
with stepping stones which we denote by a discrete index, and the action
$\bm{a}$ selects the next locations for each end-effectors and brings the
system to a new state $\bm{s}^{\prime}=f(\bm{s},\bm{a})$. Each state is
evaluated by a reward function $r(\bm{s})$ that specifies its associated
immediate reward.

To solve this problem, MCTS creates a search tree
$\mathcal{T}=(\mathcal{V}.\mathcal{E})$ where the set of nodes $\mathcal{V}$
contains the visited states and the set of edges contains the visited
transitions $(\bm{s}\overset{\bm{a}}{\rightarrow}\bm{s}^{\prime})$. Each
transition maintains the state-action value $Q(\bm{s},\bm{a})$ and the number
of visits $N(\bm{s},\bm{a})$. MCTS grows this search tree iteratively by the
following steps (see Algorithm 1 for a summary).

  1. 1.

Selection: Start from the root node (initial state) and select successively a
child until a leaf (node that has not been expanded yet or terminal state) has
been reached. When all the children of a node has been expanded, a child is
selected according to a heuristic score (1) that balances exploration and
exploitation during the search.

  2. 2.

Expansion: Unless the selected state from the previous step is a terminal
state, its successor states are added to the tree by enumerating all possible
actions. The corresponding state-action pairs are initialized with
$Q(\bm{s},\bm{a})=0$ and $N(\bm{s},\bm{a})=0$.

  3. 3.

Simulation: From one of the successor states, random actions are performed
onward for a predefined number of steps to create a simulation rollout. The
reward $r$ is evaluated at the end of the simulation.

  4. 4.

Back-propagation: The reward is then back-propagated to update the state-
action value $Q(\bm{s},\bm{a})=Q(\bm{s},\bm{a})+r$ and the number of visits
$N(\bm{s},\bm{a})=N(\bm{s},\bm{a})+1$ for all the states along the node
selected in the selection and expansion steps.

Figure 2: Top view of the environment. The stepping stones in green represent
the goal contact locations, and the stepping stones in red represent the start
contact locations. Some stepping stones are randomly removed.

As shown schematically in Fig. 2, in our problem formulation, the MCTS is
given the initial and final desired locations of all the robot end-effectors
and is asked to compute the feasible set of patches that results in a
successful motion. To balance between the exploration of unvisited and visited
states, we consider the upper confidence bound (UCB) [24] for each node. In
the Selection phase, MCTS selects the action with the highest UCB score
$\bm{a}=\operatorname*{arg\,max}_{\bm{a}}U(\bm{s},\bm{a})$ where the UCB score
is defined as

$\displaystyle
U(\bm{s},\bm{a})=\frac{Q(\bm{s},\bm{a})}{N(\bm{s},\bm{a})}+c\sqrt{\frac{\log
N(\bm{s})}{N(\bm{s},\bm{a})}}$ (1)

where $c$ is a coefficient to balance exploration against exploitation and
$N(\bm{s})=\sum_{\bm{a}}N(\bm{s},\bm{a})$ is the total number of visits for a
state. The reward function that is used to update the state-action value is

$\displaystyle
r(\bm{s})=W\sigma(\frac{1}{N_{e}}\sum_{j=1}^{n_{e}}(1-\frac{||\mathbf{c}^{j}_{W}-{}^{g}\mathbf{c}^{j}_{W}||_{2}}{d_{max}}))$
(2)

where $\mathbf{c}^{j}_{W}$ is the contact location in world frame of the $j$th
end-effector at state $\bm{s}$, ${}^{g}\mathbf{c}^{j}_{W}$ is the desired goal
location for the $j$th end-effector in world frame, and $N_{e}$ is the total
number of end-effectors. $d_{max}$ is the maximum distance between 2 contact
locations in the map. $\sigma:[0,1]\rightarrow[0,1]$ is a function that shapes
the reward. $W\in\\{-1,1\\}$ is a success indicator, by default set to $1$
when the goal state is not reached.

The following adjustments have been made to improve the efficiency of the
algorithm. We do not consider transitions to the same state
$(\bm{s}\overset{\bm{a}}{\rightarrow}\bm{s})$ to limit the redundancy of the
search. During the Expansion phase, to limit the search space, we additionally
check kinematic feasibility to prune the sequences that are either not
reachable or cause self-collision for the legs of the robot. For the
Simulation phase, actions are chosen according to an exploitation/exploration
strategy. With probability $\alpha\in[0,1]$, the next state among the feasible
states is chosen according to a uniform distribution. With probability
$1-\alpha$, the next state is chosen according to a prior distribution
$\widetilde{r}$ such as
$P_{\widetilde{r}}(\bm{s}=\hat{\bm{s}})=r(\hat{\bm{s}})/R$, with $R$ a
normalization factor. Therefore nodes with higher rewards have more chance to
be chosen. It is important to note that, similar to [25], since we can define
a reward at each state that gives us a proxy of how close this state is to the
end goal, we do not simulate a complete rollout to a terminal state as opposed
to game-play (such as chess or go). Instead, we only take one action at the
chosen state and back-propagate the reward from the resulting next state,
hence a rollout depth of one in the Simulation phase. Moreover, we perform
whole-body NMPC using the found sequence of contact locations only when a
terminal state is reached (all the end-effectors are at the desired goal
surface patches). If this contact sequence leads to dynamically feasible
whole-body motion physical simulation, $W$ is set to $1$, else to $-1$.
Ultimately, we back-propagate the final reward. Overall, the MCTS searches for
promising kinematically feasible contact sequences, and the NMPC evaluates
their dynamic feasibility at the end. This way, we avoid performing whole-body
NMPC for every rollout in the Simulation phase. This significantly reduces
search time as whole-body NMPC for a long sequence can be computationally
intensive.

We use the NMPC formulation in [3] to perform the rollouts in MCTS and
evaluate kinematically feasible contact patch candidates. Our NMPC is based on
a decomposition of the problem into a kinematic and dynamic optimizer. In this
setting, given a contact plan (which end-effector goes to which surface
patch), the dynamic optimizer takes into account the centroidal momentum
dynamics and constructs a finite-horizon optimal control problem to find a
feasible set of contact forces and centroidal trajectories. To track these
momentum trajectories, while penalizing contact constraints, a whole-body
kinematic optimal control problem is solved to output the desired whole-body
accelerations. The contact forces from the dynamics optimizer and whole-body
joint trajectories from the kinematic optimizer are then used within the full
dynamics of the system to compute joint torques. Finally, a joint low-
impedance controller is added to the feed-forward torques to account for model
errors and uncertainties. The kinematic and dynamic optimizers re-compute
trajectories at 20 Hz and the joint torques are computed at 1 kHz and are sent
to the robot.

###### Remark 1

In our NMPC formulation, we keep the contact locations fixed and consider them
to be at the center of the stepping stone. As the stepping stones we
considered in the paper are small, this simplification is not limiting in
terms of the number of feasible motions found. However, we can still optimize
for the location of the feet within the stepping stones patches in the NMPC.
We leave this modification for future work.

###### Remark 2

While we can consider all types of gaits and terrain shapes to search over, to
reduce the search space, here we only consider jumping motion over stepping
stones. Note however that the framework can be deployed to consider even a
general acyclic gait at the cost of extremely larger search space. In that
case, we would need to have different cost weights for different gaits.

##  III Learning Contact Planner

While efficient, our MCTS together with the NMPC framework cannot be run in
real-time. To enable the robot to reactively select the next feasible contact
patches given the current ones, we learn a neural network to imitate the MCTS
policy via supervised learning. While MCTS admits a natural extension of a
learnable value function and action probability prior [26, 27], we decide not
to adopt this methodology despite its success in game-play for two reasons.
First, in contrast to generic game-play, locomotion tasks on different maps
(e.g. varying locations and numbers of stepping stones) are likely to have
different states and action space; an MCTS trained on a specific map does not
generalize to other environment maps. Second, some other sensory inputs that
are not modeled in the MCTS state space may provide additional information
(eg. base velocity) on _if_ and _where_ to make the next contact. Therefore,
we take an imitation learning perspective and treat the MCTS as an algorithmic
demonstrator, whose behavior will be cloned by a neural network policy.

Figure 3: High-level contact planning network policy

We collect data of the dynamically feasible solutions discovered by the MCTS
together with the contact locations of the simulation environment. This
dataset is then used to train a neural network in a supervised fashion. A
schematic structure of the network is shown in Fig. 3. To be independent of
the global position of the robot and environment, we express all positions in
the inputs and outputs to the network with respect to the current base frame
$B$ of the robot. The inputs, $\bm{x}$ to the network are the current 3D
position of the available contact locations $\bm{x}_{\text{contact}}$, denoted
by $\bm{c}_{B}^{i}\in\mathcal{C}$ with $i\in\\{1,...,N_{c}\\}$ and
$\mathcal{C}$ the set of available contact locations ($N_{c}$ is the current
number of available contact locations), $\bm{x}_{\text{state}}$ that includes
some state variables (see Sec. III-A) and $\bm{x}_{\text{goal}}$, the final
desired end-effector locations denoted by $\bm{c}_{B}^{g_{j}}$ with
$g_{j}\in\\{1,...,N_{e}\\}$ ($N_{e}$ is the number of end-effectors). The size
of the input is equal to $3\times(N_{c}+2\times(N_{e}+1))$.

The neural network policy predicts the future $H$ contact locations for each
end-effector, denoted by $\hat{\bm{c}}_{B}^{i_{j,h}}$ for the end-effector
$j\in[1,\ldots,N_{e}]$ and the next step $h\in[1,\ldots,H]$. The size of the
output $\hat{\bm{y}}$ is equal to $3\times H\times N_{e}$. As the NMPC is
quite sensitive to small errors in the contact locations, we project the
contact locations given by the policy to the closest center of the patch. We
refer to the projection function as $p_{\mathcal{C}}$ in the following. We
show in Sec. IV-C that the distance between the contact locations given by the
policy and the projected ones is way lower than the half distance between the
two stepping stones. This means that the choice of contact patch by the policy
is not influenced by this projection.

###  III-A Problem specifications

In our setup, a quadruped robot ($N_{e}=4$) navigates in an environment with
up to $81$ stepping stones. Each stone provides the robot with a cylindrical
patch of radius $6.5\,\text{cm}$ to step onto. The stepping stones form a
regular grid so that the feet lay on $4$ stepping stones in the initial
configuration. From this, some stones are randomly removed as explained in the
next section (III-B). The simulated environment can be seen in Fig. 2.

In our experiments, $H$ encodes the NMPC horizon, which in our problem is two
jumps in the future ($H=2$). The goal is to evaluate the learned policy before
each jump and feed the selected contact locations to the NMPC. As state
variables, we consider the position of the end-effectors in the base frame $B$
and the current base linear and angular velocity.

###  III-B Dataset

To collect a diverse dataset, we run MCTS on pre-selected $N_{\text{goals}}$
goal configurations to make the evaluation of the learned policy more
systematic. Starting from the same start position at the center of an
environment, we randomly remove $N_{\text{removed}}$ stepping stones from the
environment and run MCTS for a fixed maximum number of iterations. To collect
diverse paths towards a goal, we keep up to $N_{\text{paths}}$ different
feasible paths for the same goal and environment. To cover a wider range of
robot states, for each MCTS solution, we perturb the simulation
$N_{\text{rand}}$ times and add feasible solutions to the dataset. The
randomization procedure consists of randomizing the initial state of the robot
(position and velocity) as well as the contact locations inside the selected
patches. The training data $(\bm{x},\bm{y})$ (see Fig. 3) are recorded at each
jump ($\bm{y}$ are the contact locations of the next two jumps). We repeat the
procedure on $N_{\text{env}}$ different environments (set of stepping stones).
This dataset is called multigoal in the rest of the paper.

###  III-C Neural network architectures

Our learning problem structure is a selection procedure as the policy should
ideally return contact locations that are given as input. While this can be
achieved using a projection function, some network architectures are suited to
this task such as the Pointer-Network architecture (III-C1) that we consider
as a candidate. Additionally, our dataset is multi-modal as MCTS provides
different contact sequences for the same start and goal contact locations in a
given environment. It is not possible to represent such multi-modal data
distribution with a conventional uni-modal policy class as the model could
collapse to one of the modes or an average over several modes (see Fig 7).
Therefore, we consider Denoising diffusion probabilistic models (DDPMs) [28]
(cf. III-C2), as another potential candidate, since they are theoretically
grounded to handle multi-modality [29] and is practically verified for some
robotic applications [30] [31]. We also consider multi-layer perceptron (MLP)
architecture as a baseline.

####  III-C1 Pointer-Networks

Pointer-Networks [32] take as input a sequence and output discrete indices,
called _pointers_ , that select elements from the input sequence. In this
case, the projection $p_{\mathcal{C}}$ is not performed as the model directly
outputs from the input set. The architecture is composed of $2$ recurrent
networks. The first one, the encoder, encodes the input sequence and provides
as many hidden states as input elements. The second one, the decoder, relies
on an additive attention mechanism that operates on the past decoder’s hidden
states and all the encoder’s hidden states. At each step, the output of the
decoder is the index of the encoder’s hidden state that has the maximum
attention value with the past decoder’s hidden state. This operation is
repeated for as many times as needed. In our case $8$ times (the next two
contact locations for each four legs). The attention mechanism also provides
hints about the confidence of the network in the selection procedure. However,
we did not use this feature in this study.

To make the model select only the contact locations from the input $\bm{x}$,
$\bm{x}_{\text{contact}}$ is given as the input sequence while
$[\bm{x}_{\text{state}},\bm{x}_{\text{goal}}]$ is embedded and given as the
first hidden state of the encoder. Like so, contact locations can be provided
in any order and a different number which is not the case for instance for an
MLP.

####  III-C2 Diffusion models

DDPMs are generative models that map samples from a latent random distribution
to the data distribution in $T$ steps by successive denoising of the original
noise. For each intermediate step $t\in[1,T]$, one can sample a corrupted
input $\bm{x}_{t}$ by adding noise $\bm{\epsilon}_{t}$ to a sample of the data
$\bm{x}_{0}$. A variance schedule assigns an increasing noise level at each
step $t$ so that $\bm{x}_{T}$ can be seen as a pure random noise (usually from
a Gaussian distribution). Those corrupted samples are used to train the
diffusion model $\bm{\epsilon}_{\bm{\theta}}$, parametrized by $\bm{\theta}$,
to estimate the noise added in a supervised manner. It is done by minimizing
the MSE loss between the actual sampled noise $\bm{\epsilon}_{t}$ and the
estimated one $\bm{\epsilon}_{\theta}(\bm{x}_{0}+\bm{\epsilon}_{t},t)$.
Minimizing the MSE loss leads to the minimization of the variational lower
bound of the KL-divergence between the data distribution and the distribution
of samples drawn from the DDPM [28].

To sample with the trained model, noise is successively removed from a random
sample $\bm{x}_{T}$ in the following way

$\bm{x}_{t-1}=a_{t}\left(\bm{x}_{t}-b_{t}\,\bm{\epsilon}_{\theta}(\bm{x}_{t},t)\right)+\sigma_{t}\bm{z}$

where $\bm{z}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$, and $a_{t}$, $b_{t}$,
$\sigma_{t}$ are computed according to the noise schedule.

U-Net-based architectures are widely used by DDPMs, especially for conditional
image generation. Following [30], we consider a conditional U-Net1D with 1D
convolutions applied on the input sequence length (end-effector dimension).
The conditioning is done on both the denoising iteration $t$ of the diffusion
process and the current input $\bm{x}$. Similarly to what has been done with
the Pointer-Network in Sec. III-C1, the input $\bm{x}$ is split into two. We
apply a multi-head attention layer [33] with
$[\bm{x}_{\text{state}},\bm{x}_{\text{goal}}]$ as query and
$\bm{x}_{\text{contact}}$ as key and values. Thus, contact locations can be
shuffled and provided in any number. The resulting embedding of
$\bm{x}_{\text{state}}$ is concatenated to a sinusoidal position embedding
[33] of $t$ to form the input of a feature-wise linear modulation (FiLM) [34]
layer, as proposed by [31]. The architecture is detailed in Fig. 4. Finally,
the conditioning vector is added to all layers of the encoder and the decoder
parts of the network, as done in [30].

Figure 4: High-level description of the conditional U-Net1D with multi-head
attention. sin refers to Sinusoidal positional embedding applied to $t$.

##  IV Results and Discussion

In this section, we present the results of applying the proposed framework to
control the Solo12 quadruped robot [35] on stepping stones in the Pybullet
simulation environment [36]. We start by showing the capability of the MCTS
together with NMPC to generate non-trivial and physically realizable contact
plans on stepping stones (IV-A). Then, through a simple example scenario, we
compare the capability of different neural network architectures to handle
multi-modality in the dataset (IV-B). Then, we present the result of using the
diffusion model to learn a policy that can output contact plans from the state
of the robot and the environment (IV-C). Finally, we show the effectiveness of
our learned policy in selecting the feasible contact patches as the
environment and goal change dynamically.

###  IV-A MCTS contact planner

To test the effectiveness of the MCTS with NMPC in generating feasible jumping
motions on stepping stones, $N_{\text{removed}}=12$ stepping stones are
randomly removed from the map. For each episode, there are
$N_{\text{box}}=81-12=69$ stepping stones in the environment. Given that for
jumping motion all $N_{e}$ end-effectors are in contact, the number of
configurations for each jump is equal to $N_{\text{box}}^{N_{e}}$. This rules
out the use of a brute-force search to solve the problem. Pruning drastically
reduces the number of possible configurations, decreasing our branching factor
to 324 at maximum. However, this is still high.

As hyperparameters, we found out that $c=0.03$ and $\sigma=sigmoid(T(x-1))$
with $T=5$ leads to a sufficient balance between exploration and exploitation.
$\sigma$ shapes the reward so that it is significantly increased for a state
close to the goal state leading to a better exploitation for nodes close to
the terminal node.

We evaluate our method over 100 environments by estimating the execution time,
the number of iterations, and the number of NMPC rollouts performed to find
the first feasible solutions with a maximum number of iterations equal to
$10^{4}$. MCTS was able to find 58 feasible contact plans. For those, the
median execution time is 8.35 sec, the median number of iterations is 930 and
the average number of NMPC simulations performed is 5.8. The execution time
depends mostly on the contact sequence length and the number of NMPC
simulations as shown in Fig. 5. It is likely that some of the environments
admitted no feasible solution. Indeed, by doing the same experiment without
removing any contact locations, MCTS found 94 solutions out of 100 under
$10^{4}$ iterations.

Qualitatively, contact plans are quite similar to each other as MCTS explores
nodes close to the goal more often, as they have been visited less. Multi-
modality is more pronounced at short horizons, as can be seen in Fig. 6. Some
interesting and non-trivial motions found by the MCTS are also shown in the
accompanying video.

Figure 5: performance analysis of the MCTS. The median execution time is 8.35s
(blue line on the color bar). After 4000 iterations, it is unlikely to find a
feasible contact plan (here only 4 out of 100). NMPC simulation is what takes
most of the time during the search.

Figure 6: 3 contact sequences generated by MCTS on the same environment. Each
path following the arrows is a contact sequence for an individual foot.
Starting contact locations are circled in red and goal contact locations in
green (the contact location at the center circled in green is both a starting
and goal contact location).

###  IV-B Choice of network architecture

To systematically compare the performance of different neural network
architectures, we design a simple learning problem. In this problem, the
dataset, multimod, consists of 2 different contact plans to reach the same
goal (in diagonal) starting from the same position in the same environment.
One plan goes first straight and left, the other one goes left and then
straight to reach the goal. We record 50 randomized runs for each plan for a
total of 700 samples in the dataset.

Since we perturb the states while collecting data, the dataset is not strictly
multi-modal. We would like to evaluate how this partial multi-modality affects
the trained policy to reproduce the variety of contact plans of the dataset.
Generating diverse contact plans is beneficial as some might be more relevant
in specific states.

We compare 3 different model architectures presented in III-C on this dataset.
The training parameters are detailed in the appendix.

####  IV-B1 Results

As shown in Fig. 7, only the diffusion model was able to reproduce the two
possible modes of the solution. On 20 simulations, the diffusion model reached
the goal by going 12 times first to the right and 8 times first up (we used a
different seed each time to generate the initial noise). MLP and Pointer-
Network reached the goal respectively 20 and 19 times (Pointer-Network failed
once because of the randomization procedure, all the predicted plans were the
same) but collapsed to one of the solutions, as can be seen on the base
trajectories in Fig. 7. As we would like to benefit from the multiple feasible
solutions of MCTS for a given goal, we focus on diffusion models for the
navigation task in the following section. Note that, for the diffusion model,
the projection error
$||\hat{\bm{c}}_{B}^{i_{j,h}}-p_{\mathcal{C}}(\hat{\bm{c}}_{B}^{i_{j,h}})||_{2}$
is 5.4mm on average and at a maximum 2.2cm for the next first jump, which is
approximately 3 times less than the half distance between two stepping stones.

(a) MLP

(b) PtrNet

(c) Diffusion model

Figure 7: Base position $(x,y)$ for 20 randomized runs with models trained on
the multimod dataset. Starting contact locations are circled in red and goal
contact locations in green. The base position is recorded every 20ms. The size
of the circles represents the height of the base. Only the diffusion model
7(c) is able to reproduce both two contact sequences of the training dataset.

###  IV-C Learning to navigate on stepping stones

In this subsection, we aim to qualitatively show that the learned policy can
reactively generate feasible contact plans when the environment or the goals
change on the fly. To do that, we generated $N_{\text{env}}=80$ different
environments with $N_{\text{removed}}=12$ and $N_{\text{goals}}=8$ (up, down,
left, right and diagonals) at least 2 jumps away from the starting location,
$N_{\text{paths}}=3$ and $N_{\text{repeat}}=5$. Our training dataset,
multigoal, has 52224 samples in total. Our test environment has been generated
in the same conditions on $N_{\text{env}}=20$ different environment. We tested
our policy on test goals for which MCTS found a feasible contact plan. Our
training procedure and hyperparameters are similar to [31] and are detailed in
the appendix.

As done in [31], we used a DDIM [37] approach to decouple the number of
denoising iterations in training and inference. Denoising the same number of
training iterations takes more time but usually leads to better results. For
15 denoising iterations and above, the policy achieved approximately a similar
success rate with a policy trained on $T=50$ training iterations. Therefore,
we used 15 denoising iterations in the experiments. It takes 70 ms for the
policy to be evaluated on an AMD Ryzen 5 5625u CPU. We expect a significant
decrease in the computation time on a GPU.

####  IV-C1 Static environments

Our policy has been trained to output different contact locations for the next
jump as our dataset contains multi-modal samples. Therefore, by combination,
the final full-length contact plan is likely to differ from the ones of MCTS.
This is confirmed as $40\%$ of the successful contact plans were not in the
training dataset when replaying on the training environment. However, since 3
different paths maximum are recorded for each goal in a given environment, it
could lead to a poorer success rate as the policy might end up out of the
training distribution. When replayed in the same conditions 20 times for a
goal in a corner of the environment, our trained policy reached the goal all
the times in 14 different ways. Some examples can be seen in the accompanying
video. This confirms the results of section Sec. IV-B with a policy trained on
a more diverse dataset.

####  IV-C2 Dynamic environments

Now, we proceed to evaluate how well the policy can be used in a dynamic
environment. To do that, we randomly remove one contact location from the
environment before each jump while the robot is reaching a goal. The removed
stepping stones are chosen among the ones found by MCTS. This way, it is more
likely to remove a stepping stone that would have been initially selected,
which shows the ability of the policy to replan reactively. This is a
challenging task. Removing stepping stones could make the robot jump into a
position that is bound to fail or that is out of the training data
distribution. Our policy sporadically succeeded on this task ($22\%$ of the
time on 50 trials on a goal in diagonal). Additionally, our trained policy was
also able to perform navigation to reach some user-defined changing goals in a
new environment with 12 removed stepping stones. Successful examples of those
two tasks can be seen in the accompanying video.

##  V Conclusions and future work

In this paper, we presented a framework to efficiently search for contact
plans in a known environment for legged robots. In our framework, we proposed
a customized version of MCTS together with an NMPC to search for feasible
solutions on stepping stones. We showed that we can reliably find feasible
solutions for different arrangements of the stepping stones as well as initial
and final conditions. Collecting feasible rollouts enabled us to collect a
rich dataset and learn a control policy that can generate contact plans for
the NMPC reactively. Our extensive analyses showed that, due to the multi-
modal nature of the dataset, using diffusion models is an ideal way to perform
supervised learning on the dataset.

In the future, we plan to explore alternative architectures for diffusion
models, as the one used in this study is derived from works with different
learning problem structures. Also, we would like to implement the learned
high-level contact planner together with the NMPC on the real hardware. To
further reduce the computation time of the whole control pipeline for real-
world experiments, we are also interested in replacing the NMPC with a learned
low-level policy [38]. We also aim to extend our framework to optimize the
gait as well as the surface patches, given a local environment map. Finally,
we are interested in learning the high-level policy from partial observation
of the robot on-board cameras.

TABLE I: Training parameters for the models. EMA: Exponential Moving Average. wd: weight decay. Model | Dataset | Optimizer | LR | LR scheduler | Batch Size | Epochs | Additional Parameters
---|---|---|---|---|---|---|---
MLP | multimod | Adam | $10^{-3}$ | Exp. decrease (0.996) | 32 | 500 |  for removed contacts height is set to 0
Pointer-Network | multimod | Adam | $10^{-3}$ | Exp. decrease (0.996) | 8 | 500 | shuffled contact
Diffusion Model | multimod |  AdamW (EMA: 0.75 + wd: $10^{-6}$) | $10^{-3}$ |  Cosine annealing w. linear warmup | 32 | 500 | shuffled contact
Diffusion Model | multigoal |  AdamW (EMA: 0.75 + wd: $10^{-6}$) | $10^{-3}$ |  Cosine annealing w. linear warmup | 64 | 2000 | shuffled contact

## References

  * [1] J.-P. Sleiman, F. Farshidian, M. V. Minniti, and M. Hutter, “A unified mpc framework for whole-body dynamic locomotion and manipulation,” IEEE Robotics and Automation Letters, vol. 6, no. 3, pp. 4688–4695, 2021.
  * [2] E. Daneshmand, M. Khadiv, F. Grimminger, and L. Righetti, “Variable horizon mpc with swing foot dynamics for bipedal walking control,” IEEE Robotics and Automation Letters, vol. 6, no. 2, pp. 2349–2356, 2021.
  * [3] A. Meduri, P. Shah, J. Viereck, M. Khadiv, I. Havoutis, and L. Righetti, “Biconmp: A nonlinear model predictive control framework for whole body motion planning,” arXiv preprint arXiv:2201.07601, 2022.
  * [4] C. Mastalli, W. Merkt, G. Xin, J. Shim, M. Mistry, I. Havoutis, and S. Vijayakumar, “Agile maneuvers in legged robots: a predictive control approach,” arXiv preprint arXiv:2203.07554, 2022.
  * [5] Y. Tassa, T. Erez, and E. Todorov, “Synthesis and stabilization of complex behaviors through online trajectory optimization,” in 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 4906–4913, IEEE, 2012.
  * [6] I. Mordatch, E. Todorov, and Z. Popović, “Discovery of complex behaviors through contact-invariant optimization,” ACM Transactions on Graphics (TOG), vol. 31, no. 4, pp. 1–8, 2012.
  * [7] M. Posa, C. Cantu, and R. Tedrake, “A direct method for trajectory optimization of rigid bodies through contact,” The International Journal of Robotics Research, vol. 33, no. 1, pp. 69–81, 2014.
  * [8] A. W. Winkler, C. D. Bellicoso, M. Hutter, and J. Buchli, “Gait and trajectory optimization for legged systems through phase-based end-effector parameterization,” IEEE Robotics and Automation Letters, vol. 3, no. 3, pp. 1560–1567, 2018.
  * [9] S. Tonneau, A. Del Prete, J. Pettré, C. Park, D. Manocha, and N. Mansard, “An efficient acyclic contact planner for multiped robots,” IEEE Transactions on Robotics, vol. 34, no. 3, pp. 586–601, 2018.
  * [10] R. Deits and R. Tedrake, “Footstep planning on uneven terrain with mixed-integer convex optimization,” in 2014 IEEE-RAS international conference on humanoid robots, pp. 279–286, IEEE, 2014.
  * [11] S. Tonneau, D. Song, P. Fernbach, N. Mansard, M. Taïx, and A. Del Prete, “Sl1m: Sparse l1-norm minimization for contact planning on uneven terrain,” in 2020 IEEE International Conference on Robotics and Automation (ICRA), pp. 6604–6610, IEEE, 2020.
  * [12] F. Risbourg, T. Corbères, P.-A. Léziart, T. Flayols, N. Mansard, and S. Tonneau, “Real time footstep planning and control of the solo quadruped robot in 3d environments,” 2022.
  * [13] B. Ponton, A. Herzog, S. Schaal, and L. Righetti, “A convex model of humanoid momentum dynamics for multi-contact motion generation,” in 2016 IEEE-RAS 16th International Conference on Humanoid Robots (Humanoids), pp. 842–849, IEEE, 2016.
  * [14] B. Aceituno-Cabezas, C. Mastalli, H. Dai, M. Focchi, A. Radulescu, D. G. Caldwell, J. Cappelletto, J. C. Grieco, G. Fernández-López, and C. Semini, “Simultaneous contact, gait, and motion planning for robust multilegged locomotion via mixed-integer convex optimization,” IEEE Robotics and Automation Letters, vol. 3, no. 3, pp. 2531–2538, 2017.
  * [15] Y.-C. Lin, B. Ponton, L. Righetti, and D. Berenson, “Efficient humanoid contact planning using learned centroidal dynamics prediction,” in 2019 International Conference on Robotics and Automation (ICRA), pp. 5280–5286, IEEE, 2019.
  * [16] M. Bogdanovic, M. Khadiv, and L. Righetti, “Model-free reinforcement learning for robust locomotion using demonstrations from trajectory optimization,” Frontiers in Robotics and AI, vol. 9, 2022.
  * [17] Z. Zhuang, Z. Fu, J. Wang, C. Atkeson, S. Schwertfeger, C. Finn, and H. Zhao, “Robot parkour learning,” arXiv preprint arXiv:2309.05665, 2023.
  * [18] X. Cheng, K. Shi, A. Agarwal, and D. Pathak, “Extreme parkour with legged robots,” arXiv preprint arXiv:2309.14341, 2023.
  * [19] C. Zhang, N. Rudin, D. Hoeller, and M. Hutter, “Learning agile locomotion on risky terrains,” arXiv preprint arXiv:2311.10484, 2023.
  * [20] H. Zhu, A. Meduri, and L. Righetti, “Efficient object manipulation planning with monte carlo tree search,” in 2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 10628–10635, IEEE, 2023\.
  * [21] L. Amatucci, J.-H. Kim, J. Hwangbo, and H.-W. Park, “Monte carlo tree search gait planner for non-gaited legged system control,” in 2022 International Conference on Robotics and Automation (ICRA), pp. 4701–4707, IEEE, 2022.
  * [22] J. Wang, S. Kim, S. Vijayakumar, and S. Tonneau, “Multi-fidelity receding horizon planning for multi-contact locomotion,” in 2020 IEEE-RAS 20th International Conference on Humanoid Robots (Humanoids), pp. 53–60, IEEE, 2021\.
  * [23] V. Tsounis, M. Alge, J. Lee, F. Farshidian, and M. Hutter, “Deepgait: Planning and control of quadrupedal gaits using deep reinforcement learning,” IEEE Robotics and Automation Letters, vol. 5, no. 2, pp. 3699–3706, 2020.
  * [24] L. Kocsis and C. Szepesvári, “Bandit based monte-carlo planning,” in European conference on machine learning, pp. 282–293, Springer, 2006.
  * [25] Y. Labbé, S. Zagoruyko, I. Kalevatykh, I. Laptev, J. Carpentier, M. Aubry, and J. Sivic, “Monte-carlo tree search for efficient visually guided rearrangement planning,” IEEE Robotics and Automation Letters, vol. 5, no. 2, pp. 3715–3722, 2020.
  * [26] D. Silver et al., “Mastering the game of go without human knowledge,” Nature, vol. 550, no. 7676, pp. 354–359, 2017.
  * [27] J. Schrittwieser, I. Antonoglou, T. Hubert, K. Simonyan, L. Sifre, S. Schmitt, A. Guez, E. Lockhart, D. Hassabis, T. Graepel, et al., “Mastering atari, go, chess and shogi by planning with a learned model,” Nature, vol. 588, no. 7839, pp. 604–609, 2020.
  * [28] J. Ho, A. Jain, and P. Abbeel, “Denoising diffusion probabilistic models,” 2020\.
  * [29] A. Block, A. Jadbabaie, D. Pfrommer, M. Simchowitz, and R. Tedrake, “Provable guarantees for generative behavior cloning: Bridging low-level stability and high-level behavior,” 2023.
  * [30] M. Janner, Y. Du, J. B. Tenenbaum, and S. Levine, “Planning with diffusion for flexible behavior synthesis,” 2022.
  * [31] C. Chi, S. Feng, Y. Du, Z. Xu, E. Cousineau, B. Burchfiel, and S. Song, “Diffusion policy: Visuomotor policy learning via action diffusion,” 2023.
  * [32] O. Vinyals, M. Fortunato, and N. Jaitly, “Pointer networks,” 2017.
  * [33] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, “Attention is all you need,” 2023.
  * [34] E. Perez, F. Strub, H. de Vries, V. Dumoulin, and A. Courville, “Film: Visual reasoning with a general conditioning layer,” 2017.
  * [35] F. Grimminger, A. Meduri, M. Khadiv, J. Viereck, M. Wüthrich, M. Naveau, V. Berenz, S. Heim, F. Widmaier, T. Flayols, et al., “An open torque-controlled modular robot architecture for legged locomotion research,” IEEE Robotics and Automation Letters, vol. 5, no. 2, pp. 3650–3657, 2020.
  * [36] E. Coumans and Y. Bai, “Pybullet, a python module for physics simulation for games, robotics and machine learning.” http://pybullet.org, 2016–2021.
  * [37] J. Song, C. Meng, and S. Ermon, “Denoising diffusion implicit models,” 2022.
  * [38] M. Khadiv, A. Meduri, H. Zhu, L. Righetti, and B. Schölkopf, “Learning locomotion skills from mpc in sensor space,” in Learning for Dynamics and Control Conference, pp. 1218–1230, PMLR, 2023.

## Appendix

In this appendix, we present the detailed values used for training the
netweorks. Training parameters can be seen on Tab. I.

####  V-1 MLP

We used a standard MLP with 4 hidden layers of latent dimension 64 with
LeakyReLU activation. The size input dimension is 273 (81 stepping stones
locations) and the size output dimension is 24. The total number of trainable
parameters of this model is 35736.

####  V-2 Pointer-Networks

Both the encoder and decoder are LSTMs with 2 layers and hidden dimension 32.
The attention mechanism is also of hidden dimension 32.
$[\bm{x}_{\text{state}},\bm{x}_{\text{goal}}]$ is embedded through an MLP with
2 hidden layers of hidden dimensions 16 and 32 with PReLU activation. The
total number of trainable parameters is 31521 which is comparable to the MLP
architecture considered.

When training on the large dataset with removed stepping stones we chose a
hidden dimension of 256 and 3 lstm layers.

####  V-3 U-Net1D

We chose a U-Net with 3 layers with respectively a channel width equal to 64,
128 and 256 and a kernel size of 3 for the convolutions. Convolutions are
sliding on the temporal/end-effector dimension. The sinusoidal embedding
dimension is 32. The multi-head attention layer has only one head. We used a
squared cosine noise schedule with $\beta_{1}=0.004$ and $\beta_{T}=0.02$ as
suggested by [28] for $T=50$ training iterations. The model has in total
2639207 parameters.


# Strong surface termination dependence of the electronic structure of polar
superconductor LaFeAsO revealed by nano-ARPES

Sung Won Jung<EMAIL_ADDRESS>Diamond Light Source, Harwell Campus,
Didcot, OX11 0DE, United Kingdom  Department of Physics and Research Institute
of Natural Science, Gyeongsang National University, Jinju 52828, Republic of
Korea  Luke C. Rhodes  School of Physics and Astronomy, University of St.
Andrews, St. Andrews KY16 9SS, United Kingdom  Matthew D. Watson  Diamond
Light Source, Harwell Campus, Didcot, OX11 0DE, United Kingdom  Daniil V.
Evtushinsky  Laboratory for Quantum Magnetism, Institute of Physics, École
Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland  Cephise
Cacho  Diamond Light Source, Harwell Campus, Didcot, OX11 0DE, United Kingdom
Saicharan Aswartham  Leibniz Institute for Solid State and Materials Research,
01171 Dresden, Germany  Rhea Kappenberger  Leibniz Institute for Solid State
and Materials Research, 01171 Dresden, Germany  Sabine Wurmehl  Leibniz
Institute for Solid State and Materials Research, 01171 Dresden, Germany
Bernd Büchner  Leibniz Institute for Solid State and Materials Research, 01171
Dresden, Germany  Institut für Festkörperphysik, Technische Universität
Dresden, D-01171 Dresden, Germany  Timur K. Kim<EMAIL_ADDRESS>Diamond Light Source, Harwell Campus, Didcot, OX11 0DE, United Kingdom

###### Abstract

The electronic structures of the iron-based superconductors have been
intensively studied by using angle-resolved photoemission spectroscopy
(ARPES). A considerable amount of research has been focused on the LaFeAsO
family, showing the highest transition temperatures, where previous ARPES
studies have found much larger Fermi surfaces than bulk theoretical
calculations would predict. The discrepancy has been attributed to the
presence of termination-dependent surface states. Here, using photoemission
spectroscopy with a sub-micron focused beam spot (nano-ARPES) we have
successfully measured the electronic structures of both the LaO and FeAs
terminations in LaFeAsO. Our data reveal very different band dispersions and
core-level spectra for different surface terminations, showing that previous
macro-focus ARPES measurements were incomplete. Our results give direct
evidence for the surface-driven electronic structure reconstruction in
LaFeAsO, including formation of the termination-dependent surface states at
the Fermi level. This new experimental technique, which we have shown to be
very powerful when applied to this prototypical compound, can now be used to
study various materials with different surface terminations.

##  I Introduction

LaFeAsO is the parent compound of the so-called “1111” family of iron-based
superconductors – the family in which superconductivity was first discovered
by Kamihara et al. Kamihara  _et al._ (2008) and in which the highest reported
superconducting transition temperatures have been found Zhi-An  _et al._
(2008); Wang  _et al._ (2008). Despite being among the first compounds
discovered, many questions remain as to its experimental electronic structure.
This is partly because of the lower quality of the original single crystals,
but more fundamentally is due to the fact that its structure includes
alternately charged LaO and FeAs planes. Thus upon cleaving the sample, it can
have two surface terminations, each of which is a polar surface with
uncompensated charge, and unrepresentative of the bulk.

Figure 1: Termination dependent structure of LaFeAsO surface (a) layered
crystal structure with polarised LaO and FeAs surface terminations; (b) sketch
of the cleaved crystal surface with various terminations.

Earlier angular resolved photoemission spectroscopy (ARPES) studies Lu  _et
al._ (2009); Liu  _et al._ (2010a); Yang  _et al._ (2010); Zhang  _et al._
(2016); Huh  _et al._ (2020) revealed an electronic structure that
substantially differed from the expectations from naive bulk density
functional theory (DFT) calculations. In particular, the spectra included a
bright hole-like band around the Brillouin zone center. It was quickly
appreciated that this band, and several other prominent bands, were
attributable to the surface FeAs layer which has an effective hole-doping Liu
_et al._ (2010a), with the overall spectra being somewhat reminiscent of the
hole-doped end-compound KFe2As2 but totally unlike the spectra on other
nominally charge-compensated parent compounds such as BaFe2As2 Fedorov  _et
al._ (2019); Watson  _et al._ (2019) and NaFeAs Watson  _et al._ (2018).
However, the spectra also contain bands that are not attributed to the topmost
FeAs layer, and it has been a significant challenge to establish which, if
any, features can be attributed to the bulk electronic structure of LaFeAsO.

One of the possible strategies to reveal the real bulk electronic structure of
LaFeAsO in ARPES is to deliberately suppress these surface contribution, for
example by dosing alkali metals on freshly cleaved surface of the crystal in
ultra-high vacuum Yang  _et al._ (2010); Zhang  _et al._ (2016). However, this
also leads to smearing of the photoemission spectra and therefore
significantly reduces the quality of obtained data and confidence in the
obtained values of the Fermi surface dimensions. Moreover, while single
surface terminations can be accessed using microscopic techniques such as
scanning tunneling microscopy (STM) Zhou  _et al._ (2011a), the much larger
spot-size of conventional ARPES light sources ($\sim$1 mm for the He-lamp lab-
based and $\sim$100 micron for the synchrotron radiation sources) means that
previous studied have all measured the superposition of the signals from two
different terminations.

In order to overcome these issues and directly probe the electronic structure
of different terminations of LaFeAsO, we used the nano-ARPES technique, in
which one could achieve sub-micron spatial resolution by focusing incoming
light with Fresnel zone plates (FZP). This method has been successfully
applied to other iron-based superconductors to reveal electronic structure of
the single antiferromagnetic domain as in case of BaFe2As2 Watson  _et al._
(2019), or nematic domain as in case of FeSe Rhodes  _et al._ (2020). By
analysing the chemical shift and overall intensity of As $3d$ core levels, we
identify the spatial distribution of surface terminations, and obtain
different valence band spectra from different positions. From this, we can
directly assign the bands observed near the Fermi level as originating from
either a FeAs termination, or the FeAs layer below an LaO termination. Our
results are compared with DFT calculations allowing for relaxation at the
surface Eschrig  _et al._ (2010).

##  II Methods

High quality single crystals of LaFeAsO were synthesized by solid state single
crystal growth method as described elsewhere Rhea Kappenberger and Saicharan
Aswartham and Francesco Scaravaggi and Christian G.F. Blum and Mihai I. Sturza
and Anja U.B. Wolter and Sabine Wurmehl and Bernd Büchner (2018). Nano-ARPES
measurements were performed using 70 eV photon energy in linear horizontal
polarisation at the nano-ARPES branch of I05-ARPES beamline, using a Scienta
DA30 electron analyser. Reference high resolution ARPES data at 70 eV photon
energy in linear horizontal polarisation were obtained at the HR-ARPES branch
of the same I05-ARPES beamline at Diamond Light Source using R4000 electron
energy analyser Hoesch  _et al._ (2017). All measurements were taken at low
temperatures (below 40 K for nano-ARPES and below 10 K for HR-ARPES) and in
ultra-high vacuum (better than $2\times{}10^{-10}$ mbar). Samples were cleaved
in vacuum and oriented along the $\Gamma-M$ direction of the Brillouin Zone.
The raw data for the spatial map of As $3d$ core levels data have been
normalized to take account of the MCP detector sensitivity, using reference
data measured in swept mode. DFT calculations were performed using open-source
package Quantum Espresso Giannozzi  _et al._ (2017) using Projected Augmented
Wavefunctions (PAW) and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation
functional. A kinetic energy cutoff of 80 Ry and 960 Ry were used for the
wavefunction and charge density respectively for both bulk and slab
calculations. A $k$-grid of $8\times{}8\times{}4$ and $6\times{}6\times{}1$
were used for the bulk and slab calculations respectively. For the slab
calculations, 5 layers of FeAs were interspaced with 4 layers of LaO to create
a FeAs terminated surface, and vice versa for the LaO terminated surface,
equivalent to Ref. Eschrig  _et al._ (2010). Then 15 Å of vacuum was added to
both sides, and these structures were then allowed to relax until the total
force was less than 0.001 a.u. For details, please refer to the input files
attached in the supplemental material.

##  III Results

The crystal structure of LaFeAsO consists of FeAs layers, the essential
building blocks of the Fe-based superconductor family, alternated with layers
of LaO, which can be considered as an insulating spacer layer. In order to
prepare a surface for photoemission spectroscopy, the material is cleaved
between the layers, leaving two possible surface terminations Eschrig  _et
al._ (2010). As shown schematically in Fig. 1(a), both of the terminations
result in polar surfaces, with different charge from the bulk value, leading
to unique 2D electronic states at each of these surfaces. ARPES is a surface-
sensitive technique that can access such surface states, however with a
macroscopic beam spot, it is inevitable that the signal from both terminations
is superimposed.

Figure 2: As $3d$ spatial maps of LaFeAsO (a) spatial resolved 2D map of As
$3d$ core levels total intensity; (b) spatially integrated As $3d$ core level
spectra with fit showing contribution from surface (P1), subsurface (P2) and
bulk (P3) FeAs layers; (c,d,e) spatial resolved 2D maps of As $3d$ core levels
fit contributions from surface, subsurface and bulk FeAs layers; (f,g) area
ratio of the fit components for FeAs and LaO terminations.

To establish the termination-dependent electronic structures, we used a small
beam spot ($\sim$800 nm), and scanned the sample while measuring the
photoelectrons to map out the sample both spatially and electronically. In
order to quantitatively map the spatial distribution of terminations, we
measured As $3d$ core-level spectra across a region of the sample. The total
intensity, shown in Fig.2(a), is fairly uniform except for a region of lower
intensity (red color). On the basis of the intensity alone, we could already
assign this as a region of LaO termination, since the As atoms in this case
are in the subsurface layer and therefore should be weaker in intensity
compared to FeAs terminations. However the chemical shift analysis of the As
$3d$ core states in the different layers allows a more thorough fingerprinting
of these different terminations.

The summation of the As $3d$ spectra across the whole spatial map is presented
in Fig. 2(b), clearly showing contributions from several peaks, identified in
fitting analysis as doublets P1, P2 and P3. While the spectrum is complex, we
can identify the doublet at lowest binding energy, P1, as arising from the
topmost As on an FeAs termination, since it has the largest chemical shift.
This assignment is further supported by our additional photon-energy and
geometry-dependent measurements performed with a larger beam spot on mixed-
termination samples (see Supplementary Figure A1) which show that the P1 peak
is most prominent in conditions which are most surface-sensitive, namely at
lower photon energy and at grazing incidence. Our assignment of P1 is also
consistent with previous experiments where the intensity of low binding energy
peak drastically decreased with alkali metal dosing of the surface Yang  _et
al._ (2010); Zhang  _et al._ (2016). The doublet with highest binding energy,
P3, shows the opposite dependence to P1, and thus corresponds to bulk-like As
layers from deeper in the sample, where the surface charge is fully screened.
The doublet with intermediate binding energy, P2, can be attributed to the As
atoms that are neither bulk-like nor on the topmost surface; either from the
As atoms below the Fe within the FeAs surface termination, or from the top As
in the sub-surface FeAs layer below an LaO termination.

Using such three doublet fit of the summed data to constrain parameters, we
then fit the spectra on a pixel-by-pixel basis, allowing the amplitude of the
peaks to vary, leading to the spatial maps displaying the dependence of the
intensity of the P1, P2 and P3 peaks in Fig. 2(c-e), from which we also take
two line cuts in (f,g). For cut 1 in Fig. 2(f), corresponding to a uniform
region of FeAs termination, all the peaks P1, P2 and P3 contribute to the
spectra, with very small variation. However, for cut 2 in Fig. 2(g), we find a
region approximately 10 microns wide, in which the P2 peak completely
dominates the spectral intensity, while P1 and P3 are suppressed. The
suppression of the P1 peak, which was assigned to the FeAs termination,
implies that the dominant P2 peak here comes from the LaO termination (i.e.
from the top As in the sub-surface FeAs layer). Thus, our spatial mapping
analysis not only identifies the different contributions to the As $3d$ core
level spectra, but also allows us to confidently place the beam either on a
region with FeAs termination, or to measure on the areas with LaO termination.

Figure 3: Termination dependent electronic structure of LaFeAsO (a,b,c) band
dispersions at the Fermi level from large beam spot at high-resolution branch
and FeAs and LaO terminations submicron beam spot at nano-ARPES branch;
(d,e,f) valence band dispersions from large beam spot at high-resolution
branch and FeAs and LaO terminations submicron beam spot at nano-ARPES branch;
(g,h,i) As $3d$ core level spectra from large beam spot at high-resolution
branch and FeAs and LaO terminations submicron beam spot at nano-ARPES branch.

With the knowledge of the distribution of the terminations established, we
turn to the corresponding electronic structure of LaFeAsO. First, in Fig.
3(a), we show the measured band dispersions at the Fermi level and valence
bands using high-resolution ARPES, but with a beam spot of $\sim$80 $\mu$m.
From the conventional ARPES data of LaFeAsO in Fig. 3(a), we label a number of
the most easily identifiable bands as $\alpha$,$\beta$,$\gamma$,$\delta$;
these have all been reported in previous studies Liu  _et al._ (2010a); Yang
_et al._ (2010); Zhang  _et al._ (2016). While there have been several
attempts to identify the origin of all the bands from measurements with a
macroscopic beam spot Zhang  _et al._ (2016) and comparison with DFT
calculations Eschrig  _et al._ (2010), our methodology relies on more direct
approach to measure electronic structure of various surface terminations.

In Fig. 3(b), taken on a FeAs termination, we observe most prominently the
$\alpha$ band, which confirms it as a surface FeAs band Liu  _et al._ (2010a);
Yang  _et al._ (2010); Zhang  _et al._ (2016), which is the brightest feature
in most measurement conditions. The kink-like structure of the $\alpha$ band
has been previously discussed Liu  _et al._ (2010b); Yang  _et al._ (2010) and
is evidence of strong electron-electron interactions in the hole-doped surface
layer. In addition, the $\beta$ band contributes to the spectrum on the FeAs
termination. On the other hand, these features are all suppressed when the
spot is moved onto the LaO termination, shown in Fig. 3(c), and instead the
spectrum includes an electron-like $\gamma$ band at the M point. This
electron-like band, and also the hole-like $\delta$ band, originate from Fe-
derived bands from the subsurface FeAs layer.

The electronic structure that emerges from the FeAs subsurface layer below the
LaO termination has previously been obscured in measurements on mixed
termination samples, since ARPES is highly surface sensitive and the bands
from the FeAs termination dominate the intensity. However, it is noteworthy
for several reasons. First, the $\gamma$ electron band - most likely
comprising two electron-like bands in close proximity, unresolved in the
current measurements - is reasonably large, with a $k_{F}\approx{}$ 0.25 Å-1
and $E_{B}\approx{}$-80 meV and the observed hole-like states at $\Gamma$ in
Fig. 3(c) are fully occupied, implying an uncompensated system with net
electron doping. The subsurface layer was predicted to have a net charge of
0.114 electrons per Fe compared with the bulk material Eschrig  _et al._
(2010), qualitatively consistent with our experiments. Conceptually, the LaO
layer does a relatively poor job of screening the polar surface charge, such
that there is still a substantial additional charge in the subsurface layer
compared with a bulk FeAs layer. Second, we do not observe any electronic
reconstruction, which implies that the effective electron doping of the
subsurface FeAs layer is sufficient to move beyond the antiferromagnetic
phase, which in the bulk survives to $x\approx{}$0.05 Scaravaggi  _et al._
(2021). Third, we note an fascinatingly close analogy between the electronic
structure of the subsurface FeAs and the high-$T_{c}$ superconducting
monolayer FeSe/STO: both are strictly 2D systems, involving an interface with
an oxide material, and in which a net electron doping of $\sim$0.1 electrons
per Fe occurs to yield a Fermi surface containing only electron pockets Tan
_et al._ (2013); He  _et al._ (2014).

Finally, we turn our attention to the valence bands at 4–6 eV binding energy,
which provide a complementary understanding of the observed surface
terminations. In Fig. 3(d), we present high resolution measurements of the
valence band using an $\sim$80 $\mu$m beamspot, observing two bright and sharp
bands around the M point, connected to band minima with weaker intensity at
the $\Gamma$ point. By comparing these dispersions with DFT-based
calculations, as shown in Fig. 4(a-c), we assign these as O $2p$ bands.
Although As $4p$ states are also expected at these binding energies, those
bands tend to be only weakly observed in photoemission measurements due to
high-energy interactions in the FeAs layer Evtushinsky  _et al._ (2017).

When we look at the distinct valence band data from the two terminations, we
observe that the LaO termination (Fig. 3(e)) has the O $2p$-derived valence
bands at generally higher binding energy than the FeAs termination (Fig.
3(f)). This can be naturally explained, since the heavy electron doping of the
LaO layer pushes the surface layer oxygen states to deeper binding energies.
The overall downward shift of the oxygen bands on the LaO surface is further
reproduced by DFT slab calculations on FeAs and LaO-terminated surfaces, shown
in Fig. 4(b) and (c). This shift independently confirms the heavily electron
doped nature of the LaO surface state. Meanwhile, the calculations in Fig.
4(b) show that the oxygen states in the subsurface LaO layer beneath the FeAs
termination have dispersions and binding energies much closer to the bulk
ones, since the FeAs layer screens the surface charge quite effectively.

Figure 4: DFT electronic structure of LaFeAsO \- simulations of the band
structure highlighting the oxygen (red) dominated valence bands for a) a bulk
paramagnetic crystal, b) an FeAs terminated surface, c) a LaO terminated
surface. In (b,c) we only colour the oxygen weight arising from the closest
layer to the surface.

##  IV Discussion

In this work, using nano-ARPES as a local probe method, we established a
strong dependence of the electronic structure with surface termination due to
significant charge transfer at the polar surface of the LaFeAsO crystal.
Because of this, it is not possible to examine the bulk electronic structure
using surface sensitive techniques like STM and ARPES. On the contrary, we can
only observe two very different surface electronic structures with opposite
hole or electron doping. Our results are in good qualitative agreement with
termination dependent STM data Zhou  _et al._ (2011b).

More importantly, it explains, why unlike in non-polar pnictides Rosenthal
_et al._ (2014); Zhang  _et al._ (2012); Yi  _et al._ (2012), the band folding
corresponding to the antiferromagnetic (AFM) order below $T_{N}$ $\sim$135 K
was not been detected in LaFeAsO by ARPES. Although we were able to separately
measure the electronic structure of the FeAs surface termination and the FeAs
subsurface under an LaO termination, neither of these are truly representative
of the bulk electronic structure of LaFeAsO. Moreover, direct photoemission
from deeper layers is suppressed beyond detection due to the short escape
depth of low energy photoelectrons. Therefore we conclude that the true bulk
electronic structure of LaFeAsO, including the expected electronic
reconstruction in the AFM phase, cannot be accessed by ARPES, at least in the
VUV photon energy range. Instead, based on this study, we can interpret with
confidence the previous ARPES data for LaFeAsO obtained with conventional
(large) beam spot Lu  _et al._ (2009); Liu  _et al._ (2010a); Yang  _et al._
(2010); Zhang  _et al._ (2016); Huh  _et al._ (2020) as a superposition of
signals from two different surface terminations.

As we have shown, the charge transfer induced surface doping is large enough
to suppress the AFM phase for the parent compound. So one could ask the
question, what are the implications of this effect on the superconductivity in
the “1111” oxopnictide family? Despite their surface-sensitivity, both STM
Fasano  _et al._ (2010) and ARPES Liu  _et al._ (2008) have reported an
opening of the superconducting gap in fluorine doped single crystals. While
STM measures the spectral function integrated along all pockets of the Fermi
surface, ARPES provides direct information on the momentum dependence of the
superconducting order parameter. A detailed momentum dependent ARPES study
Kondo  _et al._ (2008) revealed an almost isotropic superconducting gap for
the large hole-like band around the $\Gamma$-point. It is clear from our
results, this particular Fermi surface belongs to the FeAs termination.
Because FeAs termination electronic structure is strongly holed doped, it is
unlikely that this top layer is superconducting by itself. A much more
plausible scenario is that surface state superconductivity is induced by
proximity to the the superconducting bulk crystal, an intriguing prospect
which merits further investigation.

There are numerous analogies between our work on LaFeAsO, and recent work on
layered delafossite materials such as PdCoO2, where there has been interest in
the novel electronic states that appear at the polar surfaces after cleaving
Mazzola  _et al._ (2018). In fact, for any ‘naturally heterostructured’
material consisting of alternating layers, the problem of multiple possible
surface terminations contributing to ARPES spectra with a conventional beam
spot is widespread. However, the nano-ARPES methodology that we use here, and
in other recent experiments Mazzola  _et al._ (2018); Iwasawa  _et al._ (2018,
2019) opens the door to obtaining single-termination spectra.

##  V Conclusions

Using angular-resolved photoemission spectroscopy with a sub-micron focused
beam spot (nano-ARPES) we measured the electronic structures for both LaO and
FeAs terminations in LaFeAsO. While large bulk of the ARPES studies is
concentrated on a narrow energy window around the Fermi level, we have shown
that the detailed combined analysis including the nano-ARPES photoemission
signal from much higher binding energies allows for a direct decomposition of
the observed intensity into the contributions from different surface layers.
Our data shows different Fermi surfaces, valence band electronic structure and
core-level spectra for different surface terminations. Our experiments confirm
the theoretical calculations and explain controversy of the previous
conventional ARPES measurements. These results give direct evidence for the
surface-driven electronic structure reconstruction in LaFeAsO, including
formation of the termination-dependent surface states at the Fermi level. The
nano-ARPES technique used to study this prototypical “1111” compound, has and
can be successfully applied to study different polar materials with different
surface terminations.

###### Acknowledgements.

##  VI Acknowledgments

We thank Diamond Light Source for access to beamline I05 (proposal numbers
SI15074 and SI23890) that contributed to the results presented here. LCR
acknowledges funding from the Royal Commission for the Exhibition of 1851. The
work at IFW was supported by the Deutsche Forschungsgemeinschaft (DFG) through
the Priority Programme SPP1458. SA thanks the DFG for funding (AS 523\4-1 &
523\3-1).

## References

  * Kamihara  _et al._ (2008) Y. Kamihara, T. Watanabe, M. Hirano, and Hideo Hosono, “Iron-Based Layered Superconductor La[O1-xFx]FeAs (x = 0.05-0.12) with Tc = 26 K,” Journal of the American Chemical Society 130, 3296–3297 (2008).
  * Zhi-An  _et al._ (2008) Ren Zhi-An, Lu Wei, Yang Jie, Yi Wei, Shen Xiao-Li, Zheng-Cai, Che Guang-Can, Dong Xiao-Li, Sun Li-Ling, Zhou Fang, and Zhao Zhong-Xian, “Superconductivity at 55 K in Iron-Based F-Doped Layered Quaternary Compound Sm[O1-xFx]FeAs,” Chinese Physics Letters 25, 2215–2216 (2008).
  * Wang  _et al._ (2008) Cao Wang, Linjun Li, Shun Chi, Zengwei Zhu, Zhi Ren, Yuke Li, Yuetao Wang, Xiao Lin, Yongkang Luo, Shuai Jiang, Xiangfan Xu, Guanghan Cao, and Zhu'an Xu, “Thorium-doping–induced superconductivity up to 56 K in Gd1-xThxFeAsO,” EPL (Europhysics Letters) 83, 67006 (2008).
  * Lu  _et al._ (2009) D.H. Lu, M. Yi, S.-K. Mo, J.G. Analytis, J.-H. Chu, A.S. Erickson, D.J. Singh, Z. Hussain, T.H. Geballe, I.R. Fisher, and Z.-X. Shen, “ARPES studies of the electronic structure of LaOFe(P,As),” Physica C: Superconductivity 469, 452 – 458 (2009).
  * Liu  _et al._ (2010a) Chang Liu, Yongbin Lee, A. D. Palczewski, J.-Q. Yan, Takeshi Kondo, B. N. Harmon, R. W. McCallum, T. A. Lograsso, and A. Kaminski, “Surface-driven electronic structure in LaFeAsO studied by angle-resolved photoemission spectroscopy,” Phys. Rev. B 82, 075135 (2010a).
  * Yang  _et al._ (2010) L. X. Yang, B. P. Xie, Y. Zhang, C. He, Q. Q. Ge, X. F. Wang, X. H. Chen, M. Arita, J. Jiang, K. Shimada, M. Taniguchi, I. Vobornik, G. Rossi, J. P. Hu, D. H. Lu, Z. X. Shen, Z. Y. Lu, and D. L. Feng, “Surface and bulk electronic structures of LaFeAsO studied by angle-resolved photoemission spectroscopy,” Phys. Rev. B 82, 104519 (2010).
  * Zhang  _et al._ (2016) P. Zhang, J. Ma, T. Qian, Y. G. Shi, A. V. Fedorov, J. D. Denlinger, X. X. Wu, J. P. Hu, P. Richard, and H. Ding, “Disentangling the surface and bulk electronic structures of LaOFeAs,” Phys. Rev. B 94, 104517 (2016).
  * Huh  _et al._ (2020) S S Huh, Y S Kim, W S Kyung, J K Jung, R Kappenberger, S Aswartham, B Büchner, J M Ok, J S Kim, C Dong, J P Hu, S H Cho, D W Shen, J D Denlinger, Y K Kim, and C Kim, “Momentum dependent $d_{xz/yz}$ band splitting in LaFeAsO,” Scientific Reports 10, 19377 (2020).
  * Fedorov  _et al._ (2019) Alexander Fedorov, Alexander Yaresko, Erik Haubold, Yevhen Kushnirenko, Timur Kim, Bernd Büchner, Saicharan Aswartham, Sabine Wurmehl, and Sergey Borisenko, “Energy scale of nematic ordering in the parent iron-based superconductor ${\mathrm{BaFe}}_{2}{\mathrm{As}}_{2}$,” Phys. Rev. B 100, 024517 (2019).
  * Watson  _et al._ (2019) Matthew D Watson, Pavel Dudin, Luke C Rhodes, Daniil V Evtushinsky, Hideaki Iwasawa, Saicharan Aswartham, Sabine Wurmehl, Bernd Büchner, Moritz Hoesch, and Timur K Kim, “Probing the reconstructed Fermi surface of antiferromagnetic BaFe2As2 in one domain,” npj Quantum Materials 4, 36 (2019).
  * Watson  _et al._ (2018) Matthew D. Watson, Saicharan Aswartham, Luke C. Rhodes, Benjamin Parrett, Hideaki Iwasawa, Moritz Hoesch, Igor Morozov, Bernd Büchner, and Timur K. Kim, “Three-dimensional electronic structure of the nematic and antiferromagnetic phases of NaFeAs from detwinned angle-resolved photoemission spectroscopy,” Phys. Rev. B 97, 035134 (2018).
  * Zhou  _et al._ (2011a) Xiaodong Zhou, Cun Ye, Peng Cai, Xiangfeng Wang, Xianhui Chen, and Yayu Wang, “Quasiparticle Interference of ${C}_{2}$-Symmetric Surface States in a LaOFeAs Parent Compound,” Phys. Rev. Lett. 106, 087001 (2011a).
  * Rhodes  _et al._ (2020) Luke C. Rhodes, Matthew D. Watson, Amir A. Haghighirad, D. V. Evtushinsky, and Timur K. Kim, “Revealing the single electron pocket of FeSe in a single orthorhombic domain,” Phys. Rev. B 101, 235128 (2020).
  * Eschrig  _et al._ (2010) Helmut Eschrig, Alexander Lankau, and Klaus Koepernik, “Calculated cleavage behavior and surface states of LaOFeAs,” Physical Review B 81, 155447 (2010).
  * Rhea Kappenberger and Saicharan Aswartham and Francesco Scaravaggi and Christian G.F. Blum and Mihai I. Sturza and Anja U.B. Wolter and Sabine Wurmehl and Bernd Büchner (2018) Rhea Kappenberger and Saicharan Aswartham and Francesco Scaravaggi and Christian G.F. Blum and Mihai I. Sturza and Anja U.B. Wolter and Sabine Wurmehl and Bernd Büchner, “Solid state single crystal growth of three-dimensional faceted LaFeAsO crystals,” Journal of Crystal Growth 483, 9 – 15 (2018).
  * Hoesch  _et al._ (2017) M. Hoesch, T. K. Kim, P. Dudin, H. Wang, S. Scott, P. Harris, S. Patel, M. Matthews, D. Hawkins, S. G. Alcock, T. Richter, J. J. Mudd, M. Basham, L. Pratt, P. Leicester, E. C. Longhi, A. Tamai, and F. Baumberger, “A facility for the analysis of the electronic structures of solids and their surfaces by synchrotron radiation photoelectron spectroscopy,” Review of Scientific Instruments 88, 013106 (2017).
  * Giannozzi  _et al._ (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H-V. Nguyen, A. Otero de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, “Advanced capabilities for materials modelling with quantum ESPRESSO,” Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Liu  _et al._ (2010b) Haiyun Liu, G. F. Chen, Wentao Zhang, Lin Zhao, Guodong Liu, T.-L. Xia, Xiaowen Jia, Daixiang Mu, Shanyu Liu, Shaolong He, Yingying Peng, Junfeng He, Zhaoyu Chen, Xiaoli Dong, Jun Zhang, Guiling Wang, Yong Zhu, Zuyan Xu, Chuangtian Chen, and X. J. Zhou, “Unusual Electronic Structure and Observation of Dispersion Kink in CeFeAsO Parent Compound of FeAs-based Superconductors,” Phys. Rev. Lett. 105, 027001 (2010b).
  * Scaravaggi  _et al._ (2021) F. Scaravaggi, S. Sauerland, L. Wang, R. Kappenberger, P. Lepucki, A. P. Dioguardi, X. Hong, F. Caglieris, C. Wuttke, C. Hess, H.-J. Grafe, S. Aswartham, S. Wurmehl, R. Klingeler, A. U. B. Wolter, and B. Büchner, “Revisiting the phase diagram of ${\mathrm{LaFe}}_{1-x}{\mathrm{Co}}_{x}\mathrm{AsO}$ in single crystals by thermodynamic methods,” Phys. Rev. B 103, 174506 (2021).
  * Tan  _et al._ (2013) Shiyong Tan, Yan Zhang, Miao Xia, Zirong Ye, Fei Chen, Xin Xie, Rui Peng, Difei Xu, Qin Fan, Haichao Xu, _et al._ , “Interface-induced superconductivity and strain-dependent spin density waves in FeSe/SrTiO3 thin films,” Nature Materials 12, 634–640 (2013).
  * He  _et al._ (2014) Junfeng He, Xu Liu, Wenhao Zhang, Lin Zhao, Defa Liu, Shaolong He, Daixiang Mou, Fangsen Li, Chenjia Tang, Zhi Li, _et al._ , “Electronic evidence of an insulator–superconductor crossover in single-layer FeSe/SrTiO3 films,” Proceedings of the National Academy of Sciences 111, 18501–18506 (2014).
  * Evtushinsky  _et al._ (2017) D. V. Evtushinsky, A. N. Yaresko, V. B. Zabolotnyy, J. Maletz, T. K. Kim, A. A. Kordyuk, M. S. Viazovska, M. Roslova, I. Morozov, R. Beck, S. Aswartham, L. Harnagea, S. Wurmehl, H. Berger, V. A. Rogalev, V. N. Strocov, T. Wolf, N. D. Zhigadlo, B. Büchner, and S. V. Borisenko, “High-energy electronic interaction in the $3d$ band of high-temperature iron-based superconductors,” Phys. Rev. B 96, 060501 (2017).
  * Zhou  _et al._ (2011b) Xiaodong Zhou, Cun Ye, Peng Cai, Xiangfeng Wang, Xianhui Chen, and Yayu Wang, “Quasiparticle Interference of ${C}_{2}$-Symmetric Surface States in a LaOFeAs Parent Compound,” Phys. Rev. Lett. 106, 087001 (2011b).
  * Rosenthal  _et al._ (2014) E P Rosenthal, E F Andrade, C J Arguello, R M Fernandes, L Y Xing, X C Wang, C Q Jin, A J Millis, and A N Pasupathy, “Visualization of electron nematicity and unidirectional antiferroic fluctuations at high temperatures in NaFeAs,” Nature Physics 10, 225–232 (2014).
  * Zhang  _et al._ (2012) Y. Zhang, C. He, Z. R. Ye, J. Jiang, F. Chen, M. Xu, Q. Q. Ge, B. P. Xie, J. Wei, M. Aeschlimann, X. Y. Cui, M. Shi, J. P. Hu, and D. L. Feng, “Symmetry breaking via orbital-dependent reconstruction of electronic structure in detwinned NaFeAs,” Phys. Rev. B 85, 085121 (2012).
  * Yi  _et al._ (2012) M Yi, D H Lu, R G Moore, K Kihou, C-H Lee, A Iyo, H Eisaki, T Yoshida, A Fujimori, and Z-X Shen, “Electronic reconstruction through the structural and magnetic transitions in detwinned NaFeAs,” New Journal of Physics 14, 073019 (2012).
  * Fasano  _et al._ (2010) Y. Fasano, I. Maggio-Aprile, N. D. Zhigadlo, S. Katrych, J. Karpinski, and Ø. Fischer, “Local Quasiparticle Density of States of Superconducting ${\mathrm{SmFeAsO}}_{1-x}{\mathrm{F}}_{x}$ Single Crystals: Evidence for Spin-Mediated Pairing,” Phys. Rev. Lett. 105, 167005 (2010).
  * Liu  _et al._ (2008) C Liu, T Kondo, ME Tillman, R Gordon, GD Samolyuk, Y Lee, C Martin, JL McChesney, S Bud’ko, MA Tanatar, _et al._ , “Fermi surface and strong coupling superconductivity in single crystal ${\mathrm{NdFeAsO}}_{1-x}{\mathrm{F}}_{x}$,” arXiv preprint arXiv:0806.2147 (2008).
  * Kondo  _et al._ (2008) Takeshi Kondo, A. F. Santander-Syro, O. Copie, Chang Liu, M. E. Tillman, E. D. Mun, J. Schmalian, S. L. Bud’ko, M. A. Tanatar, P. C. Canfield, and A. Kaminski, “Momentum Dependence of the Superconducting Gap in ${\mathrm{NdFeAsO}}_{0.9}{\mathrm{F}}_{0.1}$ Single Crystals Measured by Angle Resolved Photoemission Spectroscopy,” Phys. Rev. Lett. 101, 147003 (2008).
  * Mazzola  _et al._ (2018) Federico Mazzola, Veronika Sunko, Seunghyun Khim, Helge Rosner, Pallavi Kushwaha, Oliver J. Clark, Lewis Bawden, Igor Marković, Timur K. Kim, Moritz Hoesch, Andrew P. Mackenzie, and Phil D. C. King, “Itinerant ferromagnetism of the Pd-terminated polar surface of PdCoO2,” Proceedings of the National Academy of Sciences 115, 12956–12960 (2018).
  * Iwasawa  _et al._ (2018) Hideaki Iwasawa, Niels B. M. Schröter, Takahiko Masui, Setsuko Tajima, Timur K. Kim, and Moritz Hoesch, “Surface termination and electronic reconstruction in ${\mathrm{YBa}}_{2}{\mathrm{Cu}}_{3}{\mathrm{O}}_{7-\delta}$,” Phys. Rev. B 98, 081112 (2018).
  * Iwasawa  _et al._ (2019) Hideaki Iwasawa, Pavel Dudin, Kyosuke Inui, Takahiko Masui, Timur K. Kim, Cephise Cacho, and Moritz Hoesch, “Buried double CuO chains in ${\mathrm{YBa}}_{2}{\mathrm{Cu}}_{4}{\mathrm{O}}_{8}$ uncovered by nano-ARPES,” Phys. Rev. B 99, 140510 (2019).

##  Appendix A Appendix - Reference core level data

High energy resolution LaFeAsO As $3d$ reference core level spectra were
measured at HR-ARPES end-station with 80 micron beam spot size.

Figure A1: LaFeAsO As $3d$ XPS: (a) photon energy 70 eV and normal incidence,
(b) photon energy 160 eV and normal incidence, (c) photon energy 70 eV and
grazing incidence, (d) photon energy 160 eV and grazing incidence.

Corresponding As $3d$ reference core level spectra were fitted with 3 sets of
doublets P1, P2 and P3.

| (a) 70eV NI | (b) 160eV NI | (c) 70eV GI | (d) 160eV GI
---|---|---|---|---
P1 | 0.128 | 0.084 | 0.257 | 0.183
P2 | 0.276 | 0.103 | 0.416 | 0.327
P3 | 0.596 | 0.813 | 0.327 | 0.490
Table A1: LaFeAsO As $3d$ XPS fit results