# Inclined Quadrotor Landing using Deep Reinforcement Learning∗

Jacob E. Kooi1 and Robert Babuška2 1Jacob E. Kooi is with the Departments of
Cognitive Robotics and Delft Center for Systems and Control, Delft University
of Technology, 2628 CD Delft, The Netherlands<EMAIL_ADDRESS>Babuška with the Department of Cognitive Robotics, Delft University of
Technology, 2628 CD Delft, The Netherlands<EMAIL_ADDRESS>This work has
been submitted to the IEEE for possible publication. Copyright may be
transferred without notice, after which this version may no longer be
accessible.

###### Abstract

Landing a quadrotor on an inclined surface is a challenging manoeuvre. The
final state of any inclined landing trajectory is not an equilibrium, which
precludes the use of most conventional control methods. We propose a deep
reinforcement learning approach to design an autonomous landing controller for
inclined surfaces. Using the proximal policy optimization (PPO) algorithm with
sparse rewards and a tailored curriculum learning approach, a robust policy
can be trained in simulation in less than 90 minutes on a standard laptop. The
policy then directly runs on a real Crazyflie 2.1 quadrotor and successfully
performs real inclined landings in a flying arena. A single policy evaluation
takes approximately 2.5 ms, which makes it suitable for a future embedded
implementation on the quadrotor.

##  I Introduction

Modern quadrotors are agile and can perform complex tasks in difficult-to-
reach places. Quadrotor flight and manoeuvres are commonly controlled by
proportional integral derivative (PID) control or model predictive control
(MPC). Although these methods are adequate for setpoint or trajectory
tracking, they fall short when it comes to more complicated manoeuvres that
exceed the linearization range or require long prediction horizons. One of
such manoeuvres is the landing on an inclined surface, which is relevant for
applications like delivery, maintenance or surveillance. In order to
facilitate a safe inclined landing, the final attitude of the quadrotor must
match the slope of the landing platform. The final state of the landing
trajectory is not an equilibrium, which presents a challenge for the control
design. Owing to the under-actuated nature of the system, the landing
trajectory can be long and complex, with an initial motion away from the
landing location. This complicates the use of standard control methods like
MPC with a fixed prediction horizon and quadratic cost function.

Recent advances in deep reinforcement learning (DRL) with continuous action
spaces have made this approach suitable also for quadrotor control [1, 2, 3,
4], including landing controllers [5, 6, 7, 8, 9, 10]. However, no results
have yet been reported for inclined landing. In this paper, we develop a DRL
approach to the inclined landing problem and validate it in simulations and
real lab experiments with the Crazyflie 2.1 Nano-UAV. To the best of our
knowledge, this is the first deep-learning-based controller for inclined
landing applied to a real quadcopter. More specifically, our contributions
are:

  * •

We develop two fast Gym-based [11] simulation environments for the Crazyflie
2.1 Nano-UAV.333https://github.com/Jacobkooi/InclinedDroneLander.git One
three-dimensional environment can be used with any compatible DRL algorithm to
train set-point tracking. The other two-dimensional environment, restricted to
the vertical $xz$-plane, can be used with an on-policy algorithm to train
inclined landing. The resulting policies adequately transfer to the real
Crazyflie.

  * •

Building upon the state-of-the-art model-free proximal policy optimization
(PPO) algorithm [12], we propose a powerful curriculum learning [13] approach
to speed up the convergence when using sparse rewards, without the need for
applying a multi-goal setting like in hindsight experience replay [14] or
iterated supervised learning [15].

  * •

We test the trained policy network in simulations and then deploy it onto the
real Crazyflie quadrotor to demonstrate the actual inclined landing in an
indoor flying arena.444https://youtu.be/53YaqfwUIFU

The remainder of the paper is structured as follows. We first give an overview
of the related work in Section II. The dynamic quadrotor model used for
simulation and training is described in Section III. Next, Section IV presents
the DRL simulation framework that is used to train inclined landing and set-
point tracking. Section V describes the simulation and lab setup and presents
the experimental results. Finally, in Section VI, the conclusions and
limitations of this work are given, along with proposals for future work.

##  II Related Work

Deep reinforcement learning methods have been applied to a variety of
quadrotor control problems, including hovering [3], attitude control [2], set-
point tracking and disturbance recovery [1, 4]. Specifically for landing, a
deep neural network was employed to learn higher-order interactions to
stabilize the near-ground behavior of a nonlinear quadrotor controller [5]. A
deep Q-learning network (DQN) was used to detect a marker symbol and perform a
landing by using a downward-facing low-resolution camera [7, 8]. Least-squares
policy iteration (LSPI) was employed to autonomously land on a marker [10] and
the deep deterministic policy gradient (DDPG) algorithm [16] was used to
navigate a descending quadrotor to land on a moving platform [9]. Finally, the
work in [6] involved a convolutional neural network to estimate the heading
angle to aid UAV landing in the case of a sensor failure. However, none of the
approaches considered inclined landing and none of the methods developed can
be directly applied to this problem.

Inclined landing has been the topic of several works outside the deep learning
control literature. A nonlinear hybrid controller was proposed in [17]. A
trajectory-tracking controller first guides the quadrotor above the landing
platform and then switches to an attitude-tracking controller to ensure that
the attitude of the quadrotor adjusts to the slope of the landing platform
upon touchdown. This is an ad hoc local strategy, incapable of generating
optimal landing trajectories from arbitrary initial conditions. Besides, no
real-time control experiments have been reported in this paper. The method
proposed in [18] features a nonlinear MPC to land a quadrotor on a moving
inclined surface. Real-time experimental results were reported, showing a
successful landing. The limitations of MPC are its computational complexity
and the difficulty of parameter tuning, especially of the prediction horizon,
which needs to be long for some of the landing trajectories, making the method
unsuitable for embedded implementation on the quadrotor. The approach
developed in [19] relies on splitting up the problem in the generation of
dynamically feasible trajectories, and their subsequent trajectory tracking.
Perching to slopes of up to 90 degrees has been demonstrated in lab
experiments. To keep the problem tractable, the authors break the desired
trajectory down into segments with a maximum duration of one second. The
overall approach is more complex than the nonlinear feedback policy approach
pursued in this paper.

##  III Simulation Model

The dynamic model of the Crazyflie 2.1 Nano-UAV is formed by the equations of
motion (EOM). We divide them into the Newton-Euler equations, which govern the
axial accelerations, and an approximation of the body attitude control loops.
The command input vector $u$ to the Crazyflie’s onboard controller is defined
as

$\displaystyle
u=\begin{bmatrix}\Theta_{c}&\phi_{c}&\theta_{c}&\dot{\psi}_{c}\end{bmatrix}^{T}.$
(1)

Here, $\Theta_{c}$ is the commanded pulse-width modulation (PWM) signal
representing the total thrust, $\phi_{c}$ and $\theta_{c}$ are the commanded
pitch and roll angles, respectively, and $\dot{\psi}_{c}$ is the commanded yaw
rate [20]. These inputs are bounded by

$\displaystyle u_{\min}$
$\displaystyle=\begin{bmatrix}10000&-30^{o}&-30^{o}&-200^{o}/s\end{bmatrix}^{T},$
(2) $\displaystyle u_{\max}$
$\displaystyle=\begin{bmatrix}60000&30^{o}&30^{o}&200^{o}/s\end{bmatrix}^{T}.$

###  III-A Newton-Euler Equations

The quadrotor is modeled as a rigid body, with the axial accelerations in the
inertial frame $\begin{bmatrix}x&y&z\end{bmatrix}^{T}$:

$\begin{bmatrix}m\ddot{x}\\\ m\ddot{y}\\\
m\ddot{z}\end{bmatrix}=\textbf{R}\left(\begin{bmatrix}0\\\ 0\\\
F_{t}\end{bmatrix}+F_{a}\right)+\begin{bmatrix}0\\\ 0\\\ -mg\end{bmatrix}$ (3)

Figure 1: The quadrotor coordinate system used throughout this paper.
Subscripts B and I represent the body and inertial frame, respectively, and
$F_{i}$ is the thrust due to rotor $i$. Adopted from [21].

with $m$ the quadrotor’s mass, R the rotation matrix from the body frame to
the inertial frame, $F_{t}$ the total thrust force and $F_{a}$ the drag force.
The rotation matrix corresponding to the axis frame representation in Fig. 1
is

$\textbf{R}=\begin{bmatrix}c\psi c\theta-s\phi s\psi s\theta&-c\phi
s\psi&c\psi s\theta+c\theta s\phi s\psi\\\ c\theta s\psi+c\psi s\phi
s\theta&c\phi c\theta&s\psi s\theta+c\psi c\theta s\phi\\\ -c\phi
s\theta&s\phi&c\phi c\theta\end{bmatrix}$ (4)

where $c$ is a cosine, $s$ is a sine, and $\phi$, $\theta$ and $\psi$ are the
roll, pitch and yaw angles, respectively. The relation between the commanded
PWM $\Theta_{c}$ in (1) and $F_{t}$ in (3) is modelled by using the discrete-
time transfer function found in [22] for an individual motor $i$:

$\displaystyle\frac{F_{i}(z)}{\Theta_{c,i}(z)}=\frac{7.2345374\cdot
10^{-8}}{1-0.9695404z^{-1}}\qquad(500Hz).$ (5)

Since the Crazyflie’s onboard controller only takes a single PWM signal for
all four motors, we assume that $F_{t}\approx 4F_{i}$ with
$\Theta_{c}\approx\Theta_{c,i}$. Multiplying (5) by four and digitally
converting it to continuous-time gives

$\displaystyle\begin{bmatrix}\dot{\Omega}\\\ F_{t}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}-15.467\\\ 1.425\cdot
10^{-4}\end{bmatrix}\Omega+\begin{bmatrix}1\\\ 2.894\cdot
10^{-7}\end{bmatrix}\Theta_{c}$ (6)

where $\Omega$ is an unmeasured state used for simulation purposes only. The
drag force $F_{a}$ in (3) is expressed as [23]

$\displaystyle F_{a}$ $\displaystyle=\boldsymbol{K}_{a}{\omega}_{\Sigma}v$ (7)

with $\omega_{\Sigma}$ the sum of the rotor velocities, $v$ the body-frame
velocity vector and $\boldsymbol{K}_{a}$ a diagonal matrix of drag constants
estimated in [22]. Because $\omega_{\Sigma}$ is not known during simulation,
we approximate it from $F_{t}$ with additional conversion formula’s given in
[22].

###  III-B Body Attitude Control Loops

The body attitude rates are modelled by equations that approximate the
dynamics of the attitude control loops [24]:

$\displaystyle\dot{\phi}$
$\displaystyle=\frac{1}{\tau_{\phi}}(k_{\phi}\phi_{c}-\phi),$ (8)
$\displaystyle\dot{\theta}$
$\displaystyle=\frac{1}{\tau_{\theta}}(k_{\theta}\theta_{c}-\theta),$
$\displaystyle\dot{\psi}$ $\displaystyle=\dot{\psi}_{c}.$

Here $\tau_{\phi}$, $\tau_{\theta}$ and $k_{\phi}$, $k_{\theta}$ are the time
and gain constants for roll and pitch, respectively. The yaw rate is assumed
to instantaneously track the desired yaw rate, which is a reasonable
assumption since yaw has no effect on the quadrotor’s position [24]. Because
the closed-loop dynamics are unknown, the parameters $k_{\theta}$, $k_{\phi}$,
$\tau_{\theta}$ and $\tau_{\phi}$ need to be identified. Given the quadrotor’s
symmetry, $k_{\theta}$ and $k_{\phi}$ as well as $\tau_{\theta}$ and
$\tau_{\phi}$ are assumed equal. These parameters are estimated by fitting
data gathered by a motion-capture system to the equations (8). We conducted 20
experiments using square and sine waves ranging from zero to thirty degrees,
which gave an average fit of $85.3\%$ using Matlab’s `nlgreyest` function,
with the resulting parameters $k_{\phi}=k_{\theta}=1.1094$ and
$\tau_{\phi}=\tau_{\theta}=0.1838$ s.

To simulate the quadrotor, the model equations (3) and (8) are integrated by
using the fourth-order Runge Kutta (RK4) method. The step size is fixed and
equal to the sampling period $T_{s}=0.02$ s.

##  IV Training Deep Reinforcement Learning Policies

To train the inclined landing, the quadrotor model of Section III is used as a
simulation environment for model-free DRL. Additionally, to navigate the
quadrotor, we train set-point tracking in the same fashion. Both policy
networks map quadrotor states to the desired control input in (1), which makes
them directly applicable to the real Crazyflie.

###  IV-A Preliminaries

The learning controller (agent) interacts with the model (environment) through
trials. The environment’s state-space is denoted by $\mathcal{S}$ and a
specific value of the state at time step $k$ by $s_{k}$. The agent applies an
action $a_{k}\in\mathcal{A}$ and subsequently receives a reward
$r_{k}\in\mathbb{R}$ , after which it observes the next state $s_{k+1}$. The
action $a_{k}$ is chosen by following a stochastic policy $\pi_{k}(a|s)$ or a
deterministic policy $\mu_{k}(s)$. This policy can be optimized in many
different ways. Most techniques maximize the discounted return
$\eta(\pi_{\phi})=\operatorname{\mathbb{E}}_{\tau}[\sum_{t=0}^{T}\gamma^{t}r(s_{k},a_{k})]$,
with $\tau$ a trajectory following the policy $\pi_{\phi}$.

###  IV-B Set-Point Tracking

We first train a three-dimensional set-point tracking policy network to
empirically check the simulation-to-reality performance of the DRL algorithms
and to fly to a desired starting position for inclined landing. For this task,
the states and actions are defined as follows:

$\displaystyle s_{3d}$
$\displaystyle=\begin{bmatrix}x&y&z&v_{x}&v_{y}&v_{z}&\phi&\theta\end{bmatrix}^{T},$
(9) $\displaystyle a_{3d}$
$\displaystyle=\begin{bmatrix}\Theta_{c}&\phi_{c}&\theta_{c}\end{bmatrix}^{T}.$

Here, $v$ is the velocity in the inertial frame. Note that the yaw angle is
kept constant at zero degrees, and can thus be omitted throughout all our
experiments. For set-point tracking, we use the following reward function:

$r_{k}=-e_{p}-0.2e_{v}-0.1e_{\phi,\theta}-0.1\frac{a_{\phi,\theta}^{2}}{\max(e_{p},0.001)}$
(10)

where $e_{p}$, $e_{v}$ and $e_{\phi,\theta}$ are Euclidean distance errors of
the position, velocity and orientation, respectively, with respect to the goal
state. The term $a_{\phi,\theta}^{2}$ is the sum of the squared roll and pitch
actions (normalized between $0$ and $1$). It is scaled by the reciprocal of
$e_{p}$ to minimize oscillations near the goal position.

The policy network is a fully connected neural network with two hidden layers,
with 64 neurons each, and the tanh activation function everywhere except for
the output layer which has a linear activation function. The final output is
subsequently clipped between $-1$ and $1$. We use the PPO algorithm [12] to
train the set-point tracking network. Other state-of-the-art DRL algorithms
like twin delayed deep deterministic policy gradient (TD3) [25] and soft
actor-critic (SAC) [26] converged successfully as well, but PPO was superior
in terms of the training time and the final policy performance.

###  IV-C Inclined Landing

The inclined landing is trained in the $xz$-plane, using the following states
and actions:

$\displaystyle s_{2d}$
$\displaystyle=\begin{bmatrix}x&z&v_{x}&v_{z}&\theta\end{bmatrix}^{T},$ (11)
$\displaystyle a_{2d}$
$\displaystyle=\begin{bmatrix}\Theta_{c}&\theta_{c}\end{bmatrix}^{T}.$

For the sake of brevity, in the sequel, we refer to $s_{2d}$ and $a_{2d}$ by
$s$ and $a$, respectively. Because the quadrotor is under-actuated, an initial
swinging motion away from the landing location is required for some initial
conditions. This characteristic is incompatible with the bias a Euclidean
distance based reward like the one in (10) generates. The reward function used
for inclined landing is therefore a sparse reward defined as follows:

$r_{k}=\begin{cases}0&\text{if}\quad s_{k}\in S_{g}\\\ -\beta&\text{if}\quad
s_{k}\in S_{o}\\\ -2&\text{if}\quad s_{k}\in S_{b}\\\
-1&\text{otherwise.}\end{cases}$ (12)

Here $S_{g}=\\{s\;|\;|s_{i}-s_{g,i}|<\delta_{g,i},\,\forall i\\}$ is the set
of goal states, defined as a hyperbox around the landing attitude. The goal
threshold vector $\delta_{g}$ defines the desired landing tolerance and it is
set by the user. The landing platform itself is an obstacle associated with a
set of obstacle states $S_{o}$ and a penalty $\beta$, and $S_{b}$ represents
the set of states close to the state space boundaries.

The use of a sparse reward requires extensive exploration in order to receive
a non-negative reward and leads to prolonged training. We introduce the
following curriculum learning [13] procedure to speed up the training:

  * •

The training starts without a landing platform and with a horizontal goal
state. Only once the quadrotor reliably reaches the horizontal goal, we begin
slightly tilting the goal position after each episode. Finally, the landing
platform is introduced into the environment, see Fig. 2.

  * •

We initialize simulations near the goal state and with each episode expand the
set of initial positions. This eliminates the need for exploration by letting
the non-negative rewards propagate throughout the value network at the
beginning of training.

  * •

We start with a large goal hyperbox $S_{g}$ and as the training progresses,
the hyperbox is gradually reduced to its desired size.

This learning curriculum requires an on-policy learning algorithm, such as
PPO. Off-policy replay buffers would inevitably contain samples representing
goals that are no longer relevant. In our experience, off-policy algorithms
TD3 and SAC cannot keep up with the curriculum. The policy network
architecture is similar to the one used for set-point tracking. The input and
output layers for inclined landing are the state and action vectors in (11).

Figure 2: The progress of the training curriculum for inclined landing. The
red quadrotor represents the agent and the black quadrotor represents the goal
state.

###  IV-D Simulation to Reality Transfer

To deploy the policies on the Crazyflie 2.1 Nano-UAV, the trained Pytorch
network is converted to a Robot Operating System (ROS) node, which maps the
state to the control input vector in (1). The positions and velocities come
from a Kalman filter node, which appends the coordinates from the Optitrack
motion capture system with the estimated inertial frame velocities. The
quadrotor’s orientation is taken from the Crazyflie’s default onboard
estimator. An overview of this process is given in Fig. 3.

Figure 3: Schematic overview of the Pytorch model deployment.

##  V Experiments

All simulations, DRL training and lab experiments are done on an HP Zbook
Studio G4 laptop, with the default Nvidia Quadro M1200 GPU and an Intel Core
i7-7700HQ CPU. The additional hardware used is the Optitrack motion-capture
system, the standard Crazyflie 2.1 with a small marker holder and a Crazyradio
PA for communication with the Crazyflie.

###  V-A Experimental Setup

The simulations and DRL training are implemented in Python, using the Pytorch
package. The quadrotor dynamics are simulated using the EOM of Section III,
integrated by the RK4 method. To increase the integration speed, we compile
the EOM functions using Numba’s Just-in-time (jit) package. By keeping the
simulations and rendering in our own Gym-architecture [11] environment, we can
use well-tuned implementations of existing model-free algorithms by [27]. The
simulation bounds are equal to the dimensions of the actual flying arena with
$\begin{bmatrix}x_{min}&x_{max}&y_{min}&y_{max}&z_{min}&z_{max}\end{bmatrix}=\begin{bmatrix}-3.4&3.4&-1.4&1.4&0&2.4\end{bmatrix}$
m. The actions represent the multiplication of the clipped policy network
outputs $(a\in[-1,1])$ by the control bounds in (2), with an exception for the
PWM network output which is multiplied by $16500$ and added to the estimated
hover PWM of $42000$. All simulations use the control sampling frequency of
$50$ Hz, with episode lengths of 300 time steps, i.e., 6 seconds.

Figure 4: Landing trajectories starting from above right (top) and from
directly above (bottom) of the inclined goal position at (0, 1.25).

####  V-A1 Set-Point Simulations

The goal state is taken as
$s_{g}=\begin{bmatrix}0&0&1.2&0&0&0&0&0\end{bmatrix}^{T}$, which represents
the center of our physical flying arena. After each episode, the initial
position is set randomly anywhere within the simulation bounds, with a small
margin around the edges. The PPO implementation of [27] is used, with the
discount factor $\gamma$ reduced to $0.97$ to account for more short-term
control behaviour. We train for a total of $10^{6}$ time steps representing 27
minutes of training time. For navigation with the resulting policy, a
coordinate change suffices to fly the quadrotor to an arbitrary position.

####  V-A2 Inclined Landing Simulations

For inclined landing, we operate the quadrotor in the $xz$-plane. The DRL
algorithm used is the PPO implementation of [27], where minor changes have
been made to activate rendering every 50 training iterations and to gradually
start increasing $\gamma$ from $0.97$ to $0.99$ after 300 training iterations.
An episode ends at 300 time steps or when the state $x$ is within the goal
hyperbox $S_{g}$. The landing platform is modeled as a polygon and appears
after $8\cdot 10^{5}$ time steps. The goal threshold vector is defined as
$\delta_{g}=\begin{bmatrix}\delta_{x}&\delta_{z}&\delta_{v_{x}}&\delta_{v_{z}}&\delta_{\theta}\end{bmatrix}=\begin{bmatrix}d&d&min(10d,1.5)&min(10d,1.5)&0.25d\end{bmatrix}$,
where $d$ starts at $0.25$ in the beginning of training and gradually
decreases to $0.10$ after $2500$ episodes with $0.15/5000$ per episode. The
box of possible starting positions around the goal state expands with every
episode by $1/6000$ m in the $x$-direction and by $1/8000$ m in the
$y$-direction. Additionally, the goal state stays horizontal for the first
$4\cdot 10^{5}$ time steps and then gradually tilts towards its final
inclination of $-\pi/7$ at the rate of $(-\pi/7)/6000$ radians per episode.
The final goal state is set at
$s_{g}=\begin{bmatrix}0&1.25&0&0&-\pi/7\end{bmatrix}^{T}$. The obstacle reward
constant $\beta$ is taken as $-7$, which was found empirically to be the right
trade-off between the goal of landing on top of the platform and the necessity
to avoid the landing platform base. Training is stopped anywhere between
$1.2\cdot 10^{6}$ and $3\cdot 10^{6}$ time steps (30 to 80 minutes), when the
rendering shows that the policy executes the inclined landing reliably. Note
that rather than monitoring the loss function, the aforementioned parameters
and curricula have been empirically tuned by frequently rendering. The trend
in the loss function value is quite meaningless, given the curricula and the
discount factor adaptation.

####  V-A3 Validation on a real quadrotor

The quadrotor used is the Crazyflie 2.1 Nano-UAV with its original firmware. A
Crazyflie-specific package [20] allows us to publish the control values in (1)
directly to the ROS server, which are then transmitted with low latency to the
Crazyflie over the Crazyradio PA. The trained policy network is evaluated at
$80$ Hz, even though it was trained at $50$ Hz. This is possible, as the
policy is a function of the physical states only, and can therefore be
evaluated at any frequency. A single policy evaluation takes approximately
$2.5$ ms. The coordinates from the Optitrack motion-capture system come in at
$120$ Hz and are combined with the inertial frame velocities by the Kalman
filter node [28]. The orientation is taken from the Crazyflie’s onboard
estimator, which is received through the Crazyradio at around $80$Hz. The
obtained position, velocity and orientation estimates form the policy input.
The experiment itself starts with the drone flying to a position of choice,
using the three-dimensional set-point tracking network. Once it is positioned,
the networks are switched and the drone commences the inclined landing. This
can begin from an arbitrary location in the top half of the flying arena. The
landing trial ends when the quadrotor’s state $s_{k}$ is located in the set of
goal states $S_{g}$.

###  V-B Experiment Results

Early experiments showed a vertical offset between the simulated and real
trajectories, caused by a slight inaccuracy of our motor thrust equation (6).
We were able to compensate for this by increasing the hover PWM to 48000
during real testing. The resulting behaviour was very similar to the
simulations, as can be seen in Fig. 4.555Because of a faulty Optitrack
measurement, the bottom figure shows a single misplaced red quad with a
corresponding short drop in PWM output. These trajectories originate from the
same policy network, starting from arbitrary initial positions and ending when
$s_{k}\in S_{g}$ with
$\delta_{g}=\begin{bmatrix}0.10&0.10&1.5&1.5&0.025\end{bmatrix}$.

To further evaluate the performance, we measured the landing success rate when
starting the flight from three different initial positions. Each position was
evaluated $10$ times, resulting in the success percentages reported in Table
I. Even if an experiment did not succeed, the agent would not crash but
autonomously fly back and forth, sometimes succeeding in its second or third
attempt. However, these further attempts were not counted as a successful
landing.

TABLE I: Comparison of Real-World and Simulation Experiments setting |  | success from initial $(x,z)$ |  | total
---|---|---|---|---
| $(0,2)$ | $(-1.5,1.6)$ | $(1.5,1.8)$ |
Real-World | 90% | 70% | 100% | 86.7%
Simulation | 90% | 90% | 100% | 93.3%

No simulation-to-reality transfer techniques, such as domain randomization,
have been employed, as they did not seem to improve final performance, while
they did complicate the agent’s training. Additionally, we found that using
the Crazyflie’s onboard orientation estimates rather than the Optitrack
orientation estimates resulted in a substantial increase in performance.

The results show that inclined landing controllers can be designed by means of
DRL. These controllers can transfer adequately to reality without the need for
dynamics randomization or sensor noise. Furthermore, the same policy can be
initiated from a wide variety of starting states. Results could further
improve by additional system identification of the total motor thrust, where a
slight mismatch might have been caused by making too strong assumptions about
the scaling of the single motor model in (5) to the full motor model in (6). A
perfect thrust model will however not exist, due to motor wear and tear and
the strong relation of the drone’s battery level to its thrust output, which
is also reported in [1].

##  VI CONCLUSION AND FUTURE WORK

We have presented a model-free DRL technique to facilitate autonomous
quadrotor landing on an inclined surface. We trained a control agent with PPO,
using sparse rewards and a learning curriculum. This allows the agent to
gradually progress toward a more difficult task, i.e., larger inclination
angles of the landing platform. Moreover, we have shown that the trained
policies transfer well to reality, without employing any simulation-to-reality
transfer techniques.

A limitation of this work is the fact that we restricted the landing
trajectory to the xz-plane, which may cause some drift in the y-direction. A
three-dimensional landing policy would increase precision, albeit at the cost
of longer and more complex training. Furthermore, the reliance on external
motion-capture systems makes this work hard to reproduce with onboard sensors
only. An extension to such a setting is the topic of our future work.
Additional future work could focus on larger inclination angles of the landing
platform, as long as the onboard attitude controller would allow them. In this
way, one could try to extend the findings in this paper to perching behavior
similar to [19]. The goal inclination can also be added to the quadrotor’s
state, which would enable landing on unknown platforms, using, for instance,
the laser measurement system of [17] or an onboard camera system like in [18].
Also, a form of the platform contact dynamics in [29] could be implemented in
our simulator for a more robust landing and to aid the design of an end-to-end
controller, which would eliminate the need for an external stopping signal.

## ACKNOWLEDGMENT

Robert Babuška was supported by the European Union’s H2020 project Open Deep
Learning Toolkit for Robotics (OpenDR) under grant agreement No. 871449.

## References

  * [1] J. Hwangbo, I. Sa, R. Siegwart, and M. Hutter, “Control of a quadrotor with reinforcement learning,” _IEEE Robotics and Automation Letters_ , vol. PP, pp. 1–1, 06 2017.
  * [2] W. Koch, R. Mancuso, R. West, and A. Bestavros, “Reinforcement learning for uav attitude control,” _ACM Trans. Cyber-Phys. Syst._ , vol. 3, no. 2, Feb. 2019.
  * [3] N. O. Lambert, D. S. Drew, J. Yaconelli, R. Calandra, S. Levine, and K. S. J. Pister, “Low level control of a quadrotor with deep model-based reinforcement learning,” _IEEE Robotics and Automation Letters_ , vol. 4, pp. 4224–4230, 2019.
  * [4] A. Molchanov, T. Chen, W. Hönig, J. A. Preiss, N. Ayanian, and G. S. Sukhatme, “Sim-to-(multi)-real: Transfer of low-level robust control policies to multiple quadrotors,” _arXiv preprint arXiv:1903.04628_ , 2019\.
  * [5] G. Shi, X. Shi, M. O’Connell, R. Yu, K. Azizzadenesheli, A. Anandkumar, Y. Yue, and S. Chung, “Neural lander: Stable drone landing control using learned dynamics,” in _2019 International Conference on Robotics and Automation (ICRA)_ , 2019, pp. 9784–9790.
  * [6] Y. Bicer, M. Moghadam, C. Sahin, B. Eroglu, and N. K. Üre, “Vision-based uav guidance for autonomous landing with deep neural networks,” in _AIAA Scitech 2019 Forum_ , 2019, p. 0140.
  * [7] R. Polvara, M. Patacchiola, S. Sharma, J. Wan, A. Manning, R. Sutton, and A. Cangelosi, “Toward end-to-end control for uav autonomous landing via deep reinforcement learning,” in _2018 International Conference on Unmanned Aircraft Systems (ICUAS)_ , 2018, pp. 115–123.
  * [8] R. Polvara, M. Patacchiola, M. Hanheide, and G. Neumann, “Sim-to-real quadrotor landing via sequential deep q-networks and domain randomization,” _Robotics_ , vol. 9, no. 1, 2020.
  * [9] A. Rodríguez Ramos, C. Sampedro Pérez, H. Bavle, P. de la Puente, and P. Campoy, “A deep reinforcement learning strategy for uav autonomous landing on a moving platform,” _Journal of Intelligent & Robotic Systems_, vol. 93, 02 2019.
  * [10] M. B. Vankadari, K. Das, C. Shinde, and S. Kumar, “A reinforcement learning approach for autonomous control and landing of a quadrotor,” in _2018 International Conference on Unmanned Aircraft Systems (ICUAS)_ , 2018, pp. 676–683.
  * [11] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba, “Openai gym,” _arXiv preprint arXiv:1606.01540_ , 2016.
  * [12] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017\.
  * [13] Y. Bengio, J. Louradour, R. Collobert, and J. Weston, “Curriculum learning,” in _Proceedings of the 26th Annual International Conference on Machine Learning_ , ser. ICML ’09. New York, NY, USA: Association for Computing Machinery, 2009, p. 41–48.
  * [14] M. Andrychowicz, F. Wolski, A. Ray, J. Schneider, R. Fong, P. Welinder, B. McGrew, J. Tobin, O. Pieter Abbeel, and W. Zaremba, “Hindsight experience replay,” in _Advances in Neural Information Processing Systems_ , vol. 30, 2017.
  * [15] D. Ghosh, A. Gupta, A. Reddy, C. M. Devin, B. Eysenbach, and S. Levine, “Learning to reach goals via iterated supervised learning,” in _International Conference on Learning Representations_ , 2021.
  * [16] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning.” in _ICLR_ , 2016.
  * [17] J. Dougherty, D. Lee, and T. Lee, “Laser-based guidance of a quadrotor UAV for precise landing on an inclined surface,” _2014 American Control Conference_ , pp. 1210–1215, 2014.
  * [18] P. Vlantis, P. Marantos, C. P. Bechlioulis, and K. J. Kyriakopoulos, “Quadrotor landing on an inclined platform of a moving ground vehicle,” in _2015 IEEE International Conference on Robotics and Automation (ICRA)_ , 2015, pp. 2202–2207.
  * [19] J. Thomas, M. Pope, G. Loianno, E. Hawkes, M. Estrada, H. Jiang, M. Cutkosky, and V. Kumar, “Aggressive flight for perching on inclined surfaces,” _Journal of Mechanisms and Robotics_ , vol. 8, 2015.
  * [20] W. Hönig and N. Ayanian, _Flying Multiple UAVs Using ROS_. Springer International Publishing, 2017, pp. 83–118.
  * [21] Y. Chen and N. O. Pérez-Arancibia, “Nonlinear adaptive control of quadrotor multi-flipping maneuvers in the presence of time-varying torque latency,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2018, pp. 1–9.
  * [22] J. Förster and R. D’Andrea, “System identification of the crazyflie 2.0 nano quadrocopter,” Zürich, Tech. Rep., 2015.
  * [23] M. W. Mueller, M. Hamer, and R. D’Andrea, “Fusing ultra-wideband range measurements with accelerometers and rate gyroscopes for quadrocopter state estimation,” in _2015 IEEE International Conference on Robotics and Automation (ICRA)_ , 2015, pp. 1730–1736.
  * [24] M. Kamel, J. Alonso-Mora, R. Siegwart, and J. Nieto, “Robust collision avoidance for multiple micro aerial vehicles using nonlinear model predictive control,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2017, pp. 236–243.
  * [25] S. Fujimoto, H. van Hoof, and D. Meger, “Addressing function approximation error in actor-critic methods,” in _ICML_ , 2018, pp. 1582–1591.
  * [26] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,” in _Proceedings of the 35th International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, vol. 80, 10–15 Jul 2018, pp. 1861–1870.
  * [27] A. Raffin, A. Hill, M. Ernestus, A. Gleave, A. Kanervisto, and N. Dormann, “Stable baselines 3,” https://github.com/DLR-RM/stable-baselines3, 2019\.
  * [28] H. Zhu and J. Alonso-Mora, “Chance-Constrained Collision Avoidance for MAVs in Dynamic Environments,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 776–783, 2019.
  * [29] J. Bass and A. L. Desbiens, “Improving multirotor landing performance on inclined surfaces using reverse thrust,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 5850–5857, 2020.


# Multijet topology in high-energy nuclear collisions: jet broadening

Jin-Wen Kang  Key Laboratory of Quark & Lepton Physics (MOE) and Institute of
Particle Physics, Central China Normal University, Wuhan 430079, China  Lei
Wang  Key Laboratory of Quark & Lepton Physics (MOE) and Institute of Particle
Physics, Central China Normal University, Wuhan 430079, China  Wei Dai  School
of Mathematics and Physics, China University of Geosciences, Wuhan 430074,
China  Sa Wang  Guangdong Provincial Key Laboratory of Nuclear Science,
Institute of Quantum Matter, South China Normal University, Guangzhou 510006,
China  Ben-Wei Zhang<EMAIL_ADDRESS>Key Laboratory of Quark &
Lepton Physics (MOE) and Institute of Particle Physics, Central China Normal
University, Wuhan 430079, China

###### Abstract

This work presents the first theoretical investigation of the medium
modification of jet broadening as an event-shape observable in multijet final
states due to jet quenching in high-energy nuclear collisions. The partonic
spectrum of $pp$ collisions with next-to-leading order (NLO) accuracy at
$\sqrt{s_{\mathrm{NN}}}=5.02$ TeV is provided by the Powheg+Pythia8 event
generator, and the linear Boltzmann transport (LBT) model is utilized to
investigate the energy loss of fast partons as they traverse through the hot
and dense QCD medium. Jet broadening distributions in multijet final states
for both $pp$ and PbPb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV are
calculated. We observe an enhancement at the small jet broadening region while
a suppression at the large jet broadening region in PbPb collisions relative
to that in $pp$. This suggests that medium modifications with parton energy
loss in the QGP lead to a more concentrated energy flow in all observed
multijet events in PbPb reactions. We also demonstrate that the intertwining
of two effects, the jet number reduction and the restructured contribution,
results in the novel behavior of nuclear modification of the jet broadening
observable in PbPb collisions.

##  I Introduction

In heavy-ion collisions (HICs) at the RHIC and the LHC, a new form of nuclear
matter, the quark-gluon plasma (QGP) Busza  _et al._ (2018); Yagi  _et al._
(2005) composed of de-confined quarks and gluons, may be created, and the jet
quenching (or parton energy loss effect) has long been proposed as an
excellent hard probe to unravel the fascinating properties of the QCD medium
at extreme high density and temperature Bjorken (1982); Gyulassy  _et al._
(2004); Wang and Gyulassy (1992); Qin and Wang (2015); Abelev  _et al._
(2010); Chatrchyan  _et al._ (2011); Zhang  _et al._ (2004). It has been
extensively studied that parton energy loss effect may result in variations of
leading hadron productions such as the suppression of inclusive hadron
spectra, the disappearance of away-side di-hadron correlations momentum and
other hadron correlations Vitev and Gyulassy (2002); Vitev (2005); Zhang  _et
al._ (2007, 2009); Xu and Ko (2011); Eyyubova  _et al._ (2015); Dai  _et al._
(2015); Liu  _et al._ (2016); Dai  _et al._ (2017); Cao  _et al._ (2018);
Albacete  _et al._ (2019); Ma  _et al._ (2019); Pandey and Singh (2023). As
the running of heavy-ion programs at the LHC with its unprecedented colliding
energies, the study of parton energy loss has been extended to medium
modifications of many jet observables Bass  _et al._ (2009); Vitev  _et al._
(2008); Vitev and Zhang (2010); Majumder and Van Leeuwen (2011); He  _et al._
(2012); Apolinario  _et al._ (2013); Burke  _et al._ (2014); Senzel  _et al._
(2015); Blaizot and Mehtar-Tani (2015); Casalderrey-Solana  _et al._ (2016);
Milhano and Zapp (2016); Sirunyan  _et al._ (2018); Connors  _et al._ (2018);
Tachibana  _et al._ (2017); Dai  _et al._ (2020); Zhang  _et al._ (2019);
Aaboud  _et al._ (2019); Wang  _et al._ (2019); Cao and Wang (2021); Acharya
_et al._ (2022); Cunqueiro and Sickles (2022); Wang  _et al._ (2021); Chen
_et al._ (2022); Kumar  _et al._ (2023); Dai  _et al._ (2022); Wang  _et al._
(2023a, b). However, so far most studies of jets in HICs have focused on
inclusive jets, dijet/$Z+\mathrm{jet}$/$\gamma+\mathrm{jet}$ observables, or
their substructures, while only a very limited attention has been paid to the
global geometrical properties of multijet events Chen  _et al._ (2020);
Mallick  _et al._ (2020). It is of interest to further explore how parton
energy loss influences the global geometry of multijet events and investigate
their nuclear modifications with the existence of the QGP.

The geometrical properties of the energy flow of the collisions are usually
described by event-shape variables. The event-shape variable is a general term
for a large class of observables, including the thrust, the jet broadening,
the jet mass, the third-jet resolution parameter, the sphericity, the
$\mathcal{F}$-parameter, the 1-jettiness, the Fox-Wolfram moments, and the
aplanarity, etc Sjostrand  _et al._ (2006); Banfi  _et al._ (2010); Aad  _et
al._ (2012); Khachatryan  _et al._ (2014); Aad  _et al._ (2016, 2021); Hessler
_et al._ (2022). These event-shape variables are not only used to extract the
strong coupling $\alpha_{s}$ from properties of the final-state Abe  _et al._
(1991); Catani  _et al._ (1992); Abreu  _et al._ (1997); Biebel  _et al._
(1999); Abdallah  _et al._ (2003); Wang  _et al._ (2022), but also tune/test
the parton showers and non-perturbative components of Monte Carlo event
generators Abreu  _et al._ (1996); Sjostrand  _et al._ (2006); Kundu  _et al._
(2019); Ghosh  _et al._ (2022). Previous measurements of event-shape variables
have been performed in electron-positron, deep inelastic scattering (DIS), or
proton-antiproton collisions. Recently, some experimenters have focused on the
event-shape variables at large transverse momentum transfer in $pp$ collisions
Aad  _et al._ (2012); Khachatryan  _et al._ (2014); Aad  _et al._ (2016,
2021). Meanwhile, some theorists have begun to explore the event-shape
variables in heavy-ion collisions at the LHC energies. In Refs. Chen  _et al._
(2020); Mallick  _et al._ (2021, 2020); Prasad  _et al._ (2022), the authors
studied the transverse sphericity in heavy-ion collisions at the LHC energies.
It is noted though for transverse sphericity two-jet events may give a large
contribution, for jet broadening as an event-shape observable defined to
utilize the transverse thrust axis the contribution of two-jet events turns
out to be zero. Therefore, jet broadening should provide relatively more
complete information on the geometric properties of the energy flow of
multijet events with jet number ${N_{\mathrm{jet}}}\ignorespaces\geqslant 3$.

In this paper, we present the first theoretical study on the normalized
distribution of jet broadening (denoted as ${B_{tot}}$) in PbPb collisions. We
employ Powheg+Pythia8 Nason (2004); Frixione  _et al._ (2007); Alioli  _et
al._ (2010, 2011); Sjöstrand  _et al._ (2015), a Monte Carlo model matching
NLO matrix elements with parton shower, to obtain baseline results of jet
broadening distributions in $pp$ collisions. Our model calculations of jet
broadening distribution could provide a decent description of CMS data in $pp$
collisions. Then we calculate the medium modification of the jet broadening
distribution in multijet (${N_{\mathrm{jet}}}\ignorespaces\geqslant 3$) final-
states in PbPb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02~{}\mathrm{TeV}$ for
the first time. We show the normalized distribution of jet broadening
$B_{tot}$ is enhanced at small $B_{tot}$ region in PbPb collisions (while
suppressed at larger ${B_{tot}}$) as compared to that in $pp$ collisions. It
indicates the selected events with lower jet broadening $B_{tot}$ increase due
to the jet quenching effect. In other words, the energy flow patterns of the
collision events in PbPb become more concentrated after including parton
energy loss effect.

This paper is organized as follows. Sec. II describes details of the
theoretical frameworks, including the definition of the observables, the
generation of the $pp$ baseline, and the energy loss model. In Sec. III, we
present the numerical results of the normalized jet broadening distributions
in PbPb collisions, and discuss how jet quenching impacts the jet broadening
in heavy-ion collisions. We give a brief summary in Sec. IV, and also include
an Appendix, which provides a detailed calculation of the concerned
observables in the ideal symmetric multijet limit.

##  II Framework

In this work, we consider multijet event-shape observables that measure the
properties of the energy flow in the final-states of high-energy particle
collisions. Because performing measurements near the beam are complex, one way
to address this difficulty is to define event-shapes using only jets in a
central region. Here, we utilize those jets in the central rapidity region to
calculate the event-shapes, and we select those jets that satisfy
$|\eta_{\mathrm{jet}}|<2.4$.

The transverse thrust $\tau_{\perp}$ is defined as Aad  _et al._ (2012);
Khachatryan  _et al._ (2014); Aad  _et al._ (2021)

${\tau_{\perp}}\equiv
1-\max_{\hat{n}}\frac{\sum_{i}\left|\vec{p}_{T,i}\cdot\hat{n}_{T}\right|}{\sum_{i}p_{T,i}},$
(1)

where $\hat{n}_{T}$ is the transverse thrust axis and it is the unit vector
that maximizes the sum of the projections of $\vec{p}_{T,i}$, $p_{T,i}$
represents the transverse momentum of the $i$-th jet. This variable is
sensitive to the modeling of two-jet and multijet ($\geqslant
3~{}\mathrm{jet}$) topologies. In a perfectly balanced back-to-back two-jet
event, $\tau_{\perp}$ is zero, and in an isotropic multijet event, it closes
to $1-2\pi$ Khachatryan  _et al._ (2014); Banfi  _et al._ (2010) (see more
details in the Appendix A). With a transverse thrust axis $\hat{n}_{T}$, one
can separate the region $\mathcal{C}$ into an upper (U) side $\mathcal{C}_{U}$
consisting of all jets in $\mathcal{C}$ with $\vec{p}_{T}\cdot\hat{n}_{T}>0$
and a lower (L) side $\mathcal{C}_{L}$ with $\vec{p}_{T}\cdot\hat{n}_{T}<0$.
Given the event region $X$, we can introduce the mean transverse-energy
weighted pseudo-rapidity $\eta_{X}$ and azimuthal angle $\phi_{X}$ of this
region Khachatryan  _et al._ (2014); Aad  _et al._ (2021); Banfi  _et al._
(2010),

$\eta_{X}\equiv\frac{\sum_{i\in\mathcal{C}_{X}}p_{T,i}\eta_{i}}{\sum_{i\in\mathcal{C}_{X}}p_{T,i}},\quad\phi_{X}\equiv\frac{\sum_{i\in\mathcal{C}_{X}}p_{T,i}\phi_{i}}{\sum_{i\in\mathcal{C}_{X}}p_{T,i}},$
(2)

where $X$ refers to upper (U) or lower (L) side. Then, the jet broadening
variable in two regions is defined as

$B_{X}\equiv\frac{1}{2P_{T}}\sum_{i\in\mathcal{C}_{X}}p_{T,i}\sqrt{\left(\eta_{i}-\eta_{X}\right)^{2}+\left(\phi_{i}-\phi_{X}\right)^{2}},$
(3)

where $P_{T}$ is the scalar sum of the transverse momenta of all the jets in
the region. From which one can obtain total jet broadening,

$B_{tot}\equiv B_{U}+B_{L}.$ (4)

Thus $B_{tot}=0$ for a two-jet event with definition.

Figure 1: Schematic illustrations of multijet event configuration in the
transverse plane and the values of their jet broadening observable. The red
arrow in the above figures represents the thrust axis, which is originally a
unit vector, but we have magnified it to make it easier to see in the plots.
It is easy to see that the sign of the $\hat{n}_{T}$ vector only affects the
identification of the upper (or lower) region, and does not change the value
of total jet broadening.

We may observe that the descriptions of $\eta_{X}$ and $\phi_{X}$ given in Eq.
(2) are similar to the definition of the center of mass, so we can in some
sense regard $(\eta_{X},\phi_{X})$ as the coordinates of the center of jets
energy outflow. Therefore, jet broadening can characterize the weighted
broadening of the jets relative to the center of the outgoing energy flow. We
plot three multijet event configurations in the transverse plane shown in Fig.
1 to demonstrate the physical picture of the jet broadening. In Fig. 1 we
suppress the variation of pseudo-rapidity for simplicity by imposing
$\eta_{i}=\eta_{X}=0$, in which case ${B_{tot}}$ has a maximum value of
$\pi/8$ for the circular limit (see the Appendix A). The only difference
between the events described in Fig. 1(b) and Fig. 1(a) is that the angle
between the second and third jets is larger, which makes the total jet
broadening of the event configuration shown in Fig. 1(b) larger than that
shown in Fig. 1(a). The event in Fig. 1(c) has one more jet than in Fig. 1(b),
and thus the total jet broadening of the event configuration shown in Fig.
1(c) is larger than the one in Fig. 1(b). From Fig. 1, we can find that the
weighted broadening of energy flow can also be used to distinguish the balance
of the spatial distribution of energy flow. Briefly speaking, the multijet
energy flow broadening is very small when the jet broadening variable tends to
zero; at this point, the spatial distribution of energy flow tends to be very
imbalanced. On the contrary, the multijet energy flow broadening increases
when the jet broadening variable is away from zero; the spatial distribution
of energy flow walks in a more balanced direction. By the way, we may use the
term “more concentrated” to express the meaning of “less broadening” or
“narrower” of jet broadening distributions.

A NLO$+$PS Monte Carlo event generator is employed in this research to
simulate jet productions in $pp$ collisions. Specifically, the NLO matrix
elements for the QCD dijet processes are provided by Powheg matches with the
final-state parton showering in Pythia8 Nason (2004); Frixione  _et al._
(2007); Alioli  _et al._ (2010, 2011); Sjöstrand  _et al._ (2015). For every
event, jets are reconstructed with the anti-$k_{t}$ algorithm and distance
parameter $R=0.5$ using FastJet Cacciari  _et al._ (2008, 2012). In
calculations, we make the kinematic cuts adopted in the CMS publication
Khachatryan  _et al._ (2014): the selected events are required to include at
least three jets in the central pseudo-rapidity region
$|\eta_{\mathrm{jet}}|<2.4$, the lower threshold of the reconstructed jets
$p_{T}$ is $30{\mathrm{~{}GeV}}$, and the leading jet $p_{T}$ satisfies
$110<p_{T,1}<170{\mathrm{~{}GeV}}$.

Figure 2: The jet broadening $B_{tot}$ are calculated at 7 TeV by
Powheg$+$Pythia8 in $pp$ collisions and compared to CMS experimental data.

Fig. 2 illustrates the multijet event-shape by Powheg$+$Pythia8 and the
comparison with the CMS data Khachatryan  _et al._ (2014) in $pp$ collisions
at 7 TeV. The simulation using Powheg$+$Pythia8 framework can provide a pretty
good description of the experimental data in $pp$ collisions. Consequently, it
may provide a reliable baseline for further study of medium modifications of
event-shape. In addition, we can see the effect of hadronization on the jet
broadening is not significant when $\ln(B_{tot})>-2.5$, but make give
considerable correction in small $\ln(B_{tot})$ region. In our subsequent
calculations, we reconstruct the jets at the hadron level to include
hadronization corrections.

In this study, we employ the Linear Boltzmann Transport (LBT) model to
consider both elastic and inelastic scattering processes of the initial jet
shower partons and the thermal recoil partons in the QGP medium Wang and Zhu
(2013); He  _et al._ (2015); Cao  _et al._ (2016). The LBT model has been
widely used in the study of jet quenching and successfully describes many
experimental measurements Chen  _et al._ (2020); Zhang  _et al._ (2019); Luo
_et al._ (2018).

The LBT model uses the linear Boltzmann transport equation to describe the
elastic scattering Wang and Zhu (2013); He  _et al._ (2015); Cao  _et al._
(2016),

$\displaystyle p_{1}\cdot\partial_{t}f_{1}(p_{1})$ $\displaystyle=$
$\displaystyle-\int\frac{d^{3}p_{2}}{(2\pi)^{3}2E_{2}}\int\frac{d^{3}p_{3}}{(2\pi)^{3}2E_{3}}\int\frac{d^{3}p_{4}}{(2\pi)^{3}2E_{4}}$
(5)
$\displaystyle\times\frac{1}{2}\sum_{2(3,4)}\left(f_{1}f_{2}-f_{3}f_{4}\right)\left|\mathcal{M}_{12\to
34}\right|^{2}\left(2\pi\right)^{4}$ $\displaystyle\times
S_{2}(\hat{s},\hat{t},\hat{u})\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4}),$

where $f_{i}(p_{i})(i=2,4)$ are the parton phase-space distributions in the
thermal medium with the local temperature and the fluid velocity;
$f_{i}(p_{i})(i=1,3)$ are the phase-space densities for the jet shower parton
before and after scattering, the leading-order (LO) elastic scattering matrix
elements $\left|\mathcal{M}_{12\to 34}\right|^{2}$ (Ref. Eichten  _et al._
(1984) for massless light partons and Ref. Combridge (1979) for heavy quarks)
may be collinear divergent for massless partons, and then are regulated by a
Lorentz-invariant regularization condition Auvinen  _et al._ (2010),

$S_{2}(\hat{s},\hat{t},\hat{u})\equiv\theta\left(\hat{s}\geqslant
2\mu_{D}^{2}\right)\theta\left(-\hat{s}+\mu_{D}^{2}\leqslant\hat{t}\leqslant-\mu_{D}^{2}\right),$
(6)

where $\hat{s}$, $\hat{t}$ and $\hat{u}$ are Mandelstam variables,
$\mu_{D}^{2}$ is the Debye screening mass.

The inelastic scattering is described by the higher-twist approach Guo and
Wang (2000); Zhang and Wang (2003); Zhang  _et al._ (2004); Majumder (2012),

$\frac{dN_{g}}{dxdk_{\perp}^{2}dt}=\frac{2\alpha_{s}C_{A}P(x)\hat{q}}{\pi
k_{\perp}^{4}}\left(\frac{k_{\perp}^{2}}{k_{\perp}^{2}+x^{2}m^{2}}\right)\sin^{2}\left(\frac{t-t_{i}}{2\tau_{f}}\right),$
(7)

where $x$ is the energy fraction of the radiated gluon relative to parent
parton with mass $m$, $k_{\perp}$ is the transverse momentum of the radiated
gluon, $P(x)$ is the splitting function in vacuum, $\hat{q}$ is the jet
transport coefficient He  _et al._ (2015), and $\tau_{f}$ is the formation
time of the radiated gluons in QGP medium, with
$\tau_{f}=2Ex(1-x)\big{/}\left(k_{\perp}^{2}+x^{2}m^{2}\right)$.

The 3$+$1D CLVisc hydrodynamical model Pang  _et al._ (2012, 2015) is employed
to generate a dynamically evolving bulk medium with initial conditions from
the AMPT model Lin  _et al._ (2005). Parameters used in the CLVisc are chosen
to reproduce bulk hadron spectra with experimental data.

In the LBT model, the strong coupling constant $\alpha_{s}$ is fixed and is
the only parameter that needs to be specified to control the strength of
parton interaction. Referring to the published works Zhang  _et al._ (2019);
Luo  _et al._ (2018); Chen  _et al._ (2020), we take $\alpha_{s}$ to be 0.2 in
this work.

##  III Results and Analysis

Figure 3: The normalized distributions of the jet broadening $B_{tot}$ in
0-30% PbPb and $pp$ collisions at
$\sqrt{s_{\mathrm{NN}}}=5.02~{}\mathrm{TeV}$. The ratio of the two normalized
distributions is given in the lower panel. The error bars (vertical lines in
the center of the bins) in the figure represent the statistical uncertainties
in the data.

We now present numerical results of jet production in PbPb collisions with jet
quenching. The event selection criteria for both $pp$ and PbPb collisions at
$\sqrt{s_{\mathrm{NN}}}=5.02$ TeV are the same, given in Sec. II. In order to
study the hot nuclear alteration of the jet broadening distribution, we plot
Fig. 3, which shows the normalized distributions of the jet broadening
${B_{tot}}$ in central 0-30% PbPb and $pp$ collisions at
$\sqrt{s_{\mathrm{NN}}}=5.02$ TeV, and their ratio in the lower panel. We find
wide distributions in $pp$ and PbPb collisions within the range given in the
figure. Most of the contribution of the distribution in $pp$ collisions lies
in the region with large $\ln(B_{tot})$ values. Moreover, the normalized
distribution in PbPb collisions is shifted toward smaller $\ln(B_{tot})$,
thus, leading to an enhancement at small $\ln(B_{tot})$ region and a
suppression at large $\ln(B_{tot})$ region. This modification indicates that
the proportion of survived events with the narrow jet distribution will
increase, after including parton energy loss in the QGP. Namely, Fig. 3 shows
that the nuclear modification reduces the multijet energy flow broadening, and
thus increases the imbalance of the energy flow.

In addition, we show that the jet broadening distribution in PbPb collisions
has a peak at around the ${B_{tot}}=0.082$ ($\ln(B_{tot})=-2.5$). We also
notice considerable enhancement at smaller $\ln(B_{tot})$ region in PbPb
collisions versus $pp$ collisions. Such an enhancement is much more pronounced
than some other event-shape observables, such as transverse sphericity. In the
following we will conduct a thorough investigation on the underlying reasons
for these medium modifications of jet broadening observed in Fig. 3.

Figure 4: Relative contribution fraction of total event number jet broadening
distribution of ${N_{\mathrm{jet}}}\ignorespaces=3$,
${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events in $pp$ and PbPb
collisions.

We separate the selected events into two categories: the number of jets in the
event equals to three (${N_{\mathrm{jet}}}\ignorespaces=3$) and more than
three (${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$). In Fig. 4, we plot the
production fractions as functions of $\ln(B_{tot})$ for three-jet events and
${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events respectively in $pp$ and
PbPb collisions. We find that in both $pp$ and PbPb collisions, the three-jet
events will give dominant contributions at the low $\ln(B_{tot})$ region,
whereas ${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events domonate at the
large $\ln(B_{tot})$ region .

We plot the event normalized $\ln(B_{tot})$ distributions of
${N_{\mathrm{jet}}}\ignorespaces=3$ events and
${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events in $pp$ and PbPb
collisions at $\sqrt{s_{\mathrm{NN}}}=5.02~{}\mathrm{TeV}$, as shown in Fig.
5. We can find that the $\ln(B_{tot})$ distributions of three-jet events and
multijet (${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$) events are quite
different for both $pp$ and PbPb collisions. In $pp$ collisions, the
$\ln(B_{tot})$ distribution of three-jet events has a relatively flat peak
near $-1.8$ (${B_{tot}}=0.165$), while the $\ln(B_{tot})$ distribution of
${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events has a sharp peak near
$-1.0$ (${B_{tot}}=0.368$). For PbPb collisions, the peak of three-jet events
shifts to the left near $-2.5$ (${B_{tot}}=0.082$), which is quite different
from the distribution of $pp$ collisions. The peak of
${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events is still around $-1.0$
compared with that in $pp$ collisions, but it is flatter. In a word, the
distribution of jet broadening for most three-jet events is more concertrated
than that for ${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events. After jet
quenching, these distributions become narrower, whether they are three-jet
events or events with more than three jets.

Combining the numerical results in Fig. 4 and Fig. 5, we now see that at large
${B_{tot}}$ region, the nuclear modifications of jet broadening distributions
are mainly determined by the events with
${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$. Since in PbPb the jet broadening
distribution in ${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events at large
${B_{tot}}$ region is suppressed as compared to $pp$ (see the lower panel of
Fig. 5), it is understandable that one observes a reduction of total jet
broadening (for all events) at large ${B_{tot}}$ region in PbPb as illustrated
by Fig. 3. Following the same logic, because in PbPb the jet broadening
distribution in ${N_{\mathrm{jet}}}\ignorespaces=3$ events at small
${B_{tot}}$ region is enhanced relative to $pp$ (see the upper panel of Fig.
5), we then see the increasing of total jet broadening (for all events) at
small ${B_{tot}}$ region in PbPb collisions.

We may see the jet number reduction effect plays an important role for the
medium modifications of total jet broadening in PbPb collisions. Due to the
jet quenching effect, some jets will fall below the threshold of event
selection after energy loss, which decreases the number of jets in such
events, and then leads to the modification of event-shape observables Chen
_et al._ (2020). In Table. 1, we compare the jet number fraction in $pp$ and
PbPb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02~{}\mathrm{TeV}$. The
proportion of three-jet events in $pp$ and PbPb collisions at
$\sqrt{s_{\mathrm{NN}}}=5.02~{}\mathrm{TeV}$ is 76.98% and 81.17%,
respectively, and the rest are events with four or more jets. In PbPb
collisions, the production fraction of three-jet events increases and the
fraction of ${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events decreases,
which may make jet broadening distributions more imbalanced, because
${N_{\mathrm{jet}}}\ignorespaces=3$ events tend to give small ${B_{tot}}$
while ${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events tend to contribute
more in large ${B_{tot}}$ region. It is noted that the reduction of
${N_{\mathrm{jet}}}\ignorespaces\geqslant 5$ events (shown in Table. 1) also
underlies why the $\ln(B_{tot})$ distribution of
${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events in PbPb collisions of Fig.
5 (lower) is suppressed at large $\ln(B_{tot})$.

Table 1: The relative production fraction of ${N_{\mathrm{jet}}}\ignorespaces=3$ events, ${N_{\mathrm{jet}}}\ignorespaces=4$ and ${N_{\mathrm{jet}}}\ignorespaces\geqslant 5$ events in $pp$ and PbPb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02~{}\mathrm{TeV}$. The kinematic cuts is consistent with the previous text. |  $pp$ |  PbPb
---|---|---
${N_{\mathrm{jet}}}\ignorespaces=3$ |  76.98$\pm$0.30% |  81.17$\pm$0.39%
${N_{\mathrm{jet}}}\ignorespaces=4$ |  18.59$\pm$0.13% |  15.53$\pm$0.25%
${N_{\mathrm{jet}}}\ignorespaces\geqslant 5$ |  4.43$\pm$0.05% |  3.30$\pm$0.05%

$N_{\rm jet}=3$

$N_{\rm jet}\geqslant 4$

Figure 5: Upper (Lower): Normalized jet broadening distributions of
${N_{\mathrm{jet}}}\ignorespaces=3$ (${N_{\mathrm{jet}}}\ignorespaces\geqslant
4$) events in $pp$ and PbPb collisions at
$\sqrt{s_{\mathrm{NN}}}=5.02~{}\mathrm{TeV}$. The distributions are self-
normalized. The bottom panels of the two figures give the corresponding medium
modification factors.

However, is the jet number reduction effect the only reason for the media
modifications shown in Fig. 3? To make further exploration we next track the
origin of events that satisfy the selection criteria. Due to the energy loss
effect, the selected quenched (heavy-ion) events and their unquenched version
may come from different kinematics regions. To expose the influence of these
contributions from different origins on jet broadening, we consider pairing
selected PbPb events with $pp$ events. The specific operation is selecting
quenched (heavy-ion) events satisfying the above selection condition (given in
Sec. II) and then using the matching procedure to match the unquenched ($pp$)
event corresponding to each selected quenched event. When we use the LBT model
to perform partons evolution, the information output to the HepMC3 Buckley
_et al._ (2021) file includes the final-state partons (quenched) list and the
input partons (unquenched) list. Therefore, this matching procedure is simple
and only needs to extract the partons list of the input part in the selected
event (quenched). Similar selection method has been utilized in Zhang  _et
al._ (2022); Brewer  _et al._ (2022), and we emphasize that in this work we
apply the method in the event level instead of the jet level in Ref. Zhang
_et al._ (2022); Brewer  _et al._ (2022).

Table 2: Classification of the unquenched versions of the selected quenched events. Matched Condition |  Category
---|---
Quenched |  UnQuenched |
$p_{T}^{\mathrm{min~{}jet}}>30{\mathrm{~{}GeV}}$ $110<p_{T,1}<170{\mathrm{~{}GeV}}$ ${N_{\mathrm{jet}}}\ignorespaces\geqslant 3$ |  $p_{T}^{\mathrm{min~{}jet}}>30{\mathrm{~{}GeV}}$ $110<p_{T,1}<170{\mathrm{~{}GeV}}$ ${N_{\mathrm{jet}}}\ignorespaces\geqslant 3$ (same as Quenched) |  Survival
|  $p_{T}^{\mathrm{min~{}jet}}>30{\mathrm{~{}GeV}}$ $p_{T,1}>170{\mathrm{~{}GeV}}$ ${N_{\mathrm{jet}}}\ignorespaces\geqslant 3$ |  Falldown
|  Other contribution |  Restructured

Figure 6: Comparison between the jet broadening $B_{tot}$ distributions in
three categories of quenched events. Left: the distributions of the three
types are scaled by the same value, which is the inverse of the total weights
of these three distributions. Right: the distributions of the three types are
self-normalized.

We classify the events that pass the cut conditions after energy loss into
three categories based on the patterns of the events before jet quenching, and
details are given in Table. 2. The tag of ”Survival” indicates that both the
quenched and unquenched versions of those events match the same selection
criteria. Some unquenched events satisfy other conditions for selecting
events, except that the leading jet $p_{T,1}$ is greater than the maximum
limit ($170{\mathrm{~{}GeV}}$). These events fully satisfy the selection
conditions after jet quenching, we record such events as ”Falldown”. In
addition to the above two categories, other events are classified as
”Restructured”. Detailed analysis shows that most ”Restructured” events come
from such unquenched events, that is, the transverse momentum of the leading
jet is greater than 110 GeV, and the number of jets greater than 30 GeV is
less than three. Specifically, such events account for about 90.2% of
restructured events. For these events, the number of jets with transverse
momentum greater than 30 GeV increases after medium modification. For example,
when some unquenched events that do not satisfy the event selection conditions
pass through the QGP, some jets will split in the medium due to jet-medium
interaction, thus the number of jets of these events in the final-state of
PbPb increases, and these events then satisfy the selection criteria Zhang
_et al._ (2022).

As shown in Fig. 6, we plot the distributions of these three categories of
quenched events at 5.02 TeV. According to the left plot in Fig. 6, we observe
that the contributions of these three categories are quite different, with the
”Survival” part having the most significant contribution at 57.85%, and the
contributions of the ”Falldown” and the ”Restructured” are 27.53% and 14.63%,
respectively. We can find that the broad peak of the jet broadening
$\ln(B_{tot})$ distribution for both the ”Survival” and the ”Falldown” of PbPb
events, these broad peaks spread from $\ln(B_{tot})=-2.6$ to $-1.0$, and the
”Restructured” part has a sharp peak around $\ln(B_{tot})=-2.6$
(${B_{tot}}=0.074$). In fact, the jet broadening distributions of the
”Survival” and ”Falldown” parts of PbPb events are not significantly
different, nor are they very different from the distribution of $pp$ events.
It is worth noting that the $\ln(B_{tot})$ of the ”Restructured” part is
mainly distributed in regions $-2.9$ to $-2.2$. From Fig. 6, we can see that
the new contribution from the ”Restructured” events plays an important role
for the enhancement of total jet broadening distributions at small
$\ln(B_{tot})$.

##  IV Summary

In this work, we present the first theoretical result of the medium
modifications of jet broadening distribution due to the jet quenching effect
in heavy-ion collisions at large momentum transfer. The $pp$ baseline is
provided by Powheg$+$Pythia8, where Powheg is responsible for generating NLO
precision hard process events and Pythia8 executes procedures such as parton
shower. We use the linear Boltzmann transport (LBT) model for simulating
parton energy loss to perform the evolution of jets in the medium.

The jet broadening observable ${B_{tot}}$, which belongs to event-shape
variables, can be used to characterize the broadening degree of multijet
energy flow. The decrease in this variable indicates that the jets’
distribution within the event becomes narrower. We calculate the medium
modification factor as a function of jet broadening variable ${B_{tot}}$.
Compared to $pp$ references, the results show an enhancement at the small
$B_{tot}$ region but a suppression at large $B_{tot}$ in PbPb collisions,
which implies the medium modification leads to a narrower distribution of jet
broadening relative to $pp$ references.

We further explore the possible underlying reasons for the medium modification
of jet broadening distributions. We demonstrate that jet broadening
distributions for three-jet events and for
${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events are pretty different in
both $pp$ and PbPb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02~{}\mathrm{TeV}$.
Due to parton energy loss in the QGP, the relative production fraction of
${N_{\mathrm{jet}}}\ignorespaces\geqslant 4$ events may be reduced. This jet
number reduction effect then results in the decreasing of jet broadening
distributions in the large $B_{tot}$ region as well as an increasement in the
low $B_{tot}$ region, in high-energy nuclear collisions. Furthermore, jet-
medium interactions may give rise to a new contribution from the
”Restructured” events, which further increases jet broadening distributions at
small $B_{tot}$ in PbPb collisions.

#### Acknowledgments:

We are grateful to S Y Chen and Q Zhang for providing helpful suggestions. We
acknowledge support for our work from following open-source projects, such as
Powheg Nason (2004); Frixione  _et al._ (2007); Alioli  _et al._ (2010); Nason
(2004), Pythia8 Sjöstrand  _et al._ (2015), FastJet Cacciari  _et al._ (2012),
HepMC3 Buckley  _et al._ (2021), Rivet Bierlich  _et al._ (2020), and
Matplotlib Hunter (2007). This work is supported by the Guangdong Major
Project of Basic and Applied Basic Research No. 2020B030103008, and Natural
Science Foundation of China with Project Nos. 11935007 and 12035007.

##  Appendix A The symmetric multijet limit

The main text mentioned the maximum values of transverse thrust and jet
broadening observables with given restrictions (ignore the $p_{z}$ of the
final state jets). Here we provide detailed proofs, following the discussions
in Ref. Banfi  _et al._ (2010). For a symmetrical planar transverse event with
jet numbers $N$ of equivalent momenta, the component jets can be expressed in
the $(p_{T},\eta,\phi)$-coordinate system as $p_{i}=(p_{T},0,\frac{2\pi}{N}i)$
for $i=1,2,\cdots,N$.

###  A.1 Transverse thrust

For a symmetrical multijet event with $N$ jets of equivalent momentum, the
specific direction of the transverse thrust axis is unimportant, but what
important is the angle between it and the nearest jet beside it, which we
denote here as $\Delta\phi$. So the transverse thrust ${\tau_{\perp}}$ can be
expressed as

$\tau_{\perp}$ $\displaystyle\equiv$ $\displaystyle
1-\max_{\hat{n}}\frac{\sum_{i}\left|\vec{p}_{T,i}\cdot\hat{n}_{T}\right|}{\sum_{i}p_{T,i}}$
(8) $\displaystyle=$ $\displaystyle
1-\frac{p_{T}\sum_{i=1}^{N}\left|\cos\left(\frac{2\pi
i}{N}-\Delta\phi\right)\right|}{N\cdot p_{T}},$

where $N$ is the number of jets and $N=2\pi/\phi$. We can get the
${\tau_{\perp}}$ value for perfectly circular planar events
($N\to\infty,\phi\to 0$),

$\displaystyle\left.{\tau_{\perp}}\right|_{\max}$ $\displaystyle=$
$\displaystyle\lim_{\phi\to
0}1-\frac{\sum_{i=1}^{N}\left|\cos\left(i\phi-\Delta\phi\right)\right|}{2\pi/\phi}$
(9) $\displaystyle=$ $\displaystyle 1-\frac{2}{\pi},$

notice here that when $\phi\to 0$, must also have $\Delta\phi\to 0$.

###  A.2 Jet broadening

A keypoint is to divide between the upper and lower hemispheres of the event
for the calculation of jet broadening observable. Since the interface of two
hemispheres of the event is perpendicular to the transverse thrust axis.
Therefore, in order to maximize the sum of the absolute values of the
projections of all jets on the transverse thrust axis, it is necessary to
ensure that no jets exist on the interface.

When $N$ is an even number, we choose the center of the first and $N$-th jet
as the interface, from which we can determine that $\phi_{U}$ and $\phi_{L}$
as $\pi/N+\pi/2$ and $\pi/N-\pi/2$, respectively. Then, we can get

$\displaystyle B_{U}$ $\displaystyle=$
$\displaystyle\frac{1}{2Np_{T}}\sum_{i\in\mathcal{C}_{U}}p_{T}\sqrt{\left(\frac{2\pi
i}{N}-\frac{\pi}{N}-\frac{\pi}{2}\right)^{2}}$ (10) $\displaystyle=$
$\displaystyle\frac{1}{4N^{2}}\sum_{i=1}^{N/2}\left|4\pi i-2\pi-N\pi\right|.$

Since the lower side and the upper side are completely symmetrical, and
$\phi_{U}$ and $\phi_{L}$ are on the same line, $B_{L}$ and $B_{U}$ are
exactly the same, we have

$B_{tot}$ $\displaystyle=$ $\displaystyle
2B_{U}=\frac{1}{2N^{2}}\sum_{i=1}^{N/2}\left|4\pi i-2\pi-N\pi\right|$ (11)
$\displaystyle=$ $\displaystyle\frac{\pi}{N^{2}}\sum_{i=1}^{\lfloor
N/4\rfloor}\left(2+N-4i\right),$ $\displaystyle=$
$\displaystyle\frac{\pi}{N^{2}}\left(N-2\left\lfloor\frac{N}{4}\right\rfloor\right)\left\lfloor\frac{N}{4}\right\rfloor.$

When $N$ is an odd number, the unit direction vector of any jet can be chosen
as transverse thrust axis $\hat{n}_{T}$,

$\displaystyle\mbox{$B_{tot}$}=B_{U}+B_{L}$ $\displaystyle=$
$\displaystyle\frac{1}{2Np_{T}}\sum_{i\in\mathcal{C}_{U}}p_{T}\sqrt{\left(\frac{2\pi
i}{N}\right)^{2}}+\frac{1}{2Np_{T}}\sum_{i\in\mathcal{C}_{L}}p_{T}\sqrt{\left(\frac{2\pi
i}{N}-\frac{\pi}{N}\right)^{2}}$ (12) $\displaystyle=$
$\displaystyle\frac{2\pi}{N^{2}}\sum_{i=1}^{\lfloor(N-1)/4\rfloor}i+\frac{1}{N^{2}}\sum_{i=1}^{\lceil(N-1)/4\rceil}(2\pi
i-\pi)$ $\displaystyle=$
$\displaystyle\frac{\pi}{N^{2}}\left(\left\lceil\frac{N-1}{4}\right\rceil^{2}+\left\lfloor\frac{N-1}{4}\right\rfloor^{2}+\left\lfloor\frac{N-1}{4}\right\rfloor\right).$

For perfectly circular planar events ($N\to\infty$), both Eq. (11) and Eq.
(12) tend to $\pi/8$.

## References

  * Busza  _et al._ (2018) W. Busza, K. Rajagopal, and W. van der Schee, Ann. Rev. Nucl. Part. Sci. 68, 339 (2018), arXiv:1802.04801 [hep-ph] .
  * Yagi  _et al._ (2005) K. Yagi, T. Hatsuda, and Y. Miake,  _Quark-gluon plasma: From big bang to little bang_ , Vol. 23 (2005).
  * Bjorken (1982) J. D. Bjorken,  _Energy loss of energetic partons in quark-gluon plasma: possible extinction of high p T jets in hadron-hadron collisions_, Tech. Rep. (FERMILAB, Batavia, IL, 1982).
  * Gyulassy  _et al._ (2004) M. Gyulassy, I. Vitev, X.-N. Wang, and B.-W. Zhang, (2004), arXiv:nucl-th/0302077 .
  * Wang and Gyulassy (1992) X.-N. Wang and M. Gyulassy, Phys. Rev. Lett. 68, 1480 (1992).
  * Qin and Wang (2015) G.-Y. Qin and X.-N. Wang, Int. J. Mod. Phys. E 24, 1530014 (2015), arXiv:1511.00790 [hep-ph] .
  * Abelev  _et al._ (2010) B. I. Abelev _et al._ (STAR), Phys. Rev. C 82, 034909 (2010), arXiv:0912.1871 [nucl-ex] .
  * Chatrchyan  _et al._ (2011) S. Chatrchyan _et al._ (CMS), Phys. Rev. C 84, 024906 (2011), arXiv:1102.1957 [nucl-ex] .
  * Zhang  _et al._ (2004) B.-W. Zhang, E. Wang, and X.-N. Wang, Phys. Rev. Lett. 93, 072301 (2004), arXiv:nucl-th/0309040 .
  * Vitev and Gyulassy (2002) I. Vitev and M. Gyulassy, Phys. Rev. Lett. 89, 252301 (2002), arXiv:hep-ph/0209161 .
  * Vitev (2005) I. Vitev, Phys. Lett. B 630, 78 (2005), arXiv:hep-ph/0501255 .
  * Zhang  _et al._ (2007) H. Zhang, J. F. Owens, E. Wang, and X.-N. Wang, Phys. Rev. Lett. 98, 212301 (2007), arXiv:nucl-th/0701045 .
  * Zhang  _et al._ (2009) H. Zhang, J. F. Owens, E. Wang, and X.-N. Wang, Phys. Rev. Lett. 103, 032302 (2009), arXiv:0902.4000 [nucl-th] .
  * Xu and Ko (2011) J. Xu and C. M. Ko, Phys. Rev. C 83, 021903 (2011), arXiv:1011.3750 [nucl-th] .
  * Eyyubova  _et al._ (2015) G. Eyyubova, V. L. Korotkikh, I. P. Lokhtin, S. V. Petrushanko, A. M. Snigirev, L. V. Bravina, and E. E. Zabrodin, Phys. Rev. C 91, 064907 (2015), arXiv:1411.4487 [hep-ph] .
  * Dai  _et al._ (2015) W. Dai, X.-F. Chen, B.-W. Zhang, and E. Wang, Phys. Lett. B 750, 390 (2015), arXiv:1506.00838 [nucl-th] .
  * Liu  _et al._ (2016) Z.-Q. Liu, H. Zhang, B.-W. Zhang, and E. Wang, Eur. Phys. J. C 76, 20 (2016), arXiv:1506.02840 [nucl-th] .
  * Dai  _et al._ (2017) W. Dai, X.-F. Chen, B.-W. Zhang, H.-Z. Zhang, and E. Wang, Eur. Phys. J. C 77, 571 (2017), arXiv:1702.01614 [nucl-th] .
  * Cao  _et al._ (2018) S. Cao, T. Luo, G.-Y. Qin, and X.-N. Wang, Phys. Lett. B 777, 255 (2018), arXiv:1703.00822 [nucl-th] .
  * Albacete  _et al._ (2019) J. L. Albacete, G. Giacalone, C. Marquet, and M. Matas, Phys. Rev. D 99, 014002 (2019), arXiv:1805.05711 [hep-ph] .
  * Ma  _et al._ (2019) G.-Y. Ma, W. Dai, B.-W. Zhang, and E.-K. Wang, Eur. Phys. J. C 79, 518 (2019), arXiv:1812.02033 [nucl-th] .
  * Pandey and Singh (2023) S. Pandey and B. K. Singh, Eur. Phys. J. A 59, 6 (2023), arXiv:2210.08527 [hep-ph] .
  * Bass  _et al._ (2009) S. A. Bass, C. Gale, A. Majumder, C. Nonaka, G.-Y. Qin, T. Renk, and J. Ruppert, Phys. Rev. C 79, 024901 (2009), arXiv:0808.0908 [nucl-th] .
  * Vitev  _et al._ (2008) I. Vitev, S. Wicks, and B.-W. Zhang, JHEP 11, 093 (2008), arXiv:0810.2807 [hep-ph] .
  * Vitev and Zhang (2010) I. Vitev and B.-W. Zhang, Phys. Rev. Lett. 104, 132001 (2010), arXiv:0910.1090 [hep-ph] .
  * Majumder and Van Leeuwen (2011) A. Majumder and M. Van Leeuwen, Prog. Part. Nucl. Phys. 66, 41 (2011), arXiv:1002.2206 [hep-ph] .
  * He  _et al._ (2012) Y. He, I. Vitev, and B.-W. Zhang, Phys. Lett. B 713, 224 (2012), arXiv:1105.2566 [hep-ph] .
  * Apolinario  _et al._ (2013) L. Apolinario, N. Armesto, and L. Cunqueiro, JHEP 02, 022 (2013), arXiv:1211.1161 [hep-ph] .
  * Burke  _et al._ (2014) K. M. Burke _et al._ (JET), Phys. Rev. C 90, 014909 (2014), arXiv:1312.5003 [nucl-th] .
  * Senzel  _et al._ (2015) F. Senzel, O. Fochler, J. Uphoff, Z. Xu, and C. Greiner, J. Phys. G 42, 115104 (2015), arXiv:1309.1657 [hep-ph] .
  * Blaizot and Mehtar-Tani (2015) J.-P. Blaizot and Y. Mehtar-Tani, Int. J. Mod. Phys. E 24, 1530012 (2015), arXiv:1503.05958 [hep-ph] .
  * Casalderrey-Solana  _et al._ (2016) J. Casalderrey-Solana, D. C. Gulhan, J. G. Milhano, D. Pablos, and K. Rajagopal, JHEP 03, 053 (2016), arXiv:1508.00815 [hep-ph] .
  * Milhano and Zapp (2016) J. G. Milhano and K. C. Zapp, Eur. Phys. J. C 76, 288 (2016), arXiv:1512.08107 [hep-ph] .
  * Sirunyan  _et al._ (2018) A. M. Sirunyan _et al._ (CMS), Phys. Lett. B 785, 14 (2018), arXiv:1711.09738 [nucl-ex] .
  * Connors  _et al._ (2018) M. Connors, C. Nattrass, R. Reed, and S. Salur, Rev. Mod. Phys. 90, 025005 (2018), arXiv:1705.01974 [nucl-ex] .
  * Tachibana  _et al._ (2017) Y. Tachibana, N.-B. Chang, and G.-Y. Qin, Phys. Rev. C 95, 044909 (2017), arXiv:1701.07951 [nucl-th] .
  * Dai  _et al._ (2020) W. Dai, S. Wang, S.-L. Zhang, B.-W. Zhang, and E. Wang, Chin. Phys. C 44, 104105 (2020), arXiv:1806.06332 [nucl-th] .
  * Zhang  _et al._ (2019) S.-L. Zhang, T. Luo, X.-N. Wang, and B.-W. Zhang, Nucl. Phys. A 982, 599 (2019), arXiv:1807.06197 [nucl-th] .
  * Aaboud  _et al._ (2019) M. Aaboud _et al._ (ATLAS), Phys. Rev. Lett. 123, 042001 (2019), arXiv:1902.10007 [nucl-ex] .
  * Wang  _et al._ (2019) S. Wang, W. Dai, B.-W. Zhang, and E. Wang, Eur. Phys. J. C 79, 789 (2019), arXiv:1906.01499 [nucl-th] .
  * Cao and Wang (2021) S. Cao and X.-N. Wang, Rept. Prog. Phys. 84, 024301 (2021), arXiv:2002.04028 [hep-ph] .
  * Acharya  _et al._ (2022) S. Acharya _et al._ (ALICE), Phys. Rev. Lett. 128, 102001 (2022), arXiv:2107.12984 [nucl-ex] .
  * Cunqueiro and Sickles (2022) L. Cunqueiro and A. M. Sickles, Prog. Part. Nucl. Phys. 124, 103940 (2022), arXiv:2110.14490 [nucl-ex] .
  * Wang  _et al._ (2021) S. Wang, J.-W. Kang, W. Dai, B.-W. Zhang, and E. Wang, Eur. Phys. J. A 58, 135 (2021), arXiv:2107.12000 [nucl-th] .
  * Chen  _et al._ (2022) S.-Y. Chen, J. Yan, W. Dai, B.-W. Zhang, and E. Wang, Chin. Phys. C 46, 104102 (2022), arXiv:2204.01211 [hep-ph] .
  * Kumar  _et al._ (2023) A. Kumar _et al._ (JETSCAPE), Phys. Rev. C 107, 034911 (2023), arXiv:2204.01163 [hep-ph] .
  * Dai  _et al._ (2022) W. Dai, M.-Z. Li, B.-W. Zhang, and E. Wang, (2022), arXiv:2205.14668 [hep-ph] .
  * Wang  _et al._ (2023a) L. Wang, J.-W. Kang, Q. Zhang, S. Shen, W. Dai, B.-W. Zhang, and E. Wang, Chin. Phys. Lett. 40, 032101 (2023a), arXiv:2211.13674 [nucl-th] .
  * Wang  _et al._ (2023b) S. Wang, W. Dai, E. Wang, X.-N. Wang, and B.-W. Zhang, Symmetry 15, 727 (2023b), arXiv:2303.14660 [nucl-th] .
  * Chen  _et al._ (2020) S.-Y. Chen, W. Dai, S.-L. Zhang, Q. Zhang, and B.-W. Zhang, Eur. Phys. J. C 80, 865 (2020), arXiv:2005.02892 [hep-ph] .
  * Mallick  _et al._ (2020) N. Mallick, S. Tripathy, R. Sahoo, and A. Ortiz, (2020), arXiv:2001.06849 [hep-ph] .
  * Sjostrand  _et al._ (2006) T. Sjostrand, S. Mrenna, and P. Z. Skands, JHEP 05, 026 (2006), arXiv:hep-ph/0603175 .
  * Banfi  _et al._ (2010) A. Banfi, G. P. Salam, and G. Zanderighi, JHEP 06, 038 (2010), arXiv:1001.4082 [hep-ph] .
  * Aad  _et al._ (2012) G. Aad _et al._ (ATLAS), Eur. Phys. J. C 72, 2211 (2012), arXiv:1206.2135 [hep-ex] .
  * Khachatryan  _et al._ (2014) V. Khachatryan _et al._ (CMS), JHEP 10, 087 (2014), arXiv:1407.2856 [hep-ex] .
  * Aad  _et al._ (2016) G. Aad _et al._ (ATLAS), Eur. Phys. J. C 76, 375 (2016), arXiv:1602.08980 [hep-ex] .
  * Aad  _et al._ (2021) G. Aad _et al._ (ATLAS), JHEP 01, 188 (2021), [Erratum: JHEP 12, 053 (2021)], arXiv:2007.12600 [hep-ex] .
  * Hessler  _et al._ (2022) J. Hessler, D. Britzger, and S. Lee (H1), PoS EPS-HEP2021, 367 (2022), arXiv:2111.11364 [hep-ex] .
  * Abe  _et al._ (1991) F. Abe _et al._ (CDF), Phys. Rev. D 44, 601 (1991).
  * Catani  _et al._ (1992) S. Catani, G. Turnock, and B. R. Webber, Phys. Lett. B 295, 269 (1992).
  * Abreu  _et al._ (1997) P. Abreu _et al._ (DELPHI), Z. Phys. C 73, 229 (1997).
  * Biebel  _et al._ (1999) O. Biebel, P. A. Movilla Fernandez, and S. Bethke (JADE), Phys. Lett. B 459, 326 (1999), arXiv:hep-ex/9903009 .
  * Abdallah  _et al._ (2003) J. Abdallah _et al._ (DELPHI), Eur. Phys. J. C 29, 285 (2003), arXiv:hep-ex/0307048 .
  * Wang  _et al._ (2022) S.-Q. Wang, C.-Q. Luo, X.-G. Wu, J.-M. Shen, and L. Di Giustino, JHEP 09, 137 (2022), arXiv:2112.06212 [hep-ph] .
  * Abreu  _et al._ (1996) P. Abreu _et al._ (DELPHI), Z. Phys. C 73, 11 (1996).
  * Kundu  _et al._ (2019) S. K. Kundu, T. Sarkar, and M. Maity, Int. J. Mod. Phys. A 34, 1950219 (2019), arXiv:1906.08472 [hep-ph] .
  * Ghosh  _et al._ (2022) A. Ghosh, X. Ju, B. Nachman, and A. Siodmok, (2022), arXiv:2203.12660 [hep-ph] .
  * Mallick  _et al._ (2021) N. Mallick, R. Sahoo, S. Tripathy, and A. Ortiz, J. Phys. G 48, 045104 (2021), arXiv:2008.13616 [hep-ph] .
  * Prasad  _et al._ (2022) S. Prasad, N. Mallick, D. Behera, R. Sahoo, and S. Tripathy, Sci. Rep. 12, 3917 (2022), arXiv:2112.03892 [hep-ph] .
  * Nason (2004) P. Nason, JHEP 11, 040 (2004), arXiv:hep-ph/0409146 .
  * Frixione  _et al._ (2007) S. Frixione, P. Nason, and C. Oleari, JHEP 11, 070 (2007), arXiv:0709.2092 [hep-ph] .
  * Alioli  _et al._ (2010) S. Alioli, P. Nason, C. Oleari, and E. Re, JHEP 06, 043 (2010), arXiv:1002.2581 [hep-ph] .
  * Alioli  _et al._ (2011) S. Alioli, K. Hamilton, P. Nason, C. Oleari, and E. Re, JHEP 04, 081 (2011), arXiv:1012.3380 [hep-ph] .
  * Sjöstrand  _et al._ (2015) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, Comput. Phys. Commun. 191, 159 (2015), arXiv:1410.3012 [hep-ph] .
  * Cacciari  _et al._ (2008) M. Cacciari, G. P. Salam, and G. Soyez, JHEP 04, 063 (2008), arXiv:0802.1189 [hep-ph] .
  * Cacciari  _et al._ (2012) M. Cacciari, G. P. Salam, and G. Soyez, Eur. Phys. J. C 72, 1896 (2012), arXiv:1111.6097 [hep-ph] .
  * Wang and Zhu (2013) X.-N. Wang and Y. Zhu, Phys. Rev. Lett. 111, 062301 (2013), arXiv:1302.5874 [hep-ph] .
  * He  _et al._ (2015) Y. He, T. Luo, X.-N. Wang, and Y. Zhu, Phys. Rev. C 91, 054908 (2015), [Erratum: Phys.Rev.C 97, 019902 (2018)], arXiv:1503.03313 [nucl-th] .
  * Cao  _et al._ (2016) S. Cao, T. Luo, G.-Y. Qin, and X.-N. Wang, Phys. Rev. C 94, 014909 (2016), arXiv:1605.06447 [nucl-th] .
  * Luo  _et al._ (2018) T. Luo, S. Cao, Y. He, and X.-N. Wang, Phys. Lett. B 782, 707 (2018), arXiv:1803.06785 [hep-ph] .
  * Eichten  _et al._ (1984) E. Eichten, I. Hinchliffe, K. D. Lane, and C. Quigg, Rev. Mod. Phys. 56, 579 (1984), [Addendum: Rev.Mod.Phys. 58, 1065–1073 (1986)].
  * Combridge (1979) B. L. Combridge, Nucl. Phys. B 151, 429 (1979).
  * Auvinen  _et al._ (2010) J. Auvinen, K. J. Eskola, and T. Renk, Phys. Rev. C 82, 024906 (2010), arXiv:0912.2265 [hep-ph] .
  * Guo and Wang (2000) X.-f. Guo and X.-N. Wang, Phys. Rev. Lett. 85, 3591 (2000), arXiv:hep-ph/0005044 .
  * Zhang and Wang (2003) B.-W. Zhang and X.-N. Wang, Nucl. Phys. A 720, 429 (2003), arXiv:hep-ph/0301195 .
  * Majumder (2012) A. Majumder, Phys. Rev. D 85, 014023 (2012), arXiv:0912.2987 [nucl-th] .
  * Pang  _et al._ (2012) L. Pang, Q. Wang, and X.-N. Wang, Phys. Rev. C 86, 024911 (2012), arXiv:1205.5019 [nucl-th] .
  * Pang  _et al._ (2015) L.-G. Pang, Y. Hatta, X.-N. Wang, and B.-W. Xiao, Phys. Rev. D 91, 074027 (2015), arXiv:1411.7767 [hep-ph] .
  * Lin  _et al._ (2005) Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, Phys. Rev. C 72, 064901 (2005), arXiv:nucl-th/0411110 .
  * Buckley  _et al._ (2021) A. Buckley, P. Ilten, D. Konstantinov, L. Lönnblad, J. Monk, W. Pokorski, T. Przedzinski, and A. Verbytskyi, Comput. Phys. Commun. 260, 107310 (2021), arXiv:1912.08005 [hep-ph] .
  * Zhang  _et al._ (2022) S.-L. Zhang, M.-Q. Yang, and B.-W. Zhang, Eur. Phys. J. C 82, 414 (2022), arXiv:2105.04955 [hep-ph] .
  * Brewer  _et al._ (2022) J. Brewer, Q. Brodsky, and K. Rajagopal, JHEP 02, 175 (2022), arXiv:2110.13159 [hep-ph] .
  * Bierlich  _et al._ (2020) C. Bierlich _et al._ , SciPost Phys. 8, 026 (2020), arXiv:1912.05451 [hep-ph] .
  * Hunter (2007) J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).