printacmref=false ifaamas [AAMAS ’23]Proc. of the 22nd International
Conference on Autonomous Agents and Multiagent Systems (AAMAS 2023)May 29 –
June 2, 2023 London, United KingdomA. Ricci, W. Yeoh, N. Agmon, B. An (eds.)
2023 2023 636 Shanghai AI Laboratory Shanghai China

#  Bringing Diversity to Autonomous Vehicles: An Interpretable Multi-vehicle
Decision-making and Planning Framework

Licheng Wen, Pinlong Cai, Daocheng Fu, Song Mao and Yikang Li
<EMAIL_ADDRESS>
###### Abstract.

With the development of autonomous driving, it is becoming increasingly common
for autonomous vehicles (AVs) and human-driven vehicles (HVs) to travel on the
same roads. Existing single-vehicle planning algorithms on board struggle to
handle sophisticated social interactions in the real world. Decisions made by
these methods are difficult to understand for humans, raising the risk of
crashes and making them unlikely to be applied in practice. Moreover, vehicle
flows produced by open-source traffic simulators suffer from being overly
conservative and lacking behavioral diversity. We propose a hierarchical
multi-vehicle decision-making and planning framework with several advantages.
The framework jointly makes decisions for all vehicles within the flow and
reacts promptly to the dynamic environment through a high-frequency planning
module. The decision module produces interpretable action sequences that can
explicitly communicate self-intent to the surrounding HVs. We also present the
cooperation factor and trajectory weight set, bringing diversity to autonomous
vehicles in traffic at both the social and individual levels. The superiority
of our proposed framework is validated through experiments with multiple
scenarios, and the diverse behaviors in the generated vehicle trajectories are
demonstrated through closed-loop simulations.

###### Key words and phrases:

Autonomous Driving, Driving Behavior, Trajectory Generation, Vehicle Flow

##  1\. Introduction

With the combined efforts of academia and industry, research on autonomous
driving has flourished over the past decade. Nowadays, a growing number of
companies are testing their autonomous vehicles (AVs) on the road, becoming
new participants in the traffic flow besides human-driven vehicles (HVs). The
real-world vehicle flow encompasses a diversity of behaviors at both the
social and individual levels, known as social behavior and driving habit.
Social behavior Wang et al. (2022a); Schwarting et al. (2019) implies how a
vehicle interacts with others. When another vehicle changes lanes, overtakes
or merges, the driver chooses whether to persist in the current movement or
yield according to experience. Drivers also have their own characteristics
when driving. Therefore driving habits Ranney (1994); Wang et al. (2018); Wang
et al. (2022b) is introduced to describe the individual difference, especially
on cruising comfort and driving safety.

Figure 1. Schematic diagram of the proposed decision-making and planning
framework

A framework of the proposed decision-making and planning structure.

The behavioral diversity of AVs is mainly expressed through the decision-
making and planning modules, which are core components of the automated
driving system. The prevailing single-vehicle decision-making and planning
approaches can be divided into two categories: heavy-decision-based frameworks
and light-decision-based ones Katrakazas et al. (2015); Claussmann et al.
(2019). The former methods typically separate decision-making and motion
planning to reduce the computational burden and complexity Urmson et al.
(2008); Montemerlo et al. (2008), whereas the light-decision-based ones
improve the algorithm’s adaptability in corner cases by weakening the
decision-making module and increasing the planner sampling size Howard and
Kelly (2007); Fan et al. (2018); Ding et al. (2019). However, these approaches
struggle to game with other traffic participants in the real world and
therefore exhibit an insufficient understanding of social interactions.

As vehicle-to-everything (V2X) technologies take hold in self-driving
vehicles, centralized decision-making and trajectory planning for multiple
vehicles in traffic become possible. The centralized planner addresses the
complex interactions in the vehicle flow effectively, including competition
and collaboration. There are several types of methods for multi-vehicle
decision and planning. Optimization-based methods Sadigh et al. (2016); Fisac
et al. (2019); Hang et al. (2022) provides good modelling of game-theoretic
interactions between vehicles. Methods based on the Monte Carlo tree search
Browne et al. (2012); Lenz et al. (2016); Li et al. (2022) employ a tree
structure to manage the interactions and considerably improve the algorithm’s
efficiency by sampling vehicles’ action space. However, these methods output
the long-term control variables directly, which can be risky when there are
HVs on the road. Learning-based approaches Casas et al. (2021); Peng et al.
(2021); Chen and Krähenbühl (2022); Sun et al. (2022) obtain vehicles’ driving
strategies from open datasets and avoid collisions by predicting other cars’
trajectories. The trajectories generated by such methods fail to express
driving intentions explicitly and are therefore poorly understood by
passengers or other road users. They also rely heavily on datasets.

In addition, popular open-source autonomous driving simulators provide
solutions for background traffic generation. SUMO Lopez et al. (2018) applies
a car-following model and a lane-changing model to each vehicle. The
trajectories produced by this method can avoid collisions but lack vehicle
kinematic constraints. CARLA Dosovitskiy et al. (2017) adopts a traffic
manager module to generate commands for vehicles according to the current
simulation state. However, this approach applies rigid behavioral models,
which leads to excessively conservative and non-diverse trajectories.

We believe that modeling behavioral diversity provides a better understanding
of real-world traffic flows and leads to a more credible autonomous-driving
simulation. We propose a hierarchical multi-vehicle decision and planning
framework, where the high level makes joint decisions for all vehicles in the
scenario, and the low level generates kinematically trajectories for each
controlled AV. The key features of our framework can be summarized as DII:

  * •

Diversity: By adding a vehicle cooperation factor in the decision module and
introducing the trajectory weight set to the planning module, our framework
produces solutions with social and individual-level diversity.

  * •

Interpretability: The decision module in our framework produces interpretable
action sequences for each AV to keep the passenger on board aware of the
current situation and explicitly communicate intentions to the surrounding HVs
through brake lights and turn signals.

  * •

Interactivity: We include all vehicles in the traffic flow when making
decisions and can therefore explore complex interactions between them.
Moreover, the framework ensures a real-time response to the dynamic
environment through high-frequency replanning.

The remainder of this paper is organized as follows. Section 2 describes the
framework of our proposed methodology. Section 3 introduces the details of the
framework, including multi-vehicle Monte Carlo tree search considering social
behavior and parallel trajectory planner incorporating driving habits. Section
4 presents the performance test along with closed-loop simulations and
discusses the results. Finally, the conclusion and future work are presented
in Section 5.

##  2\. Framework

Consider a set $\mathcal{V}=\\{V_{1},V_{2},\cdots,V_{N}\\}$ containing $N$
vehicles including AVs and HVs. There are $K$ AVs available for centralized
planning, denoted as
$\mathcal{V_{C}}=\\{V_{1},\cdots,V_{k}\\}\subseteq\mathcal{V}$. The remaining
are AVs and all HVs are considered as uncontrolled vehicles that make
decisions and motion planning independently.

Our task is to generate feasible trajectories with behavioral diversity for
each controlled vehicle $V_{i}\in\mathcal{V_{C}}$ and collaborating with
uncontrolled vehicles $V_{j}\in\\{\mathcal{V}\setminus\mathcal{V_{C}}\\}$ in
the environment. We introduce a two-stage multi-vehicle framework containing
both decision-making and trajectory planning modules. The input to the
framework is a perceivable environment which contains road conditions, route
information, vehicle status and predictions of uncontrolled vehicles. The
schematic diagram of our proposed framework is illustrated in Figure 1, where
vehicles with IDs 2 and 3 are AVs controlled by our framework while vehicle 1
is an uncontrolled vehicle. The details of the two stages in our proposed
framework are described below.

In the first stage, we propose a decision module to address the social
interactions among vehicles. This module generates long-term, coarse decisions
for all vehicles jointly and simultaneously. The environment is constructed on
the Frenét frame Rovenski (2006) for decision-making, which allows the
decision module to focus on inter-vehicle interactions without considering
road constraints. Since actions like continuous lane-changing and overtaking
take a relatively long time to complete, the decision module is required to be
forward-looking and consistent. The decision module solves the interaction
game between vehicles by making centralized joint decisions on the vehicle
flow. It can generate diverse social behaviors in the flow via setting
different cooperation factors in the reward function. The decision module
generates a temporal sequence of discrete actions for each controlled vehicle.
Since each action in the sequence has a physical meaning, such as changing
lanes, accelerating, maintaining speed, etc., the output of our decision
module is interpretable.

In the second stage, the planning module receives discrete action sequences
from the first stage and generates continuous kinematic-feasible trajectories
for each controlled vehicle. Thanks to the forward-looking nature of our
front-end decision module, the trajectory planning module only needs to
generate trajectories within a short planning horizon. Our planning module
adopts a distributed parallel architecture, i.e. running a planner
independently for each vehicle, which more closely resembles human driving in
the real world. For each parallel planner, the leading part of its received
action sequence is selected as the guidance to generate a continuous
trajectory that conforms to the vehicle kinematic constraints in Cartesian
space. This module reflects the driving habits of different drivers through
planner’s selection strategy for various trajectories. It is worth emphasizing
that the planner would maintain a high replanning frequency and continuously
predicts surrounding vehicles to avoid potential collisions.

Our proposed two-stage framework provides an interpreted result with physical
meaning for each module,allowing the AV to inform passengers of the current
situation and to explicitly communicate its own intentions to the HV
surrounding them. Meanwhile, it guarantees bottom-line security. The planning
module realizes collision risks by its prediction of other vehicles and
generates collision-free trajectories, even if the decision module produces
some impractical action sequences.

##  3\. Methodology

###  3.1. Joint Multi-vehicle Decision-making

####  3.1.1. Basic Monte-Carlo Tree Search Algorithm

Monte-Carlo Tree Search, also known as MCTS, is a heuristic search algorithm
for the decision process which analyses the search space for the most
promising actions based on the Monte Carlo sampling method. The general MCTS
algorithm Browne et al. (2012) includes four phases in one search iteration:
selection, expansion, simulation, and back-propagation, as shown in Figure 2.

Figure 2. Procedure of the Monte-Carlo Tree Search

A demo of the structure of the MCTS algorithm.

The selection phase is performed based on an existing search tree. The
algorithm will recursively select the most valuable child node using a tree
policy until reaching a non-fully expanded node. Then, the selected node is
expanded with one of the remaining possible actions to get a new leaf node.
From this node, the algorithm will use a default policy (also called roll-out)
to perform the simulation until it reaches a terminal node. At the back-
propagated phase, the algorithm will extract the rewards from the simulation
sample and then use them to update the average rewards through all nodes that
have been traversed. After reaching predefined maximum iterations, the
solution can be obtained by greedily selecting the child node with a maximum
reward at each layer.

The standard UCT algorithm Browne et al. (2012) is applied in order to address
the exploration-exploitation dilemma in the selection phase. The UCT value of
child node $j$ can be calculated by:

$\displaystyle UCT$ $\displaystyle=$
$\displaystyle\overline{X}_{j}+2C_{p}\sqrt{\frac{2\ln n}{n_{j}}}$ (1)

where $\overline{X}_{j}$ is the average reward of node $j$, $n$ is the number
of times the parent node of $j$ has been visited, $n_{j}$ is the number of
times child j has been visited, and $C_{p}>0$ is a constant to balance
exploration and exploitation. The reward distribution of each node should lie
between 0 and 1.

The vanilla MCTS algorithm typically makes decisions for a single agent.
However in our framework, the decisions are required for all vehicles in the
flow at the same time. A trivial approach would be to perform a separate MCTS
search for each vehicle. However performing MCTS searches separately for each
vehicle requires accurate predictions of other vehicles’ behavior before
making a decision,. This violates Vapnik’s Principle Vapnik (2006) that”When
solving a problem, don’t solve a harder problem as an intermediate step.”, as
modeling and predicting the behavior of other agents is a harder problem.

####  3.1.2. Metanode in the MCTS

Several multi-vehicle decision methods based on the Monte Carlo tree search
have been developed. Lenz et al. Lenz et al. (2016) propose a method for
generating anticipatory and cooperative behavior of AVs, but its direct output
of motion via MCTS is poorly robust to dynamic environments, especially with
HVs mixing in the road. Another approach Li et al. (2022) sets priorities for
the vehicles and makes decisions from highest priority to lowest, with the
subsequent agent required to satisfy the decisions completed by the preceding
ones. However, such approaches fail to express the different social behaviors
in the flow, as the low-priority vehicles only passively respond to the
decisions made by the high-priority ones. Besides, it’s never easy to assign
priority to every single vehicle in the traffic flow.

Instead, our approach takes the HVs in the traffic flow into account and there
is no priority between vehicles in the decision-making process. We replace the
node in the general MCTS method with a metanode that can generate multi-
vehicle actions simultaneously. The metanode’s extension is shown in Figure 3,
where vehicle 2 and 3 are vehicles controlled by our framework while vehicle 1
is an uncontrolled vehicle.

Figure 3. A metanode at time step $t$ in the MCTS

A metanode at time step $t$ in the MCTS search.

For a metanode extended to time step $t$, it receives all controlled vehicle’s
state at the current moment. In particular, the state for vehicle $V_{i}$ at
time step $t$ can be described as $\mathbf{x}_{i}^{t}=[s,d,v]^{T}$, $s$ and
$d$ are the Frenét coordinates of the vehicle and $v$ is the longitudinal
velocity. Similar to Lenz et al. (2016), we define a set of possible actions
for each controlled vehicle to take in one time step. As listed in Table 1,
there are five actions - maintain speed (KS), accelerate (AC), decelerate
(DC), change lanes to the left (LCL) and to the right (LCR). By selecting
different actions, the vehicle can then calculate its corresponding state at
the next time step. As for uncontrolled vehicles in the flow, we assume they
performing a lane-keeping (KL) action. This action dynamically adjusts their
velocity and always maintains an appropriate distance from the leading vehicle
during the decision-making process.

Table 1. Actions available for vehicles in one time step

Action | State change | Description
---|---|---
KS | $\left[v\Delta t,0,0\right]$ |  Maintain the current velocity
AC | $\left[v\Delta t+\frac{a_{acc}}{2}\Delta t^{2},0,a_{acc}\Delta t\right]$ |  Increase the current velocity with a fixed acceleration $a_{acc}$
DC | $\left[v\Delta t-\frac{a_{dec}}{2}\Delta t^{2},0,-a_{dec}\Delta t\right]$ |  Decrease the current velocity with a fixed deceleration $a_{dec}$
LCL | $\left[v\Delta t,-\Delta d,0\right]$ |  Make a partial left lane change with a width of $\Delta d$
LCR | $\left[v\Delta t,\Delta d,0\right]$ |  Make a partial right lane change with a width of $\Delta d$
KL | - |  Keep the current lane (Only for uncontrolled vehicles)

Since decisions need to be made for all vehicles simultaneously within a
metanode, the output of the metanode would combine each vehicle’s possible
behavior. As a result, the search tree grows exponentially when the MCTS
gradually evaluates the whole interaction decision space. Finding the best
solution that is computationally feasible becomes impossible in practice.
Therefore pruning the expansion of the metanode proves necessary.

####  3.1.3. Pruning in MCTS

Although each vehicle has five optional actions, not every one of them is
feasible and some may lead to a potential collision. Obviously, vehicles on
either side of the road cannot continue to make lane-changing actions to over
borderlines. According to Swaroop and Rajagopal (2001), vehicles in traffic
should maintain a shortest safe distance $D_{s}$ from its leading vehicle.
$D_{s}$ consisted of two parts: reaction distance and minimum headway, which
are defined as:

$\displaystyle D_{s}$ $\displaystyle=$ $\displaystyle
v\cdot\tau+MTH\cdot\Delta v$ (2)

where $\tau>0$ and $MTH>0$ are both time constants, representing reaction time
and minimum time headway respectively. $v$ denotes vehicle longitude velocity
and $\Delta v=v-v_{l}$ denotes the velocity difference between the self and
lead vehicles. Further, the velocity limit that any single vehicle in the flow
should satisfy is:

$\displaystyle v$ $\displaystyle\in$ $\displaystyle\left[v_{f}-\frac{\Delta
s_{f}-\tau\cdot v_{f}}{{MTH}},\min(\frac{{MTH}\cdot v_{l}+\Delta
s_{l}}{\tau+{MTH}},\frac{\Delta s_{l}}{\tau})\right]$ (3)

where $\Delta s_{l}$ and $\Delta s_{f}$ indicate the gap between the current
vehicle and its leading/following vehicle, respectively. $v_{l}$ and $v_{f}$
denote the velocity of the leading and following vehicle. It is worth noting
that for the upper limit, even if the vehicle’s current velocity is slower
than the leading vehicle, it still needs to maintain a minimum reaction
distance.

In the metanode at time step $t$, when vehicle $i$ takes an action $a$ from
Table 1 and the resulting state $\mathbf{x}_{i}^{t+1}$ does not satisfy
Equation (3), the action $a$ becomes invalid and the corresponding state
$\mathbf{x}_{i}^{t+1}$ will be abandoned. Also for all valid state
combinations within the metanode, the algorithm checks for collisions between
new states in each combination and for conflicts in the process of reaching
the states (e.g. two vehicles changing into each other’s lanes). If a conflict
occurs, the action combination should not be extended for the metanode at time
step $t+1$. Through these techniques, the decision space is significantly
reduced and the search efficiency is guaranteed. As Figure 3 shows, these
three vehicles originally generated $5^{2}=25$ possible action combinations,
and only 4 valid combinations remain after pruning. The efficiency of the
algorithm improves more significantly when facing complex traffic situations.

####  3.1.4. Reward Function with Social Behavior

During the simulation phase of MCTS, the algorithm continues to construct the
metanode, but only one valid action combination is selected to extend the next
time step until the simulation terminates. Each simulation can be terminated
either when all vehicles have traveled a certain distance in their target
lane, or when the simulation reaches the maximum time step with any vehicles
not completing their intentions. Then a reward for this simulation should be
calculated from the result.

We first calculate the reward $R_{i}$ separately for each vehicle $i$ in the
traffic flow. The reward $R_{i}$ has three parts, namely driving in the target
lane, driving in the lane’s center line and maintaining the consistency of the
actions. Besides, the reward should meet in the range of 0 to 1.

For the vehicle flow to exhibit a diversity of social behaviors, the
cooperative tendency of the vehicle needs to be characterized by the reward.
Similar to Schwarting et al. (2019), we introduce a behavior reward
considering both reward to self and reward to others:

$\displaystyle R$ $\displaystyle=$ $\displaystyle
R_{\rm{self}}\leavevmode\nobreak\ +\leavevmode\nobreak\
\gamma\leavevmode\nobreak\ R_{\rm{other}}$ (4)

where $\gamma\in[0,1]$ being the cooperation factor. $\gamma=0$ implies the
vehicle is egoistic and takes no account of the behavior of other vehicles,
whereas $\gamma=1$ denotes the vehicle treats other vehicles equally important
to itself when making decisions. Each vehicle $i$ in the flow possesses its
unique factor $\gamma_{i}$.

Since each metanode in the search tree handles the actions of the whole
vehicle flow, the reward of the simulation is obtained by combining all
vehicles’ behavior rewards:

$\displaystyle X_{\rm{flow}}$ $\displaystyle=$
$\displaystyle\frac{1}{K}\sum_{1\leq i\leq K}\frac{R_{i}+\gamma_{i}\sum_{j\neq
i}{R_{j}}}{1+(K-1)\gamma_{i}}$ (5)

Since $R_{i}\in[0,1]$, the distribution of $X_{\rm{flow}}$ is also guaranteed
to stay in $[0,1]$. Finally, reward $X_{\rm{flow}}$ is used to update the
average rewards of all selected metanodes in the back-propagation phase.

###  3.2. Parallel Vehicle Trajectory Planning

The decision module gives a coarse-grained action sequence for each vehicle,
but these actions take no account of vehicle kinematics and road curvature.
They thus cannot be performed by the vehicles directly. Meanwhile, the
generated discrete actions may still expose the vehicles to potential
collision risk. We employ a trajectory planning module after decision-making
to generate continuous, kinematic feasible trajectories and ensure the bottom-
line safe of the whole framework.

####  3.2.1. Conversion between Cartesian and Frenét Frame

As mentioned earlier, we perform decision-making in the Frenét coordinate
system to simplify the calculations. However, trajectories generated by the
planning module require converting to the global Cartesian frame for verifying
the vehicle’s kinematic constraints. Here we briefly describe the conversion
from the Frenét frame to the Cartesian coordinates.

A Frenét coordinate system is constructed via a reference path generated by
the center-line of the current lane, it is denoted by the arc-length parameter
$s$ along the reference path and the lateral offset $d$. $\dot{s}$ and
$\dot{d}$ are the derivative with respect to time. For a vehicle with the
Frenét coordinates $\left[s,\dot{s},d,\dot{d}\right]$, we can calculate its
position $(x,y)$, velocity $v$, orientation $\theta$ in the Cartesian frame
according to the transformation proposed in Werling et al. (2010):

$\displaystyle\left\\{\begin{array}[]{lc}x=&x_{r}-d\sin\theta_{r}\\\
y=&y_{r}+d\cos\theta_{r}\\\
v=&\sqrt{\left[\left(1-\kappa_{r}d\right)\dot{s}\right]^{2}+\dot{d}^{2}}\\\
\theta=&\arcsin\left(\frac{\dot{d}}{v}\right)+\theta_{r}\\\
\end{array}\right.$ (10)

where $(x_{r},y_{r})$, $\theta_{r}$ and $\kappa_{r}$ denote the position, the
orientation and curvature of the reference path at length $s$.

####  3.2.2. Feasible Trajectory Planner

As mentioned in Section 2, in the planning module we employ a parallel
architecture for trajectory planning of each vehicle, which not only improves
the runtime efficiency of the framework but also resembles the way humans
drive in the real world. The pipeline of the single-vehicle trajectory planner
is shown in Figure 4.

Figure 4. The pipeline of a single-vehicle trajectory planner

This diagram demonstrates the pipeline of a vehicle trajectory planner.

For an AV $V_{i}\in\mathcal{V_{C}}$ under our control, it receives an action
sequence from the decision module. Then one of the parallel planners takes
over. The planner first predicts the future trajectories of other vehicles in
the flow. Then it splits the different actions in the sequence and plans a
sub-trajectory for each action segment. Finally, the planner connects all the
sub-trajectories to obtain the final feasible and continuous vehicle
trajectory.

Our sub-trajectory planner is similar to the one proposed in Werling et al.
(2010), but with some differences in the pipeline. First, the algorithm
determines the range of sub-trajectory’s target state based on the current
action segment. Although the state change for each action is a constant when
making decisions as shown in Table 1, we extend it to a continuous region of
target states for sub-trajectory planning. Additionally, for vehicles with
lane-keeping actions, the planner creates a larger target state region to make
rational in-lane decisions. The planner then uniformly samples a number of
target states within the predefined region. We employ quintic polynomial
curves in the Frenét frame to generate the jerk-optimal connections between
the start state and all sampled target states of the sub-trajectory. After
that, we convert each sub-trajectory to Cartesian coordinate by using Equation
10, including its position, velocity and heading angle. All alternative sub-
trajectories are calculated with a cost function, which is detailed in Section
3.2.3. We check whether each sub-trajectory meets the kinematic constraints
and is collision-free with other vehicles in the flow. The sub-trajectory that
satisfies all the constraints with the smallest cost is selected and the final
state of this sub-trajectory is set as the initial state for the next sub-
trajectory.

####  3.2.3. Cost Function considering Driving Habits

For each alternative sub-trajectory, a cost function should be calculated to
evaluate its safety and comfort. In particular, there are several cost terms
in the function:

  * •

Curve energy cost: $J_{i,cur}=\sum\kappa_{i}^{2}$, $\kappa_{i}$ denotes the
curvature of a trajectory point. This penalty trajectory with sharp turns.

  * •

Heading difference cost:
$J_{i,\phi}=\sum\left|\phi_{i}-\phi_{r}\right|^{2}$,$\phi_{i}$ denotes the
orientation angle of the vehicle and $\phi_{r}$ indicates the road’s direction
at the point.

  * •

Out of line cost: $J_{i,out}=\sum\left|d_{i}\right|^{2}$, $d_{i}$ denotes the
lateral offset in Frenét frame. This term is only active for actions except
for lane-changing ones.

  * •

Acceleration cost: $J_{i,acc}=\sum a_{i}^{2}$, $a_{i}$ denotes vehicle’s
acceleration at a trajectory point.

  * •

Jerk cost: $J_{i,jerk}=\sum j_{i}^{2}$, $j_{i}$ denotes vehicle’s jerk
(derivative of acceleration) at a trajectory point. This term and the
acceleration cost both contribute to the comfort of the trajectory.

  * •

Dynamic obstacle cost: $J_{i,obs}$. This term serves to avoid the vehicle from
being too close to the surrounding vehicles, which is described in detail
below.

Dynamic obstacle cost takes the distance between a vehicle following a
trajectory and dynamic obstacles on the road (i.e., surrounding vehicles) into
account. The cost for vehicle $V_{i}$ is defined as:

$\displaystyle J_{i,obs}$ $\displaystyle=$
$\displaystyle\sum_{t=0}^{T}\sum_{\forall j\neq i}J_{i,obs}(t,j)$ (11)

where $J_{i,obs}(t,j)$ characterizes the position of vehicle $V_{j}$ with
respect to the current vehicle $V_{i}$ at time step $t$. As shown in Figure 5,
we introduce an alert zone whose shape is a rectangle surrounding the body of
the vehicle $V_{i}$. The zone has a length of the shortest safe distance
$D_{s}$ defined by Equation 2 in front of the vehicle and a length of 1.5
times the body length behind the vehicle. The width of the alert zone is 1.5
times the width of the body.

Figure 5. The alert zone of vehicle $V_{i}$

This diagram demonstrates three different scenarios when calculating dynamic
obstacle cost.

The position of the vehicle $V_{j}$ in the flow relative to $V_{i}$ can be
divided into three different types: outside the alert zone, inside the alert
zone without collisions, and having collisions with the current position of
$V_{i}$. The vehicle $V_{j}$ is safe outside the alert zone and therefore
cost-free, whereas its direct collision with the trajectory is unacceptable.
In the alert zone, the closer vehicle $V_{j}$ is to $V_{i}$, the higher the
cost. In particular, $J_{i,obs}(t,j)$ is defined as follows:

$J_{i,obs}(t,j)\leavevmode\nobreak\ =\leavevmode\nobreak\
\begin{cases}0,&V_{j}\leavevmode\nobreak\ \rm{outside\leavevmode\nobreak\
alert\leavevmode\nobreak\ zone}\\\
C_{z}(2-\frac{|D|_{x}}{D_{s}+1.5l}-\frac{|D|_{y}}{1.5d}),&V_{j}\leavevmode\nobreak\
\rm{in\leavevmode\nobreak\ alert\leavevmode\nobreak\ zone}\\\
\infty,&V_{j}\leavevmode\nobreak\ \rm{has\leavevmode\nobreak\
collision}\end{cases}$ (12)

where $|D|_{x}$ and $|D|_{y}$ represent the component of the shortest distance
from the vehicle $V_{j}$ to vehicle $V_{i}$ along the length and width of the
alert zone, respectively. $l$ and $d$ denote the length and width of vehicle
body and $C_{z}>0$ is a cost constant.

All the aforesaid terms are summed up as the cost function:

$\displaystyle J_{i}$ $\displaystyle=$ $\displaystyle\sum\leavevmode\nobreak\
w_{(\cdot)}\leavevmode\nobreak\ J_{i,(\cdot)}$ (13)

where $J_{i,(\cdot)}\in\\{J_{i,cur},\leavevmode\nobreak\
J_{i,\phi},\leavevmode\nobreak\ J_{i,out},\leavevmode\nobreak\
J_{i,acc},\leavevmode\nobreak\ J_{i,jerk},\leavevmode\nobreak\ J_{i,obs}\\}$
represents all the cost items mentioned above, and $w_{(\cdot)}$ is the weight
parameter corresponding to each cost item. We observe different drivers in
traffic flow may choose various trajectories depending on their driving
habits. In other words, drivers weigh up these trajectory’s cost terms
differently. For instance, some drivers may be more concerned with driving
comfort, while others prefer to keep a longer safe distance from other
vehicles. Therefore, we define a trajectory weight set $\mathcal{W}$
containing $M$ weight vectors based on the driver’s driving habits:

$\displaystyle\mathcal{W}=\\{\mathbf{w}_{1},\mathbf{w}_{2},\cdots\mathbf{w}_{M}\\}$
(14)

Each weight vector $\mathbf{w}$ contains a unique combination of weights for
all cost terms, $\mathbf{w}=[w_{cur},\leavevmode\nobreak\
w_{\phi},\leavevmode\nobreak\ w_{out},\leavevmode\nobreak\
w_{acc},\leavevmode\nobreak\ w_{jerk},\leavevmode\nobreak\ w_{obs}]$. AVs in
the flow pick different weight vectors to generate their trajectories, in
conjunction with the cooperation factor considering social behavior in the
decision module, bringing diversity to the whole vehicle flow.

##  4\. Experiment

Figure 6. Vehicle trajectories in the different scenarios

.

We utilize real-world driving records from the CitySim dataset Zheng et al.
(2022) to generate the experimental environment, which consisted of three
different scenarios: Freeway, Ramp, and Roundabout, as shown in Figure 6. The
Freeway scenario is the simplest one, containing a four-lane expressway. The
ramp scenario includes a three-lane main road with an on-ramp merging into it.
The roundabout scenario measures the performance of vehicle flows at an
unsignalled two-lane roundabout entrance. We take the vehicle information from
the dataset as our experiment’s initial state of the vehicle flow. The
vehicles about to make a lane change (including entering and exiting the ramp)
and the vehicles directly interacting with them are set as AVs controlled by
the centralized planner. In contrast, other unaffected vehicles are set as
uncontrolled HVs. The vehicles’ lane change intentions are fed into the
framework as input. The proposed framework is implemented entirely in Python
for fast prototyping.

The vehicle size in the flow is set to 5m $\times$ 2m. The constant
$C_{p}=1/\sqrt{2}$ in Equation 1. We set reaction time $\tau=0.5$s and minimum
time headway $MTH=3$s in Equation 2. The fixed acceleration $a_{acc}$ and
deceleration $a_{dec}$ in available actions of Table 1 are set to 0.6m/s2. The
time step is set to 1.5s in the decision module and 0.1s in the planning
module. To generate vehicle flows with diversity, three different types of
driving behavior are defined: aggressive, normal and conservative. The
aggressive vehicle takes less account of the surrounding vehicles when making
decisions and thus has a smaller cooperation factor $\gamma$, along with a
higher tolerance to both obstacle avoidance and acceleration changes. The
opposite is true for the conservative vehicles. The specific style-related
parameters are defined in Table 2.

Table 2. Three types of driving behaviors Type | Aggressive | Normal | Conservative
---|---|---|---
$\gamma$ | 0.1 | 0.5 | 0.9
$w_{i,obs}$ | 2.0 | 4.0 | 6.0
$w_{i,out}$ | 6.5 | 5.0 | 4.0
$w_{i,jerk}$ | 0.7 | 1.0 | 1.2
Table 3. Comparative experiments in the Freeway scenario | 2 Controlled AVs | 3 Controlled AVs | 4 Controlled AVs
---|---|---|---
|  Sequential MCTS |  Ours w/o Planning |  Ours |  Sequential MCTS |  Ours w/o Planning |  Ours |  Sequential MCTS |  Ours w/o Planning |  Ours
Sucess Rate $\uparrow$ | 100% | 100% | 100% | 100% | 100% | 100% | 70% | 90% | 100%
Avg Expand Nodes $\downarrow$ | 1304.1 | 1196.6 | 1148.7 | 1387.6 | 1127.8 | 1037.6 | 1706.6 | 1785.0 | 1527.3
Avg Finish Time (s) $\downarrow$ | 3.57 | 5.75 | 5.68 | 3.57 | 4.95 | 3.50 | 4.58 | 9.29 | 8.93
Min Distance (m) $\uparrow$ | 4.34 | 8.41 | 9.13 | 4.74 | 7.02 | 10.78 | 5.02 | 7.54 | 10.94

Our proposed framework generates promising planning solutions within all three
scenarios. The generated traffic trajectories are shown in Figure 6, where the
controlled AVs in the flow is denoted by colored trajectories and uncontrolled
HVs by grey trajectories. The interpretable actions generated by the decision
module are shown in the figure through the brake lights and turn signals.

We further compare and analyze the performance of our proposed framework and
the diversity of the solutions it generates in the following sections.

###  4.1. Framework Performance

In this subsection, we test the performance of the proposed framework. To
evaluate the effectiveness of our proposed method, we first implemented a
sequential MCTS approach as a comparison, which is used in Li et al. (2022)
(marked as Sequential MCTS). This method assigns a right-of-way priority to
each vehicle in the flow and makes decisions for vehicles in order of priority
at each decision step. In this experiment, vehicles with greater longitudinal
distance possess a higher decision priority. As the previous study did not
consider uncontrolled HVs in the scenario, we employ a car-following model
Treiber et al. (2000) to predict all HVs’ trajectories and feed them into the
sequential MCTS as dynamic obstacles. We also include a decision-only
framework (marked as Ours w/o Planning) in the comparison experiments to
demonstrate the necessity and effectiveness of the planning module in our
proposed framework.

According to Treiber et al. (2000), the longitudinal acceleration in
sequential MCTS is chosen to be [$-1,0,1$]m/s2, and the lateral velocity is
chosen to be [$-1.2,0,1.2$]m/s. We set the decision time step for sequential
MCTS and decision-only framework to 1.0s. The planning horizon for all three
methods is 10s and the maximum iteration number of MCTS is 3000.

The experiments are conducted in the four-lane freeway scenario with the
number of controllable AVs set to 2, 3, and 4. Each method is repeated ten
times under the same experimental setup. We measure each method’s success
rate, the average number of expansion nodes in MCTS, the average time to
finish the vehicle’s intention and the minimum inter-vehicle distance in the
generated solutions. The results are shown in Table 3.

All three methods achieve a 100% success rate in the freeway scenarios with
two and three AVs. However, when the number of controlled AVs reaches four,
sequential MCTS can only attain a success rate of 70%. The decision-only
approach has a 90% success rate, while our framework still manages to pass all
tests. As the sequential MCTS defines the priority for vehicles in the flow,
the higher-priority vehicles do not yield to the lower-priority ones when
making decisions. This results in the high-priority vehicles taking locally
optimal decisions that make it impossible for the subsequent vehicles to
complete their intentions.

The number of extended nodes in the search tree measures the search efficiency
of each method. Our proposed framework has fewer available actions for each
vehicle and adopts various pruning strategies as we mentioned in Section
3.1.3. Experimental results show that our approach expands fewer nodes than
the comparison methods, which means that our framework searches for the
optimal solution faster without wasting actions on infeasible actions.

Figure 7. Three different cases in the closed-loop simulation

This diagram demonstrates the pipeline of a vehicle trajectory planner.

Figure 8. Vehicles’ initial states in closed-loop simulation

This diagram demonstrates the pipeline of a vehicle trajectory planner.

In terms of the quality of the solutions, although the average time for
vehicles to complete their intentions is shorter in the solutions generated by
the sequential MCTS, the minimum inter-vehicle distances of these solutions
are shorter and hence carry more risk. The shortest distance between vehicles
in the solution planned by sequential MCTS is generally less than 5 meters
(the length of the vehicle), which can easily lead to rear-end collisions. Our
decision module, in contrast, achieves a minimum distance of seven to eight
meters thanks to introducing the shortest safe distance $D_{s}$. The planning
module refines the trajectory further, bringing the whole approach to a
minimum distance of around 10 meters, i.e., twice the vehicle length. Besides,
the whole framework has a shorter average intention finishing time compared to
the decision-only approach, even outperforming the sequential MCTS in the case
of 3 vehicles.

###  4.2. Closed-loop Simulation

We conducted several closed-loop simulations in the ramp scenario to
investigate the diversity of the vehicle trajectories regarding social
behavior and driving habits. We also demonstrate the proposed framework’s
responsiveness when facing human drivers’ unanticipated actions.

The closed-loop simulation runs at 10Hz. Our decision-making module has a
horizon of 10.5 seconds, with 1.5 seconds for each action. The horizon of the
planning module is 3 seconds with a time step of 0.1 seconds, and the
replanning cycles are 0.5 seconds. Our test scenario is shown in Figure 8.
Vehicles with the number 2, 4, and 5 are AVs controlled by our framework,
while vehicles with the number 1, 3, and 6 are HVs in the environment.
Vehicles 2 and 5 are driving side by side in the adjacent lanes of the main
road, while vehicle 4 is on the ramp about to merge onto the main road. In
addition, vehicle 2 intends to change its lane to the right, creating a
complex situation where three cars in the flow compete for the right-of-way of
the rightmost lane simultaneously. Figure 7 shows three different cases
arising from this same initial state. The velocity profiles of controlled AVs
in each case are shown on the right of the figure.

In Case 1, vehicle 2 is assigned an aggressive driving behavior, vehicle 4
adopts a conservative driving behavior and vehicle 5 is assigned the normal
behavior. Their style-related parameters are set according to Table 2. The
closed-loop experimental results show that due to its backward initial
position, vehicle 5 slows down to give way to vehicle 4, who maintains its
initial velocity and merges onto the main road. At the moment $T_{1}$, as
vehicle 2 has a more aggressive driving behavior, it is unlikely to yield to
vehicle 4. Vehicle 2 continues to drive in the left front of vehicle 4 and
indicates its intention to change lanes. Finally, vehicle 2 slows down to make
space for vehicle 4 to complete its lane change at $T_{2}$.

We change the driving behavior of the Avs in Case 2. The driving behavior of
vehicle 4 is altered to aggressive and vehicle 2 is adjusted to conservative.
The driving behavior of vehicle 5 remains the same. Vehicle 5 still slows down
to let vehicle 4 merge into the main road at $T_{1}$, but this time vehicle 4
does not yield to vehicle 2’s lane change and chooses to accelerate through
it. Vehicle 2, aware of vehicle 4’s acceleration, chooses to slow down and
completes the lane change action after vehicle 2, which is quite different
from Case 1.

In Case 3 we explore what happens when an uncontrollable HV in the environment
does not cooperate with AV’s behavior. We change vehicle 5 from AV to HV that
is not controlled by our planner and set a larger initial speed for it.
Instead of slowing down to give way to vehicle 2, as our framework expected,
it accelerates into the ramp entrance, creating a quite dangerous situation.
Thanks to the high-frequency replanning in our planning module, vehicle 2 and
vehicle 4 abort their original actions from the decision module when they
realize that vehicle 5 is accelerating in order to avoid the potential
collision. Vehicle 4 performs a deceleration and vehicle 2 keeps its original
lane. Once vehicle 5 has passed, vehicle 2 and 4 then complete their lane-
changing and merging behavior safely. This case demonstrates the importance of
integrating the planning module into the framework, which can quickly respond
to changes in the road and resolve potential collisions between vehicles,
safeguarding the entire vehicle flow.

##  5\. Conclusions

This paper presents a novel hierarchical multi-vehicle decision-making and
planning framework that produces trajectories with diversity, interactivity,
and interpretability in mixed AV and HV scenarios. By introducing the
vehicle’s cooperation factor in the decision module and the trajectory weight
set in the planning module, the trajectories generated by our framework are
diverse at both the social interaction and individual habit levels. Our
framework makes decisions jointly for all vehicles in the flow and ensures a
real-time response to the dynamic environment through a high-frequency re-
planning method. Also, the decision module generates interpretable action
sequences, which are leveraged to inform passengers of the current situation
and explicitly communicate self-intents to the surrounding HVs.

In the further work, our framework can be used to generate background traffic
around the ego car in the autonomous driving simulator. It can also be
utilized to perform data augmentation of the existing traffic datasets with
different driving styles. In addition, our framework is generic in nature,
meaning each module’s internal algorithms can be replaced, provided their
respective inputs and outputs are satisfied.

## References

  * (1)
  * Browne et al. (2012) Cameron B Browne, Edward Powley, Daniel Whitehouse, Simon M Lucas, Peter I Cowling, Philipp Rohlfshagen, Stephen Tavener, Diego Perez, Spyridon Samothrakis, and Simon Colton. 2012.  A Survey of Monte Carlo Tree Search Methods.  _IEEE Transactions on Computational Intelligence and AI in Games_ 4, 1 (mar 2012), 1–43.  https://doi.org/10.1109/TCIAIG.2012.2186810
  * Casas et al. (2021) Sergio Casas, Wenjie Luo, and Raquel Urtasun. 2021.  IntentNet: Learning to Predict Intention from Raw Sensor Data.  CoRL (jan 2021).  arXiv:2101.07907 http://arxiv.org/abs/2101.07907
  * Chen and Krähenbühl (2022) Dian Chen and Philipp Krähenbühl. 2022.  Learning From All Vehicles. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_. 17222–17231.
  * Claussmann et al. (2019) Laurene Claussmann, Marc Revilloud, Dominique Gruyer, and Sébastien Glaser. 2019.  A review of motion planning for highway autonomous driving.  _IEEE Transactions on Intelligent Transportation Systems_ 21, 5 (2019), 1826–1848.
  * Ding et al. (2019) Wenchao Ding, Lu Zhang, Jing Chen, and Shaojie Shen. 2019\.  Safe trajectory generation for complex urban environments using spatio-temporal semantic corridor.  _IEEE Robotics and Automation Letters_ 4, 3 (2019), 2997–3004.  https://doi.org/10.1109/LRA.2019.2923954
  * Dosovitskiy et al. (2017) Alexey Dosovitskiy, German Ros, Felipe Codevilla, Antonio Lopez, and Vladlen Koltun. 2017.  CARLA: An Open Urban Driving Simulator. In _Proceedings of the 1st Annual Conference on Robot Learning_. 1–16.
  * Fan et al. (2018) Haoyang Fan, Fan Zhu, Changchun Liu, Liangliang Zhang, Li Zhuang, Dong Li, Weicheng Zhu, Jiangtao Hu, Hongye Li, and Qi Kong. 2018\.  Baidu Apollo EM Motion Planner.  (2018).  arXiv:1807.08048
  * Fisac et al. (2019) Jaime F Fisac, Eli Bronstein, Elis Stefansson, Dorsa Sadigh, S Shankar Sastry, and Anca D Dragan. 2019.  Hierarchical game-theoretic planning for autonomous vehicles. In _2019 International conference on robotics and automation (ICRA)_. IEEE, 9590–9596.
  * Hang et al. (2022) Peng Hang, Chao Huang, Zhongxu Hu, and Chen Lv. 2022\.  Driving conflict resolution of autonomous vehicles at unsignalized intersections: A differential game approach.  _IEEE/ASME Transactions on Mechatronics_ 27, 6 (2022), 5136–5146.
  * Howard and Kelly (2007) Thomas M. Howard and Alonzo Kelly. 2007.  Optimal Rough Terrain Trajectory Generation for Wheeled Mobile Robots.  _The International Journal of Robotics Research_ 26, 2 (2007), 141–166.  https://doi.org/10.1177/0278364906075328
  * Katrakazas et al. (2015) Christos Katrakazas, Mohammed Quddus, Wen-Hua Chen, and Lipika Deka. 2015.  Real-time motion planning methods for autonomous on-road driving: State-of-the-art and future research directions.  _Transportation Research Part C: Emerging Technologies_ 60 (2015), 416–442.  https://doi.org/10.1016/j.trc.2015.09.011
  * Lenz et al. (2016) David Lenz, Tobias Kessler, and Alois Knoll. 2016\.  Tactical cooperative planning for autonomous highway driving using Monte-Carlo Tree Search.  _IEEE Intelligent Vehicles Symposium, Proceedings_ 2016-Augus, Iv (2016), 447–453.  https://doi.org/10.1109/IVS.2016.7535424
  * Li et al. (2022) Chenran Li, Tu Trinh, Letian Wang, Changliu Liu, Masayoshi Tomizuka, and Wei Zhan. 2022\.  Efficient Game-Theoretic Planning With Prediction Heuristic for Socially-Compliant Autonomous Driving.  _IEEE Robotics and Automation Letters_ 7, 4 (2022), 10248–10255.  https://doi.org/10.1109/LRA.2022.3191241
  * Lopez et al. (2018) Pablo Alvarez Lopez, Michael Behrisch, Laura Bieker-Walz, Jakob Erdmann, Yun-Pang Flötteröd, Robert Hilbrich, Leonhard Lücken, Johannes Rummel, Peter Wagner, and Evamarie Wießner. 2018.  Microscopic Traffic Simulation using SUMO, In The 21st IEEE International Conference on Intelligent Transportation Systems.  _IEEE Intelligent Transportation Systems Conference (ITSC)_.  https://elib.dlr.de/124092/
  * Montemerlo et al. (2008) Michael Montemerlo, Jan Becker, Suhrid Bhat, Hendrik Dahlkamp, Dmitri Dolgov, Scott Ettinger, Dirk Haehnel, Tim Hilden, Gabe Hoffmann, Burkhard Huhnke, et al. 2008\.  Junior: The stanford entry in the urban challenge.  _Journal of field Robotics_ 25, 9 (2008), 569–597.
  * Peng et al. (2021) Zhenghao Peng, Quanyi Li, Ka Ming Hui, Chunxiao Liu, and Bolei Zhou. 2021.  Learning to Simulate Self-Driven Particles System with Coordinated Policy Optimization.  _Advances in Neural Information Processing Systems_ 34 (2021), 10784–10797.
  * Ranney (1994) Thomas A Ranney. 1994\.  Models of driving behavior: a review of their evolution.  _Accident analysis & prevention_ 26, 6 (1994), 733–750.
  * Rovenski (2006) Vladimir Y. Rovenski. 2006\.  _Differential geometry of curves and surfaces: A concise guide_.  Birkhäuser, Boston, MA.  https://doi.org/10.1007/b137116
  * Sadigh et al. (2016) Dorsa Sadigh, Shankar Sastry, Sanjit A Seshia, and Anca D Dragan. 2016.  Planning for autonomous cars that leverage effects on human actions. In _Robotics: Science and systems_ , Vol. 2. Ann Arbor, MI, USA, 1–9.
  * Schwarting et al. (2019) Wilko Schwarting, Alyssa Pierson, Javier Alonso-Mora, Sertac Karaman, and Daniela Rus. 2019\.  Social behavior for autonomous vehicles.  _Proceedings of the National Academy of Sciences of the United States of America_ 116, 50 (2019), 2492–24978.  https://doi.org/10.1073/pnas.1820676116
  * Sun et al. (2022) Qiao Sun, Xin Huang, Junru Gu, Brian C Williams, and Hang Zhao. 2022.  M2I: From factored marginal trajectory prediction to interactive prediction. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 6543–6552.
  * Swaroop and Rajagopal (2001) D. Swaroop and K. R. Rajagopal. 2001.  A review of constant time headway policy for automatic vehicle following. In _IEEE Conference on Intelligent Transportation Systems, Proceedings, ITSC_. IEEE, Oakland, CA, 65–69.  https://doi.org/10.1109/itsc.2001.948631
  * Treiber et al. (2000) Martin Treiber, Ansgar Hennecke, and Dirk Helbing. 2000\.  Congested traffic states in empirical observations and microscopic simulations.  _Physical review E_ 62, 2 (2000), 1805.
  * Urmson et al. (2008) Chris Urmson, Joshua Anhalt, Drew Bagnell, Christopher Baker, Robert Bittner, MN Clark, John Dolan, Dave Duggins, Tugrul Galatali, Chris Geyer, et al. 2008\.  Autonomous driving in urban environments: Boss and the urban challenge.  _Journal of field Robotics_ 25, 8 (2008), 425–466.
  * Vapnik (2006) Vladimir Vapnik. 2006\.  _Estimation of Dependences Based on Empirical Data_.  Springer, New York, NY.  https://doi.org/10.1007/0-387-34239-7
  * Wang et al. (2022b) Junhua Wang, Wenxiang Xu, Ting Fu, Hongren Gong, Qiangqiang Shangguan, and Anae Sobhani. 2022b.  Modeling aggressive driving behavior based on graph construction.  _Transportation Research Part C: Emerging Technologies_ 138 (2022), 103654\.
  * Wang et al. (2022a) Wenshuo Wang, Letian Wang, Chengyuan Zhang, Changliu Liu, and Lijun Sun. 2022a.  Social Interactions for Autonomous Driving: A Review and Perspectives.  _Foundations and Trends in Robotics_ 10, 3-4 (2022), 198–376.  https://doi.org/10.1561/2300000078
  * Wang et al. (2018) Yunpeng Wang, Junjie Zhang, and Guangquan Lu. 2018\.  Influence of driving behaviors on the stability in car following.  _IEEE Transactions on Intelligent Transportation Systems_ 20, 3 (2018), 1081–1098.
  * Werling et al. (2010) Moritz Werling, Julius Ziegler, Sören Kammel, and Sebastian Thrun. 2010.  Optimal trajectory generation for dynamic street scenarios in a frenét frame. In _2010 IEEE International Conference on Robotics and Automation_. IEEE, Anchorage, Alaska, USA, 987–993.  https://doi.org/10.1109/ROBOT.2010.5509799
  * Zheng et al. (2022) Ou Zheng, Mohamed Abdel-Aty, Lishengsa Yue, Amr Abdelraouf, Zijin Wang, and Nada Mahmoud. 2022.  CitySim: A Drone-Based Vehicle Trajectory Dataset for Safety Oriented Research and Digital Twins.  https://doi.org/10.48550/ARXIV.2208.11036


# Spectral Analysis of Proton Eigenfunctions
in Crystalline Environments

L. Gamberale 111Email address<EMAIL_ADDRESS>Quantumatter Inc.,
Dover DE USA
LEDA srl, Università Milano Bicocca
I-20126 Milano, Italy  G. Modanese 222Email address:
<EMAIL_ADDRESS>Free University of Bozen-Bolzano
Faculty of Engineering
I-39100 Bolzano, Italy

###### Abstract

A system with protons confined in a periodic potential is examined using the
Schrödinger equation and the Bloch theorem. Real-world conditions introduce
deviations from an ideal harmonic system due to truncated and non-quadratic
potentials. The energy spectrum of bound eigenfunctions is derived,
independent from the crystal momentum vector $\vec{K}$. Numerical solutions
reveal intriguing structures in the system, including bound states with
degeneracy equal to site number $N_{w}$ resembling a finite harmonic
oscillator spectrum. Notably, the outcomes for protons differ from electronic
energy bands due to proton mass, which results in a greater number of possible
bound states. Building upon previous work, this study investigates the
conditions for electromagnetic coherence onset in the system. The energy gap
and oscillation amplitude constraints are demonstrated. A perturbation term
accounts for deviations in energy level spacing, resulting in minor
electromagnetic mode splitting. Maintaining a phase differential between modes
below $\pi/2$ over an oscillation period is the condition for generating a
negative energy coherent contribution in the interaction with the
electromagnetic field. By analyzing the frequency dependence of the number of
bound states, it is deduced that in a concrete example concerning hydrogen
loading in metals within the octahedral phase the number of bound states is
equal to 12, with an average frequency of 0.4 eV. This resolves the issue of
unlimited energy gap and unbounded oscillation amplitude related to the
coherent transition.

##  I Introduction

The past two decades have witnessed extensive research focused on the
potential to achieve and manipulate coherent and collective quantum optical
phenomena at the nano-scale level, such as superradiance and coherent
population trapping [1, 2, 3, 4, 5, 6, 7, 8].

Superradiance, also known as the Dicke effect, manifests when an ensemble of
molecules or quantum oscillators confined within a sub-wavelength region
collaboratively emit and absorb coherent radiation. Coherence among these
emitters is mediated by the electric radiation field. In some cases, the
field-matter coupling can be significantly enhanced by the presence of surface
plasmons, resulting in what is known as the plasmonic Dicke effect [9, 10, 11,
12, 13, 14, 15, 16, 17].

This study delves into an idealized system comprising $N$ charges oscillating
within a lattice structure, exhibiting a superradiant transition above a
specific threshold density. The dynamics of this system involve bulk plasmons,
a topic we previously analyzed analytically [18]. In our prior work, we
demonstrated the existence of an energy gap for a specific set of
wavefunctions in the large $N$ limit, displaying coherence in both the matter
and field sectors.

In [19], we recall the Hamiltonian for this model and employ canonical
transformations of the photon field operators to diagonalize the photonic
term. That work succinctly summarizes key findings from [18] in a self-
consistent manner, and then describes some numerical simulations showing the
occurrence of a quantum phase transition for small values of $N$, also in the
presence of dissipation.

One of the fundamental characteristics of the model, as extensively discussed
in reference [18], lies in its uniform definition of the proper frequency
denoted as $\omega$ for the oscillators within the material. This uniformity
is achieved through the utilization of a periodic electrostatic potential
referred to as the ”jellium crystal”. Importantly, this potential also serves
the purpose of constraining the amplitude of oscillations within the material.
During the diagonalization process of the total Hamiltonian, the frequency
$\omega$ is combined with the plasma frequency $\omega_{p}$ associated with
the oscillating charges. This combination yields a modified or ”dressed”
frequency denoted as $\omega^{\prime}=\sqrt{\omega^{2}+\omega_{p}^{2}}$. It is
noteworthy that the photon momentum, represented as $\mathbf{k}$, remains
unchanged throughout this process. Consequently, within this model, there
exist no states in which electromagnetic energy generated within the material
can propagate into the vacuum, thereby forming a naturally occurring
resonating cavity.

Another crucial aspect is that the coherent electromagnetic modes that emerge,
oscillating with a fixed phase relation to the matter oscillators, possess
wave vectors pointing in all possible spatial directions. This enhances the
coupling between matter and field, enabling the formation of the energy gap.

In the subsequent sections, we will focus on the bound states of protons bound
to crystal sites. In Sect. II we find the spectrum through a numerical
solution of the Schrödinger equation and taking into account the Bloch
theorem. In Sect. III we consider a perturbation term which accounts for
deviations in energy level spacing, resulting in minor e.m. mode splitting.

Finally, we rigorously calculate the maximum amplitude of the coherent
oscillation, which had been manually set using a heuristic argument in [18],
and summarize our findings.

##  II Single-particle bound states in a finite quadratic periodic potential

In reference [18], we tackled the issue of one-particle bound states within an
ideal periodic infinite potential well. In idealized scenarios, such wells
provide a simplified framework for understanding quantum systems. However,
when considering real-world conditions, oscillators, or particles, deviate
from ideality due to truncated and non-quadratic potentials. This departure
from ideal behavior carries significant consequences for the structure of
bound energy levels. In the subsequent analysis, we will delve into the
profound implications of these deviations, examining how they impact the
characteristics of bound energy states and, more importantly, its consequences
in the properties of their coherent states.

Consider an ensemble of protons with mass $m$ immersed in a periodic harmonic
potential in 3D of the form

$V(\vec{x})=\omega\sum_{ijk}{\mathcal{V}}(\sqrt{m\omega}|\vec{x}-\vec{x}_{ijk}|)$
(1)

where $\omega$ is the free vibration frequency, $\vec{x}_{ijk}$ are the 3D
positions of the lattice sites and

${\mathcal{V}}(q)=\frac{1}{2}(q^{2}-q_{max}^{2})\theta(q_{max}^{2}-q^{2})$ (2)

where $q_{max}$ is a parameter that controls the width and depth of the
potential well. The Heaviside (step) function $\theta$ guarantees that the
potential is always non-positive, as represented in Fig. 1 in the
$x-$direction with $y=z=0$.

The Schrödinger equation that determines the eigenstates of the single-
particle problem is

$\left[-\frac{1}{2m}\nabla^{2}+V(\vec{x})\right]\psi(\vec{x})=E\psi(\vec{x})$
(3)

whose solution can be afforded using the Bloch theorem. By setting

$\psi_{\vec{K},\vec{n}}(\vec{x})=e^{i\vec{K}\cdot\vec{x}}u_{\vec{n}}(\vec{x})$
(4)

and working in the tight binding approximation [20] we can set
$u_{\vec{n}}(\vec{x})=\sum_{ijk}\phi_{\vec{n}}(\vec{x}-x_{ijk})$ where
$\phi_{\vec{n}}(\vec{x})$ are the eigenfunctions of the Schrödinger equation
for the isolated potential $\omega\mathcal{V}(\sqrt{m\omega}|\vec{x}|)$.

The validity of this approximation stems from the observation that the
particles under examination have a mass much larger than that of the
electrons, resulting in their wave-functions displaying significant
localization around the centers of the potential wells. This localization, in
turn, guarantees that the level of overlap between the wave-functions
associated with neighboring sites can be considered practically negligible to
a considerable degree at least up to a certain level of excitation of the
bound states.

The energy spectrum of the bound eigenfunctions described by Equation (3) can,
in fact, be derived from the spectrum of the equation

$\left[-\frac{1}{2m}\nabla^{2}+\omega\mathcal{V}(\sqrt{m\omega}|\vec{x}|)\right]\phi_{\vec{n}}(\vec{x})=E_{\vec{n}}\phi_{\vec{n}}(\vec{x}).$
(5)

independent of the momentum vector $\vec{K}$, that can take any value.

Figure 1:  Periodic 1D potential with 40 sites with lattice spacing
$d=2.5\cdot 10^{-8}\text{ cm}=1.27\cdot 10^{-3}\text{ eV}^{-1}$.  Figure 2:
One particle energy spectrum calculated using the parameters of Table 1. Note
the harmonic-oscillator-like energy spacing of the bound states and their
degeneracy, equal to the number of lattice sites $N_{w}$. The free eigenvalues
have a Brillouin-like structure. Different bans correspond to different
colors. The finite number of wells causes intermediate eigenvalues to appear
between adjacent bands (see sub-figure).

To substantiate the above observations, we now study numerically the 1D
version of the problem, described by the Schröedinger equation

$\left[-\frac{1}{2m}\frac{\partial^{2}}{\partial
x^{2}}+V(x)\right]\psi(x)=E\psi(x)$ (6)

where the 1D potential is given by ($\tilde{x}_{i}=d(i+\frac{1}{2})$,
$N_{w}/2\leq i<N_{w}$, where $N_{w}$ is the number of wells considered in the
calculation)

$V(x)=\omega\sum_{i}\mathcal{V}\left(\sqrt{m\omega}(x-\tilde{x}_{i})\right).$
(7)

Through a rotational symmetry argument, it can be demonstrated that the 1D
problem is analogous to the 3D problem and shares the same energy spectrum.
Specifically, it can be shown that Eq. (6) corresponds to the radial component
of the Schrödinger equation in three dimensions (3D) for a single potential
well (Eq. (5)), focusing on the reduced radial solutions characterized by zero
angular momentum (”$s$-wave” solutions). Consequently, the numerical results
we aim to present remain valid even within the framework of the complete 3D
scenario.

By setting the discretization points $x_{j}=\varepsilon j$, $-(M-1)/2\leq
j\leq(M-1)/2$ ($M$ odd) where the discretization step is
$\varepsilon=\frac{d}{N_{c}}$ where $N_{c}$ is the number of discretization
points per cell, Eq. (6) is rewritten in discretized form as

$\left[-\mathbf{T}+2m\varepsilon^{2}\mathbf{V}\right]\bm{\psi}=2m\varepsilon^{2}E\bm{\psi}$
(8)

where $\mathbf{T}$ is the second derivative operator in discretized form and
is defined by $\mathbf{T}_{i,i\pm 1}=1$, $\mathbf{T}_{i,i}=-2$ for $-M/2\leq
i\leq M/2$, $T_{i,j}=0$ in all other cases,
$\mathbf{V}_{ij}=V(x_{i})\delta_{ij}$ and $\bm{\psi}_{i}=\psi(x_{i})$.

Table 1: Chosen parameters for the numerical calculation parameter | symbol | value
---|---|---
oscillation frequency | $\omega^{\prime}$ | 0.41 eV
proton mass | $m$ | 938 MeV
lattice spacing | $d$ | 2.5 Angstrom
number of points of a single cell | $N_{c}$ | 257
number of points of the potential wells | $N_{w}$ | 101
spatial discretization | $\varepsilon$ | 0.0097 Angstrom
number of crystal sites | $N_{s}$ | 41
number of points inside the crystal | $M_{c}$ | 10537
total number of points | $M$ | 12537
Figure 3:  Eigenfunctions of the bound states computed by direct
diagonalization of Eq. (6) with the parameters of Table 1.  Figure 4:  Linear
fit of the photon energies associated to the dipolar transitions between
adjacent bound energy levels (the last two points to the right are excluded
from the fit).

The final parameter requiring specification is the selection of the value for
$q_{max}$. In the context of addressing the challenge of hydrogen loading in
metals within the octahedral phase, prior literature indicates that the
available radius within the octahedral voids is approximately 0.41 $d/2$.
Thus, our decision is to set $q_{max}=0.41\sqrt{m\omega}\frac{d}{2}$.

The eigenvalue problem of Eq. (8) is now completely set and can be solved
numerically. The numerical solution with the parameters indicated in Tab. 1 is
shown in Fig. 2 and reveals a very interesting structure. The bound states
with negative energy have a degeneracy equal to the site number $N_{w}$ as is
expected for the tight binding case and their spectrum is similar to that of
the harmonic oscillator. The free states (positive energy) follow a typical
Brillouin dispersion relation, since the degeneracy is removed by the fact
that the tight binding approximation breaks down, the wave-functions not being
confined in their respective lattice sites. The eigenfunctions of the bound
states are very similar to the eigenfunctions of a harmonic oscillator
centered in a specific site of the lattice (see Fig. 3).

One might inquire about the difference between the well-established outcomes
regarding electronic energy bands, which exhibit a significantly distinct
structure compared to the solutions obtained here for protons. The answer lies
in the fact that, given that protons are approximately 1800 times more massive
than electrons, the lower-energy wave functions are considerably more
localized than those of electrons. This localization arises from the spatial
scaling factor being $\frac{1}{\sqrt{m\omega}}$, consequently leading to a
significantly greater number of possible bound states. In fact, the same
calculation conducted for electrons reveals that the maximum number of bound
states is at most 2.

A closer inspection reveals that the negative energy levels are not exactly
equally spaced, resulting in several frequencies of the electromagnetic field
being involved. Fig. 4 shows the dependence of the photon energies for the
various couplings between adjacent energy levels. Such an issue will be dealt
with in the next section.

##  III Quasi-harmonic case

In [18] we have analyzed the conditions for the onset of electromagnetic
coherence for a system of $N$ charged oscillators immersed in a periodic
harmonic potential and we have heuristically inferred the conditions for the
non-divergence of the energy gap and oscillation amplitude. More specifically,
the divergence of the cited quantities depends directly on the unlimited
number of the bounded energy levels of the harmonic oscillator.

We now demonstrate that when considering the more realistic scenario in which
the periodic potential is confined as outlined in Section II, both the energy
gap and the oscillation amplitude are naturally constrained without any
additional assumption.

Another complicating factor in our analysis arises from the non-uniform
spacing of energy levels within the oscillators, leading to the coupling of
several electromagnetic field modes with the oscillators, as opposed to just a
single mode.

To address this complexity, we begin by examining the idealized case in which
all oscillator levels, denoted as $E_{i}=(i+\frac{1}{2})\omega_{0}$, are
uniformly spaced and we introduce a small perturbation term of
$\frac{1}{2}\delta_{\omega}i(i-1)$ to account for deviations. This
modification results in an energy difference between adjacent levels of
$E_{i+1}-E_{i}=\omega_{0}+\delta_{\omega}i$.

Commencing with $\delta_{\omega}=0$, we initially observe that the coherent
state comprises a solitary frequency, denoted as $\omega^{\prime}$.
Subsequently, we incrementally increase $\delta_{\omega}$ by an infinitesimal
amount, resulting in an exceedingly minor splitting of the various
electromagnetic modes. For the sake of continuity, we anticipate that the
state characterized by the minimum energy still consists of a solitary
electromagnetic mode.

Upon conducting a thorough examination of the phase transition mechanism, we
ascertain that, in order to generate a negative contribution of the
interaction term throughout a single oscillation period of the oscillators, it
is imperative that the phase differential between the modes with frequencies
$\omega_{i}$ originating from the $i$-th energy level and the frequency
$\omega^{\prime}$ associated with the coherent mode remains confined to a
value less than $\frac{\pi}{2}$.

The above condition can be implemented by requiring

$|\omega_{i}-\bar{\omega}|<\frac{\bar{\omega}}{4}$ (9)

where $\bar{\omega}$ is the averaged oscillation frequency of the em field and

$\omega_{i}=\omega_{0}-|\delta_{\omega}|i\ \ \ \ \ \ 0\leq i\leq i_{max}$ (10)

and where $\delta_{\omega}$ is the coefficient of variation of $\omega_{i}$
with the level index (see Fig. 4). Moreover, we choose the average frequency
$\bar{\omega}$ as the mean position with respect to the involved frequencies
so that

$\bar{\omega}=\omega_{0}-|\delta_{\omega}|\frac{i_{max}}{2}.$ (11)

Solving Eq. (9) for the variable $i$, we obtain the maximum value as
$i_{max}=\frac{2\omega_{0}}{5|\delta_{\omega}|}$, and the average frequency is
$\bar{\omega}=\frac{4\omega_{0}}{5}$. Plugging in the numerical values from
Table 1 and referring to Fig. 4, we readily deduce that $i_{\text{max}}\approx
200$. Utilizing the specified parameters in Table 1, we determine that the
number of bound states is $N_{b}=i_{\text{max}}=12$ (see Fig. 2).
Consequently, we can infer that all these bound states contribute to the
establishment of coherence, resulting in an averaged frequency of
$\bar{\omega}=0.4$ eV.

Finally, we can compute the maximum amplitude of the coherent oscillation that
is the parameter that ultimately controls the size of the energy gap (see Eqs.
(30) and (43b) of [18] for the up-to-1st and up-to-2nd order calculation). An
estimate of the squared amplitude of the coherent state $\ket{\alpha}$ is the
number $N_{exc}$ of the highest excitation contributing to the harmonic
oscillator that, in our case, is $N_{b}$. We can roughly say that
$|\alpha|^{2}\simeq N_{b}$. By making use of Eq. (34) of [18] and, being for
protons $\omega^{\prime}=0.41$ eV (Eq.(A24) of [18]), we find an estimate for
the oscillation amplitude

$f=\frac{2}{d}\sqrt{\frac{N_{b}}{m\omega^{\prime}}}\simeq 0.28,$ (12)

to be compared with the guessed value of 0.4 in Ref. [18], found on purely
heuristic grounds. Also the value of the energy gap changes from the estimated
value of 1 eV to the more accurate value 0.46 eV.

To summarize, the limited number of energy levels involved in the coherent
transition immediately resolves the problem related to the unlimited energy
gap foreseen in Ref. [18] and the related unboundedness of the oscillation
amplitude. In Ref. [18] we imposed ’by hand’ the limitedness of the coherent
oscillation through a heuristic physical argument. In the present paper we
theoretically substantiate such an ansatz through an explicit calculation and
we find a result essentially consistent with the previous result.

## References

  * [1] T. Brandes, “Coherent and collective quantum optical effects in mesoscopic systems,” Physics Reports, vol. 408, no. 5-6, pp. 315–474, 2005.
  * [2] M. Heyl, “Dynamical quantum phase transitions: a review,” Reports on Progress in Physics, vol. 81, no. 5, p. 054001, 2018.
  * [3] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, “Ultrastrong coupling regimes of light-matter interaction,” Reviews of Modern Physics, vol. 91, no. 2, p. 025005, 2019.
  * [4] V. Link and W. T. Strunz, “Dynamical phase transitions in dissipative quantum dynamics with quantum optical realization,” Physical Review Letters, vol. 125, no. 14, p. 143602, 2020.
  * [5] S. I. Azzam, A. V. Kildishev, R.-M. Ma, C.-Z. Ning, R. Oulton, V. M. Shalaev, M. I. Stockman, J.-L. Xu, and X. Zhang, “Ten years of spasers and plasmonic nanolasers,” Light: Science & Applications, vol. 9, no. 1, p. 90, 2020.
  * [6] V. Bordo, “Cooperative effects in spherical spasers: Ab initio analytical model,” Physical Review B, vol. 95, no. 23, p. 235412, 2017.
  * [7] D. J. Bergman and M. I. Stockman, “Surface plasmon amplification by stimulated emission of radiation: quantum generation of coherent surface plasmons in nanosystems,” Physical Review Letters, vol. 90, no. 2, p. 027402, 2003.
  * [8] M. I. Stockman, “The spaser as a nanoscale quantum generator and ultrafast amplifier,” Journal of Optics, vol. 12, no. 2, p. 024004, 2010.
  * [9] R. H. Dicke, “Coherence in spontaneous radiation processes,” Physical Review, vol. 93, no. 1, p. 99, 1954.
  * [10] A. V. Andreev, V. I. Emel’yanov, and Y. A. Il’inskiĭ, “Collective spontaneous emission (Dicke superradiance),” Soviet Physics Uspekhi, vol. 23, no. 8, p. 493, 1980.
  * [11] S. Sivasubramanian, A. Widom, and Y. Srivastava, “Gauge invariant formulations of Dicke-Preparata super-radiant models,” Physica A: Statistical Mechanics and its Applications, vol. 301, no. 1-4, pp. 241–254, 2001.
  * [12] S. Sivasubramanian, A. Widom, and Y. Srivastava, “Super-radiance and the unstable photon oscillator,” International Journal of Modern Physics B, vol. 15, no. 05, pp. 537–548, 2001.
  * [13] V. N. Pustovit and T. V. Shahbazyan, “Cooperative emission of light by an ensemble of dipoles near a metal nanoparticle: The plasmonic Dicke effect,” Physical Review Letters, vol. 102, no. 7, p. 077401, 2009.
  * [14] V. N. Pustovit and T. V. Shahbazyan, “Plasmon-mediated superradiance near metal nanostructures,” Physical Review B, vol. 82, no. 7, p. 075429, 2010.
  * [15] M. Scheibner, T. Schmidt, L. Worschech, A. Forchel, G. Bacher, T. Passow, and D. Hommel, “Superradiance of quantum dots,” Nature Physics, vol. 3, no. 2, pp. 106–110, 2007.
  * [16] V. V. Temnov and U. Woggon, “Superradiance and subradiance in an inhomogeneously broadened ensemble of two-level systems coupled to a low-q cavity,” Physical Review Letters, vol. 95, no. 24, p. 243602, 2005.
  * [17] J. A. Greenberg and D. J. Gauthier, “Steady-state, cavityless, multimode superradiance in a cold vapor,” Physical Review A, vol. 86, no. 1, p. 013823, 2012.
  * [18] L. Gamberale and G. Modanese, “Coherent plasma in a lattice,” Symmetry, vol. 15, no. 2, p. 454, 2023.
  * [19] L. Gamberale and G. Modanese, “Numerical simulations unveil superradiant coherence in a lattice of charged quantum oscillators,” Physica B: Condensed Matter, vol. 671, p. 415406, 2023.
  * [20] J. C. Slater and G. F. Koster, “Simplified LCAO method for the periodic potential problem,” Physical review, vol. 94, no. 6, p. 1498, 1954.


# Direct observation of a Fermi surface and superconducting gap in LuNi2B2C

P. Starowicz  Ames Laboratory and Department of Physics and Astronomy, Iowa
State University, Ames, IA 50011, USA  M. Smoluchowski Institute of Physics,
Jagiellonian University, Reymonta 4, 30-059 Kraków, Poland  C. Liu  Ames
Laboratory and Department of Physics and Astronomy, Iowa State University,
Ames, IA 50011, USA  R. Khasanov  Ames Laboratory and Department of Physics
and Astronomy, Iowa State University, Ames, IA 50011, USA  Physik-Institut der
Universität Zürich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland  T.
Kondo  Ames Laboratory and Department of Physics and Astronomy, Iowa State
University, Ames, IA 50011, USA  G. Samolyuk  Ames Laboratory and Department
of Physics and Astronomy, Iowa State University, Ames, IA 50011, USA  D.
Gardenghi  Ames Laboratory and Department of Physics and Astronomy, Iowa State
University, Ames, IA 50011, USA  Bob Jones University, Greenville, SC 29614,
USA  Y. Lee  Ames Laboratory and Department of Physics and Astronomy, Iowa
State University, Ames, IA 50011, USA  T. Ohta  Advanced Light Source,
Berkeley National Laboratory, Berkeley, CA 94720, USA  B. Harmon  Ames
Laboratory and Department of Physics and Astronomy, Iowa State University,
Ames, IA 50011, USA  P. Canfield  Ames Laboratory and Department of Physics
and Astronomy, Iowa State University, Ames, IA 50011, USA  S. Bud’ko  Ames
Laboratory and Department of Physics and Astronomy, Iowa State University,
Ames, IA 50011, USA  E. Rotenberg  Advanced Light Source, Berkeley National
Laboratory, Berkeley, CA 94720, USA  A. Kaminski  Ames Laboratory and
Department of Physics and Astronomy, Iowa State University, Ames, IA 50011,
USA

###### Abstract

We measured the Fermi surface (FS), band dispersion and superconducting gap in
LuNi2B2C using Angle Resolved Photoemission Spectroscopy. Experimental data
were compared with the tight-binding version of the Linear Muffin-Tin Orbital
(LMTO) method and Linearized Augmented Plane-Wave (LAPW) calculations. We
found reasonable agreement between the two calculations and experimental data.
The measured FS exhibits large parallel regions with a nesting vector that
agrees with a previous positron annihilation study and calculations of the
generalized susceptibility. The measured dispersion curves also agree
reasonably well with the TB-LMTO calculations, albeit with some differences in
the strength of the hybridization. In addition, the spectrum in the
superconducting state revealed a 2meV superconducting gap. The data also
clearly shows the presence of a coherent peak above the chemical potential,
$\mu$ that originates from thermally excited electrons above the energy of
2$\Delta$. This feature was not previously observed in the Lu-based material.

###### pacs:

74.70.Dd, 71.18.+y, 71.20.-b, 71.27.+a

##  I Introduction

Figure 1: (Color online) Comparison between the Fermi surface maps measured by
ARPES and the linear muffin-tin orbital (TB-LMTO) calculation. (a) Sketch of
the first Brillouin zone for LuNi2B2C. (b-d) ARPES mapping at the chemical
potential for incident photon energies of 128.13 eV, 119.44 eV, and 102.98 eV,
respectively. (e-g) The Fermi surface maps obtained by TB-LMTO calculations
for constant $k_{z}$ values equal to (e) 0.2, (f) 0.15 and (g) 0.8 expressed
in the units of $\Gamma$ \- Z distance.

Rare earth nickel borocarbides RNi2B2C (R - rare earth) constitute an
interesting class of materials Cava1994 ; Nagarajan1994 ; Canfield1998 ;
MullerNarozhnyi2001 ; MazumdarNagarajan2005 , in which there is a competition
and coexistence between superconductivity and magnetism. Amongst these
compounds, nonmagnetic LuNi2B2C has the highest superconducting critical
temperature of 16.6 K Cava1994 . The borocarbides exhibit a peculiar
anisotropy of the superconducting gap, the character of which is still under
debate. It is believed that the gap is highly anisotropic in the two non-
magnetic compounds LuNi2B2C and YNi2B2C Boaknin2001 ; Bobrov2005 ; Maki2002 ;
Raychaudhuri2004 ; MartinezSamper2003 ; Izawa2002 ; Yokoya2000 . Its symmetry
was proposed to be s \+ g Maki2002 , which is consistent with certain
experimental results Raychaudhuri2004  but an anisotropic s-wave symmetry has
also been considered MartinezSamper2003 . Other experimental data indicate
that the gap in YNi2B2C has point nodes along the (100) and (010) directions
Izawa2002 . LuNi2B2C crystallizes in a body-centered tetragonal structure with
lattice parameters a = 3.4639 Å, and c = 10.6313 ÅSiegrist1994 . Its crystal
structure consists of Lu-C layers with Ni2B2 sheets in between. Previously
calculations reveal that LuNi2B2C is characterised by a large density of
states (DOS) at the Fermi energy ($E_{F}$) originating mainly, but not
exclusively, from Ni d electrons Mattheiss1994 ; PickettSingh1994 ;
Coehoorn1994 . Another interesting feature is a flat band along the $\Gamma$-X
direction just above $E_{F}$. The Fermi surface (FS) topography of LuNi2B2C
was studied by ab-initio calculations Rhee1995 ; Kim1995 ; Dugdale1999 . Band
structure calculations Rhee1995  revealed a pronounced maximum in the
generalized electronic susceptibility at ($\sim$0.6a*, 0, 0), where a* $\equiv
2\pi/a$ and most likely arising from large nested regions of the FS. Moreover,
phonon softening was observed in LuNi2B2C by means of inelastic neutron
scattering for a range of wave vectors around (0.5a*, 0, 0) Dervenagas1995 .
Interestingly enough the magnetic ordering, which was found in RNi2B2C
compounds with magnetic atoms R = Er, Ho, Tb and Gd manifest similar
modulation vector usually close to (0.55a*,0,0)MullerNarozhnyi2001 . The first
experimental studies of the LuNi2B2C (RNi2B2C) Fermi surface were performed by
means of two-dimensional angular correlation of electron-positron annihilation
radiation (2D-ACAR) and the data were compared to the Linear Muffin-Tin
Orbital (LMTO), local density approximation (LDA) calculations Dugdale1999 .
Nested parts of the FS were found with a nesting vector corresponding to both
the phonon softening and the magnetic modulation vectors. The fraction of the
FS participating in nesting was determined to be 4.4 $\pm$ 0.5% Dugdale1999 .
That study was however limited only to a rough “calipering” of the Fermi
surface. Knowledge of the experimental band structure, Fermi surface and
quasiparticle properties is deemed essential to understand the interplay of
the various interactions in these materials, as it may shed new light on other
phenomena such as anisotropic superconductivity, the role of phonon softening
and the relationship between the superconductivity and magnetic ordering in
the borocarbides. It is also a pre-requisit for direct determination of the
alleged anisotropy of the superconducting gap in these materials. In this
report we present angle resolved photoelectron spectroscopy (ARPES)
measurements of the band dispersion, Fermi surface and supercondcuting gap in
the borocarbide with the highest $T_{c}$, LuNi2B2C. The experimental results
were compared with the tight-bounding LMTO (TB-LMTO) method and the full
potential LAPW calculations. We found reasonable agreement with theory. The
most significant difference between the calculations and experimental data is
the strength of the hybridization. We also determined the superconducting gap
to be 2.58 meV (extrapolated for T=0), in good agreement with the gap expected
from the superconducting transition temperature
(${{2\Delta}/{k_{B}T_{C}}}=2.78$).

##  II Experimental

LuNi2B2C single crystals were grown at Ames Laboratory by means of a high-
temperature flux techniqueCanfield1998 ; Fisher2001 . The plate-like crystals
were cleaved in situ at pressures better than 3 x 10-11 Tr to reveal and
maintain fresh a-b surfaces. The Fermi surface and band structure mapping were
performed at the 7.0.1 beamline at the Advanced Light Source, using a Scienta
R4000 analyzer. The energy and angle resolution were set at $\sim$ 30 meV and
$\sim$ 0.5 deg, respectively. The energy gap was measured with a Scienta 2002
analyser and He-I photon source (h$\nu$ = 21.2 eV), in which the overall
energy resolution was set at 2 meV. The normal state data were measured at the
Synchrotron Radiation Center using the PGM beamline and Scienta 2002
endstation, with the energy and angular resolution set at $\sim$ 13 meV and
0.25 deg, respectively. Tight-binding linear muffin-tin orbital calculations
were performed by the TB-LMTO program, version 47 LMTO , and the Full-
Potential Linearized Augmented Plane-Wave (LAPW) calculations were performed
using the Wien2k package Wien .

Figure 2: Dispersion of the conduction bands obtained with a photon energy of
102.98 eV (same as Fig. 1d), compared with the TB-LMTO calculation. a) Fermi
surface map with momentum cuts indicated by the solid lines, in which cut (c)
and (i) pass through the $\Gamma$ and X points, respectively. b) Fermi surface
contours obtained by the TB-LMTO calculation for the value of $k_{z}$
corresponding to the data in panel (a). Panels (c-i): measured band dispersion
along the cuts indicated in (a). Panels (j-p): calculated band dispersion
along the cuts marked in panel (b).

##  III Results and discussion

Band structure and semi-planar Fermi Surface cuts were determined for incident
photon energies 128.13 eV, 119.44 eV and 102.98 eV (Fig. 1b-d), where the
$\Gamma$ point in the Brillouin zone (Fig. 1a) corresponds to normal emission
of electrons along the (001) direction. A $\sin r/r$ correction term ($r$ is
the distance from the $\Gamma$ point) was used to account for mapping of the
momentum space onto the angular distribution of photoelectrons. The ARPES
process in 3D materials leaves some ambiguity as to the $k_{z}$ component of
the momentum (perpendicular to the sample surface), because it is not
conserved in the photoemission process. This is due to jump of the potential
at the sample surface. From conservation of energy and remaining components of
the momentum one can calculate the relative changes of the $k_{z}$ for various
photon energies. To obtain the offset one needs to seek guidance from the band
structure calculations and identify the high symmetry points in the
data.HUFNER  This allows estimation of the $k_{z}$ offset. The change of the
wave vector component $k_{z}$ (parallel to the c axis) between scans in Fig.
1b ($h\nu$=128.13) and 1c ($h\nu$=119.44) was calculated from momentum and
energy conservation to be 0.33 of the $\Gamma$-Z distance. Similarly the
change of the wave vector component $k_{z}$ between scans in Fig. 1c
($h\nu$=119.44) and 1d ($h\nu$=102.98) was 0.66 of the $\Gamma$-Z. The
calculated Fermi surfaces were obtained for constant $k_{z}$ values by means
of the TB-LMTO method and are shown in Fig. 1 panels e-g. We estimated the
values of the inner potential, $V_{0}$ = 9.4 eV and the work function $\phi$ =
4.6 eV, by comparing the high symmetry points between the calculated Fermi
surfaces and the experimental data. This allowed us to determine the offset of
the photon energy that corresponds to $k_{z}$ = 0.

The Fermi surface maps for the incident photon energies of 128.13 eV and
119.44 eV reveal large parallel parts of the FS with essentially the same
nesting vector (spacing between the linear sections): $k_{n}$ = (0.59 $\pm$
0.04)a* for Fig. 1b and $k_{n}$ = (0.58 $\pm$ 0.04)a* for Fig. 1c. $k_{z}$ is
expressed in the units of the $\Gamma$-Z distance, where the $\Gamma$ point
corresponds to $k_{z}$ = 0. Although the full three dimensional FS was not
determined in great detail, a similar nesting vector was found for different
$k_{z}$ values which indicates that the FS likely has considerable nesting
properties for a wide range of $k_{z}$ values. The constancy of the value of
the nesting vector between $k_{z}$=0.15 and -0.2 is also consistent with
results of calculations. The spacing between the parallel segments of the
Fermi surface predicted by TB-LMTO calculation is between 0.54a* and 0.55a*,
the LAPW calculation results (not shown) are 0.50a* and 0.57a*, respectively.
The detected $k_{n}$ is very close to the theoretically predicted value
obtained from the generalized susceptibility Rhee1995 . Our results also agree
reasonably well with the nesting vector previously determined via 2D-ACAR
Dugdale1999 .

The Fermi surface map obtained at 102.98 eV very closely resembles the
calculated Fermi surface for $k_{z}$ = 0.8a*. The overall shapes of the
measured and calculated Fermi surface sheets (Fig. 1d and 1g) are very
similar, however there is one significant difference. In the calculations the
four oval parts of the Fermi surface centered about $\Gamma$-Z are well
separated in momentum space (Fig. 1g), while the data reveals that they
actually are connected at the edges (Fig. 1d). A lack of separation in the
experimental data may indicate that the hybridization gap is overestimated in
the calculations. These oval parts arise from the intersection of the electron
and hole-like bands. Interestingly enough at the edges along the diagonal
directions (e. g. 110) the bottom of the electron band and the top of the hole
band appear to be pinned at the chemical potential, resulting in a
characteristic “flower” shape.

In Fig. 2 we plot the band dispersion data along a few selected cuts in
momentum space obtained at an incident photon energy of 102.98 eV (Figs.
2c-i), along with a calculated (TB-LMTO method) band dispersion for
$k_{z}\sim$ 0.8a* (Figs. 2j-p). The agreement between the measured and
calculated band dispersion is rather good, especially in the proximity of the
chemical potential. In the corresponding TB-LMTO calculations (Figs. 2j-p),
the same overall features are well reproduced, which shows the validity of the
calculation in this material to a certain extent. This agreement also
validates the assignment of $k_{z}$ values to the cuts measured at various
photon energies, which is very important when studying 3D materials with
ARPES. The most significant difference is the hybridization gap, which is
quite large in the calculations but its signatures are for the most part
absent in the measured data. For example in Fig. 2 panels (k) and (l) the high
and low energy branches form hybridization gap of about 200 meV at E = -0.3
eV, while in the corresponding measured data (panels (d) and (e)) the bands
appear to disperse without a signature of the hybridization gap. One should
consider if the observed features have any relation to the superconducting gap
asymmetry and observed phonon softening in LuNi2B2C. It should be noted that
band structure calculations PickettSingh1994 ; Kim1995  show a flat band lies
very close to, but slightly above, the Fermi level. This feature was
unfortunately not observed in our data due to the Fermi function cut-off.
However, a higher DOS near the Fermi level would explain the large number of
scattered electrons observed with k vectors along (110) and phonon softening
for the discussed wave vectors. Consequently this may lead to an anisotropy of
the superconducting order parameter. This is in agreement with the results
proposing for YNi2B2C that the superconducting gap is larger just at (110) and
diminishes or even has nodes along the (100) and (010) directions Izawa2002 .

Given the above concern, we measured the energy gap in LuNi2B2C by partial
angle-integrated photoelectron spectroscopy and compared with the normal state
Fermi surface. The opening of the superconducting gap is clearly shown in Fig.
3. In order to determine the magnitude of the gap, the Dynes function
Dynes1978  was fitted to the symmetrized Norman1998  spectrum (Fig. 3b). The
fitted function yields the gap value of $\Delta$ = 1.5 meV for the sample at T
= 11 K with the $\Gamma$ parameter equal to 0.05 meV.

The striking feature in Fig. 3c is the pronounced peak above the chemical
potential. This peak arises from thermal excitation of electrons above the
2$\Delta$. This points to high DOS just above $\mu$ which is consistent with
the idea that the flat band along (110) direction, a large part of which is
slightly over the Fermi energy playing an important role in this anisotropic
superconductivity. Similar peaks were recently reported in Y based
borocarbides YOKOYA2007 . According to BCS theory the energy gap value at zero
temperature ($\Delta_{0}$) is 2.45 meV for $T_{c}$ = 16 K superconductor. Our
$\Delta$ value of 1.5 meV obtained at T=11K corresponds to ($\Delta_{0}$)=2.6
meV, in excellent agreement with the BCS predictions.

Figure 3: (Color online) Superconducting gap of LuNi2B2C a) measured at T =
11$\pm$1 K, compared with the normal state at T = 40 K. b) The Dynes function
(solid red line) with the parameters $\Delta$ = 1.5 meV and $\Gamma$ = 0.05
meV fitted to the symmetrised spectrum (solid black circles). c) enlarged
portion of superconducting spectra from (a) close to the chemical potential.

##  IV Conclusions

We have performed measurements of the Fermi surface, band dispersion and
supercondcuting gap for highest $T_{c}$ rare earth nickel borocarbide
superconductor LuNi2B2C. The experimental data were compared with two
different density functional calculations. The overall agreement between
theory and measurement is good. In the experiment, large parallel FS parts
spaced with the vector $k_{n}$ = 0.59a* have been found for two different
incident photon energies, which is a confirmation of the previous theoretical
predictions Rhee1995  and earlier experimental studies Dugdale1999 . The
calculated FS confirms the existence of large nested parts, with a nesting
vector in good agreement with the ARPES results presented here. The
superconducting gap was measured and we also observed a coherent peak above
the chemical potential. This peak arises due to electrons being thermally
excited above the energy of 2$\Delta$.

##  V Acknowledgments

This work was supported by Director Office for Basic Energy Sciences, US DOE.
Work at Ames Laboratory was supported by the Department of Energy - Basic
Energy Sciences under Contract No. DE-AC02-07CH11358. Advanced Light Source is
operated by the U.S. DOE under Contract No. DE-AC03-76SF00098. Synchrotron
Radiation Center is supported by the National Science Foundation under award
No. DMR-0537588. R. K. gratefully acknowledges support of K. Alex Müller
Foundation.

## References

  * (1) R. J. Cava, H. Takagi, H. W. Zandbergen, J. J. Krajewski, W. F. Peck Jr, T. Siegrist, B. Batlogg, R. B. van Dover, R. J. Felder, K. Mizuhashi, J. O. Lee, H. Eisaki, and S. Uchida, Nature 367, 252 (1994).
  * (2) R. Nagarajan, C. Mazumdar, Z. Hossain, S.K. Dhar, K. V. Gopalakrishnan, L. C. Gupta, C. Godart, B. D. Padalia, and R. Vijayaraghavan, Phys. Rev. Lett. 72, 274 (1994).
  * (3) P. C. Canfield, P. L. Gammel, and D. J. Bishop, Physics Today 51, 40 (1998).
  * (4) K.-H. Müller and V. N. Narozhnyi, Rep. Prog. Phys. 64, 943 (2001).
  * (5) C. Mazumdar and R. Nagarajan, Current Science 88, 83 (2005).
  * (6) T. Yokoya, T. Kiss, T. Watanabe, S. Shin, M. Nohara, H. Takagi, and T. Oguchi, Phys. Rev. Lett. 85, 4952 (2000).
  * (7) E. Boaknin, R. W. Hill, C. Proust, C. Lupien, and L. Taillefer, Phys. Rev. Lett. 87, 237001 (2001).
  * (8) N. L. Bobrov, S. I. Beloborod’ko, L. V. Tyutrina, I. K. Yanson, D. G. Naugle, and K. D. D. Rathnayaka, Phys. Rev. B 71, 014512 (2005).
  * (9) K. Maki, P. Thalmeier, and H. Won, Phys. Rev. B 65, 140502 (2002).
  * (10) P. Raychaudhuri, D. Jaiswal-Nagar, Goutam Sheet, S. Ramakrishnan, and H. Takeya, Phys. Rev. Lett. 93, 156802 (2004).
  * (11) P. Martinez-Samper, H. Suderow, S. Vieira, J. P. Brison, N. Luchier, P. Lejay, and P. C. Canfield, Phys. Rev. B 67, 014526 (2003).
  * (12) K. Izawa, K. Kamata, Y. Nakajima, Y. Matsuda, T. Watanabe, M. Nohara, H. Takagi, P. Thalmeier, and K. Maki, Phys. Rev. Lett. 89, 137006 (2002).
  * (13) T. Siegrist, H. W. Zandbergen, R. J. Cava, J. J. Krajewski, and W. F. Peck Jr, Nature 367, 254 (1994).
  * (14) L. F. Mattheiss, Phys. Rev. B 49, 13279 (1994).
  * (15) W. E. Pickett and D. J. Singh, Phys. Rev. Lett. 72 , 3702 (1994).
  * (16) R. Coehoorn, Physica C 228, 5671 (1994).
  * (17) J. Y. Rhee, X. Wang, and B. N. Harmon, Phys. Rev. B 51, 15585 (1995).
  * (18) H. Kim, C.-D. Hwang, and J. Ihm, Phys. Rev. B 52, 4592 (1995).
  * (19) S. B. Dugdale, M. A. Alam, I. Wilkinson, R. J. Hughes, I. R. Fisher, P. C. Canfield, T. Jarlborg, and G. Santi, Phys. Rev. Lett. 83, 4824 (1999).
  * (20) P. Dervenagas, M. Bullock, J. Zarestky, P. Canfield, B. K. Cho, B. Harmon, A. I. Goldman, and C. Stassis, Phys. Rev. B 52, 9839 (1995).
  * (21) P. C. Canfield, I. R. Fisher, J. Crystal Growth, 225, 155-161 (2001) .
  * (22) O. Jepsen and O. K. Anderson, Solid State Commun. 9, 1763 (1971) and Phys. Rev. B 29, 5965 (1984); P. Blöchl, O. Jepsen and O. K. Anderson, Phys. Rev. B 49, 16223 (1994).
  * (23) P. Blaha, K. Schwarz, G. Madsen, D. Kvasnicka and J. Luitz, (2001) _WIEN2k, An Augmented Plane Wave_ \+ _Local Orbitals Program for Calculating Crystal Properties_ (Karlheinz Schwarz, Tech. Univ. Wien, Austria).
  * (24) S. Hufner, Photoelectron Spectroscopy, pp 268-270 (Springer, Berlin 1995).
  * (25) R. C. Dynes, V. Narayanamurti, and J. P. Garno, Phys. Rev. Lett. 41, 1509 (1978).
  * (26) M. R. Norman, H. Ding, M. Randeria, J. C. Campuzano, T. Yokoya, T. Takeuchi, T. Takahashi, T. Mochiku, K. Kadowaki, P. Guptasarma, and D. G. Hinks, Nature 392, 157 (1998).
  * (27) T. Baba, T. Yokoya, S. Tsuda, T. Kiss, T. Shimojima, K. Ishizaka, H. Takeya, K. Hirata, T. Watanabe, M. Nohara, H. Takagi, N. Nakai, K. Machida, T. Togashi, S. Watanabe, X.-Y. Wang, C. T. Chen, and S. Shin, Phys. Rev. Lett. 100, 017003 (2008).