# Nonlinear Modes as a Tool for Comparing the Mathematical Structure of
Dynamic Models of Soft Robots

Pietro Pustina1,4, Davide Calzolari2,3, Alin Albu-Schäffer2,3, Alessandro De
Luca1, Cosimo Della Santina3,4 The reserach of Pietro Pustina was supported by
the Sapienza SoftRobot project. The research of Alessandro De Luca was
supported by the PNRR MUR project PE0000013-FAIR. The research of Cosimo Della
Santina was supported by the Horizon Europe Program from Project EMERGE under
Grant 101070918. The research of Davide Calzolari and Alin Albu-Schäffer was
supported by the Advanced Grant M-Runners (ID: 835284) by the European
Research Council. 1Department of Computer, Control and Management Engineering,
Sapienza University of Rome, Rome, Italy<EMAIL_ADDRESS><EMAIL_ADDRESS>of Robotics and Mechatronics, German
Aerospace Center, 82234 Oberpfaffenhofen, Germany<EMAIL_ADDRESS><EMAIL_ADDRESS>of Informatics, Technical University of
Munich, 85748 Garching, Germany. 4Department of Cognitive Robotics, Delft
University of Technology, Delft, The Netherlands<EMAIL_ADDRESS>
###### Abstract

Continuum soft robots are nonlinear mechanical systems with theoretically
infinite degrees of freedom (DoFs) that exhibit complex behaviors. Achieving
motor intelligence under dynamic conditions necessitates the development of
control-oriented reduced-order models (ROMs), which employ as few DoFs as
possible while still accurately capturing the core characteristics of the
theoretically infinite-dimensional dynamics. However, there is no quantitative
way to measure if the ROM of a soft robot has succeeded in this task. In other
fields, like structural dynamics or flexible link robotics, linear normal
modes are routinely used to this end. Yet, this theory is not applicable to
soft robots due to their nonlinearities. In this work, we propose to use the
recent nonlinear extension in modal theory –called eigenmanifolds– as a means
to evaluate control-oriented models for soft robots and compare them. To
achieve this, we propose three similarity metrics relying on the projection of
the nonlinear modes of the system into a task space of interest. We use this
approach to compare quantitatively, for the first time, ROMs of increasing
order generated under the piecewise constant curvature (PCC) hypothesis with a
high-dimensional finite element (FE)-like model of a soft arm. Results show
that by increasing the order of the discretization, the eigenmanifolds of the
PCC model converge to those of the FE model.

##  I Introduction

Nonlinear modal analysis is a powerful tool to characterize the behavior of
highly nonlinear mechanical systems [1], such as aerospace structures [2],
offshore towers [3] and micro/nano-electromechanical systems [4]. Despite the
fact that its use has been historically for system analysis and design of
structures undergoing small deformations, recent works [5, 6, 7] have used
nonlinear modal theory in analyzing robots and in generating energy-efficient
motions that asymptotically require no control effort to be sustained.

Figure 1: Schematic of our approach using nonlinear modal analysis and
eigenmanifold theory to compare different control-oriented models for the same
soft robot. Models differ in assumptions and number of DoFs, making systematic
comparison challenging. In the figure, $\sim$ indicates a similar structure,
$\nsim$ indicates dissimilarity and a question mark denotes incomparability.

Nonlinear modal analysis has found extensive application in flexible
mechanical systems under the hypothesis of small deformations - especially in
model reduction in structural mechanics have been studied [8, 9, 10, 11, 12].
In [10], energy-frequency plots are used to compare ROMs of thin shells and
plates. Finally, [11, 12] use nonlinear modes to assess ROM convergence to FE
models of specific structures.

Continuum soft robots [13], due to their high degree of deformability and
inherent softness, are expected to exhibit a rich spectrum of modal behaviors.
However, none of these works are directly applicable to soft robots because
they rely on nonlinear extensions of modal analysis built on assumptions that
do not extend to soft robots and other multi-body systems. Thus, the use of
nonlinear modes in continuum soft robotics is still unexplored.

This work proposes the first application of nonlinear modal analysis to soft
robotics by employing the recently developed concept of eigenmanifolds [14] \-
which, contrary to classic theories, builds on hypotheses compatible with soft
robotics. More specifically, with this work, we tackle the challenge of
quantitatively comparing control-oriented reduced order models (ROMs) for
continuum soft robots. In this way, it is also possible to compare control-
oriented ROMs obtained from different discretization hypotheses, including
those of [15, 16, 17, 18]. The goal is to have a mathematical tool that can be
used to quantify if two models have a similar mathematical structure. This is
hard to do by inspection of the models because even their configuration spaces
are going to be incompatible. Instead, we show in this work that modal theory
can be used to this end, as illustrated in Figure 1.

We propose a two-stage strategy to achieve this objective. First, we operate a
modal characterization of all the models to be compared. Then, we project the
modes onto task space, where we introduce similarity measures to compare their
shape and frequency content. To summarize, the main contributions of the paper
are as follows:

  * •

eigenmanifolds are proposed as a means of comparing control-oriented ROMs of
soft robots derived from different discretization hypotheses;

  * •

we propose a new energy-based continuation algorithm for eigenmanifold
computation;

  * •

a systematic procedure is proposed for comparing eigenmanifolds of different
ROMs, possibly obtained with different discretization techniques;

  * •

we compute for the first time nonlinear modes of continuum soft robots.

###  I-A Notation

Vectors and matrices are denoted in bold, while scalars are denoted in
lowercase normal font. The $n\times n$ identity matrix is represented by
$\boldsymbol{I}_{n}$. For a symmetric matrix
$\boldsymbol{A}=\boldsymbol{A}^{T}$ with $\lambda_{\min}(\boldsymbol{A})>0$,
we write $\boldsymbol{A}>0$ for positive definiteness. To simplify the
notation, arguments of functions are omitted when clear from the context.
Furthermore, given a continuously differentiable function
$f(\boldsymbol{x}):\mathbb{R}^{h}\rightarrow\mathbb{R}$,
$\boldsymbol{\nabla}_{\boldsymbol{x}}{f}\in\mathbb{R}^{h}$ and
$\boldsymbol{\nabla}_{\boldsymbol{x}}^{2}{f}\in\mathbb{R}^{h\times h}$ denote,
respectively, the gradient and the Hessian of $f$. We express the soft robot’s
dynamics as111As discussed in [19], common and advanced modeling techniques
used in soft robotics all yield equations of this form. Actuation and
dissipation terms do not appear here because nonlinear modes are defined on
the undamped and unactuated model.

$\boldsymbol{M}({\boldsymbol{q}})\ddot{\boldsymbol{q}}+\boldsymbol{c}({\boldsymbol{q}},\dot{\boldsymbol{q}})+\boldsymbol{\nabla}_{{\boldsymbol{q}}}{V({\boldsymbol{q}})}=\boldsymbol{0},$
(1)

where
${\boldsymbol{q}},\dot{\boldsymbol{q}},\ddot{\boldsymbol{q}}\in\mathbb{R}^{n}$
denote the configuration variables and their time derivatives,
$\boldsymbol{M}({\boldsymbol{q}})\in\mathbb{R}^{n\times n}$ is the mass matrix
and $\boldsymbol{c}({\boldsymbol{q}},\dot{\boldsymbol{q}})\in\mathbb{R}^{n}$
models Coriolis and centrifugal terms. Furthermore, $V({\boldsymbol{q}})$
represents the potential energy, including, e.g., the effects of gravitational
and elastic forces. We also denote the system energy as

$E({\boldsymbol{q}},\dot{\boldsymbol{q}})=\frac{1}{2}\dot{\boldsymbol{q}}^{T}\boldsymbol{M}({\boldsymbol{q}})\dot{\boldsymbol{q}}+V({\boldsymbol{q}}),$
(2)

and the solution of (1) at time $t$ from the initial condition
$({\boldsymbol{q}}_{0},\dot{\boldsymbol{q}}_{0})$ as
${\boldsymbol{q}}(t,{\boldsymbol{q}}_{0},\dot{\boldsymbol{q}}_{0})$.

##  II Background: Nonlinear Modes

Nonlinear modal theory has been traditionally developed for mechanical systems
with constant inertia matrix [20], but it has recently been extended to multi-
body mechanical systems [14]. Since continuum soft robots fall into the latter
category, in the following, we briefly summarize the main concepts from [14],
which we refer the reader to for an exhaustive treatment of the topic.

Let ${\boldsymbol{q}}_{eq}$ be a stable equilibrium configuration of (1),
i.e.,

$\boldsymbol{\nabla}_{{\boldsymbol{q}}}{V({\boldsymbol{q}}_{eq})}=\boldsymbol{0},\quad\boldsymbol{K}({\boldsymbol{q}}_{eq})=\boldsymbol{\nabla}_{{\boldsymbol{q}}}^{2}{V({\boldsymbol{q}}_{eq})}>0,$

which is equivalent to saying that ${\boldsymbol{q}}_{eq}$ is a local
minimizer of $E$ (for $\dot{\boldsymbol{q}}=\boldsymbol{0}$). The
linearization of (1) at
$({\boldsymbol{q}},\dot{\boldsymbol{q}})=({\boldsymbol{q}}_{eq},\boldsymbol{0})$
yields the second-order dynamics

$\Delta\ddot{\boldsymbol{q}}+\boldsymbol{M}^{-1\\!}({\boldsymbol{q}}_{eq})\boldsymbol{K}({\boldsymbol{q}}_{eq})\Delta{\boldsymbol{q}}=\boldsymbol{0},$
(3)

where $\Delta{\boldsymbol{q}}$ denotes a small increment of the configuration
with respect to ${\boldsymbol{q}}_{eq}$. Each eigenvector
$\boldsymbol{c}_{i},\ i\in\\{1,\cdots,n\\}$ of
$\boldsymbol{M}^{-1\\!}({\boldsymbol{q}}_{eq})\boldsymbol{K}({\boldsymbol{q}}_{eq})$
generates a local two-dimensional eigenspace

$ES_{i}=\mathrm{span}\left\\{\left(\begin{array}[]{c}\boldsymbol{c}_{i}\\\
\boldsymbol{0}\end{array}\right),\left(\begin{array}[]{c}\boldsymbol{0}\\\
\boldsymbol{c}_{i}\end{array}\right)\right\\}$

with the property that when
$(\Delta{\boldsymbol{q}},\Delta\dot{\boldsymbol{q}})$ is initialized in
$ES_{i}$, the dynamics evolves according to the harmonic oscillator equations,
namely

$\Delta\ddot{\boldsymbol{q}}+\omega_{i}^{2}\Delta{\boldsymbol{q}}=\boldsymbol{0},$
(4)

where $\omega_{i}^{2}$ is the eigenvalue associated with $\boldsymbol{c}_{i}$.
From the above equation, it can be shown that the solutions starting from
$ES_{i}$ satisfy two important properties: (i) they are periodic with period
$T_{i}=2\pi/\omega_{i}$ and (ii) they are invariant, i.e., will remain in
$ES_{i}$ at all times. These trajectories are called linear modes (LMs).

When moving to the nonlinear setting, each eigenspace $ES_{i}$ bends into a
two-dimensional sub-manifold of the state space, called eigenmanifold and
denoted as $\mathfrak{M}_{i}$. In analogy with the linear case,
$\mathfrak{M}_{i}$ can be seen as a collection of invariant periodic
trajectories of (1), called hereinafter nonlinear modes (NMs). However, NMs
can have different periods, while all LMs have the same period. For the
purpose of this paper, we need a further notion, namely that of eigenmanifold
generator. Given $\mathfrak{M}_{i}$, its generator is defined as

$\mathcal{G}_{\mathfrak{M}_{i}}=\left\\{({\boldsymbol{q}},\dot{\boldsymbol{q}})\in\mathfrak{M}_{i}\,\,|\,\,\dot{\boldsymbol{q}}=\boldsymbol{0}\right\\}.$

In a nutshell, $\mathcal{G}_{\mathfrak{M}_{i}}$ is the subset of
$\mathfrak{M}_{i}$ containing all points with zero velocity. Each element of
$\mathcal{G}_{\mathfrak{M}_{i}}$ is an initial rest configuration associated
with a periodic motion that remains in $\mathfrak{M}_{i}$ at all times. The
elements of $\mathcal{G}_{\mathfrak{M}_{i}}$ can be seen as the nonlinear
counterparts of the eigenvector $\boldsymbol{c}_{i}$. Differently from the
linear case, where we can generate $ES_{i}$ with only one vector, the
nonlinear case needs an entire set of points, i.e.,
$\mathcal{G}_{\mathfrak{M}_{i}}$, to generate $\mathfrak{M}_{i}$. An important
property of $\mathcal{G}_{\mathfrak{M}_{i}}$ is that it can be parameterized
by the energy, implying that every NM can itself be parameterized by its
energy. To simplify the notation, we will henceforth drop the subscript $i$
and denote the NM at $E$ as $\mathcal{G}_{\mathfrak{M}}(E)$.

##  III Eigenmanifold computation

For the computation of NMs, analytical and numerical approaches have been
proposed in the literature, see [20] for a thorough review on possible
strategies.

Remarkably, following a two-stage procedure, we can easily extend numerical
approaches used for the computation of NMs also for that of eigenmanifolds.
First, we compute the generator of $\mathfrak{M}$. Then, we take each point of
$\mathcal{G}_{\mathfrak{M}}$ as an initial condition to simulate the dynamics
forward in time and compute the corresponding NM, thus obtaining a slice of
$\mathfrak{M}$.

In this work, for the computation of $\mathcal{G}_{\mathfrak{M}}$ we use a
modified version of the numerical continuation algorithm presented in [21].
Superscripts are used to denote the iteration index. Each generator is
computed using a predictor-corrector algorithm that starts from the direction
given by the eigenvector $\boldsymbol{c}_{i}$ of the linearized dynamics at
the equilibrium. During the first iteration only, the linearized dynamics (3)
is used to compute a point close to $\mathcal{G}_{\mathfrak{M}}$ and
${\boldsymbol{q}}_{eq}$. In particular, an initial configuration increment
$\Delta{\boldsymbol{q}}_{0}$ is computed such that the linearized energy takes
a small incremental step. In this way, the trajectory of (1) from
${\boldsymbol{q}}_{0}^{1}={\boldsymbol{q}}_{eq}+\Delta{\boldsymbol{q}}_{0}$
with zero velocity stays close to $ES$ and has approximately the same period
$T^{1}$ of that of $\Delta{\boldsymbol{q}}_{0}$. Then,
$(T^{1},{\boldsymbol{q}}_{0}^{1},\boldsymbol{0})$ is refined with a Newton-
Raphson iteration scheme that seeks a solution of the periodicity constraint

${\boldsymbol{r}}({\boldsymbol{q}}_{0},T)=\left(\begin{array}[]{c}{\boldsymbol{q}}_{0}-{\boldsymbol{q}}(T,{\boldsymbol{q}}_{0},\boldsymbol{0})\\\\[2.0pt]
-\dot{\boldsymbol{q}}(T,{\boldsymbol{q}}_{0},\boldsymbol{0})\end{array}\right)=\boldsymbol{0},$
(5)

in the unknowns
$({\boldsymbol{q}}_{0},\boldsymbol{0})\in\mathcal{G}_{\mathfrak{M}}$ and
$T\in(0,\infty)$. At step $j\geq 2$, a prediction for the current solution
$(T^{j},{\boldsymbol{q}}_{0}^{j},\boldsymbol{0})$ is computed along the
direction given by the tangent vector to the previous solution
$(T^{j-1},{\boldsymbol{q}}_{0}^{j-1},\boldsymbol{0})$. The prediction is then
corrected with a shooting procedure. In [21], the shooting procedure aims at
solving ${\boldsymbol{r}}({\boldsymbol{q}}_{0}^{j},T^{j})=\boldsymbol{0}$,
which is an over-constrained system of $2n$ equations in $n+1$ unknowns. Since
the predictor step requires fairly small energy increments to guarantee
convergence of the procedure, the computation of the generator can be quite
time consuming. To speed up the evaluation, we propose to constrain the energy
increment between consecutive iterations. In particular, during the correction
stage we solve the augmented system

$\left(\begin{array}[]{c}{\boldsymbol{r}}({\boldsymbol{q}}_{0}^{j},T^{j})\\\\[2.0pt]
E({\boldsymbol{q}}_{0}^{j},\boldsymbol{0})-E({\boldsymbol{q}}_{0}^{j-1},\boldsymbol{0})-\Delta
E\end{array}\right)=\boldsymbol{0},$ (6)

being $\Delta E>0$ a desired energy step. This simple modification has proven
to significantly reduce the computation time, without affecting stability of
the procedure. Since convergence of the scheme is guaranteed only in the
vicinity of $\mathcal{G}_{\mathfrak{M}}$, $\Delta E$ must be properly selected
to prevent numerical instability. To this end, a simple adaptive rule has been
adopted. In particular, $\Delta E$ is initialized to a reference value
$\Delta\bar{E}$ and a maximum number $N_{\max}$ of iterations is defined for
the correction. If convergence is not achieved within $N_{\max}$ iterations,
$\Delta E$ is halved and the correction repeated.

##  IV Metrics

In this work, we aim at comparing how the eigenmanifolds of a continuum soft
robot change with the order and type of the discretization adopted and thus
with the number DoFs. A direct comparison in configuration space is not
possible in general, unless the discretizations have the same number of DoFs.

Figure 2: Schematic representation of the eigenmanifold projection. Given
three control-oriented ROMs, their eigenmanifolds $\mathfrak{M}_{1}$,
$\mathfrak{M}_{2}$ and $\mathfrak{M}_{3}$ cannot be directly compared in
configuration space. To overcome this issue, we introduce a constant
dimensional task space $\boldsymbol{h}({\boldsymbol{q}},\boldsymbol{s})$ in
which $\mathfrak{M}_{1}$, $\mathfrak{M}_{2}$ and $\mathfrak{M}_{3}$ are
projected and here define similarity measures.

Even though two models have the same number of DoFs, their configuration
spaces may not be directly comparable. For example, consider the dynamic model
of a soft robot obtained using a piecewise constant curvature and an affine
approximation of the strain. Even if the discretization order is chosen such
that the models have the same number of DoFs, the configuration variables will
still belong to different spaces, making a direct comparison impossible.

We propose to overcome this issue by introducing a task space function
$\boldsymbol{h}({\boldsymbol{q}},\boldsymbol{s}):\mathbb{R}^{n}\times\mathcal{S}\rightarrow\mathcal{T}$,
where $\boldsymbol{s}\in\mathcal{S}\subseteq\mathbb{R}^{q}$ denotes a set of
hyper-parameters for $\mathcal{T}\subseteq\mathbb{R}^{m}$ and $m$ has always
the same dimension, independently of the dimension $n$ of ${\boldsymbol{q}}$.
Given $\boldsymbol{h}({\boldsymbol{q}},\boldsymbol{s})$, any eigenmanifold can
be projected into $\mathcal{T}$, resulting in a collection of periodic orbits.
Because all trajectories of an eigenmanifold are periodic by definition, their
projection into $\mathcal{T}$ is a periodic trajectory in $\mathcal{T}$ with
same period. Therefore, the representation of $\mathfrak{M}$ into
$\mathcal{T}$ is a sub-manifold of $\mathcal{T}$ with properties similar to
$\mathfrak{M}$. The approach is illustrated in Fig. 2.

For example, consider the direct kinematics of a slender soft arm. Regardless
of the discretization technique used, it must always be possible to
reconstruct the robot backbone (using a single hyper-parameter $s$). In other
words, we know the direct kinematics mapping

$\boldsymbol{h}({\boldsymbol{q}},s)=\boldsymbol{p}({\boldsymbol{q}},s)=\left(\begin{array}[]{c}p_{x}({\boldsymbol{q}},s)\\\
p_{y}({\boldsymbol{q}},s)\\\ p_{z}({\boldsymbol{q}},s)\\\
\end{array}\right):\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{R}^{3},$
(7)

where $\boldsymbol{p}({\boldsymbol{q}},s)$ is the position of the point at
distance $s\in[0,L]$ from the base, being $L$ the robot length at rest.

Because all the level sets of the energy

$L_{\bar{E}}=\left\\{({\boldsymbol{q}},\dot{\boldsymbol{q}})\in\mathfrak{M}\,\,|\,\,E({\boldsymbol{q}},\dot{\boldsymbol{q}})=\bar{E}\right\\}$

give raise to evolutions that trace the same path in configuration space and
differ only in a phase shift, a single point of $L_{\bar{E}}$ should be
projected. The generator points are good candidates for this purpose because
they are all characterized by an initial zero velocity, which allows to get
rid of any phase shift due to the initialization in $\mathfrak{M}$.

Given the task manifold, it is now possible to define similarity measures
between two eigenmanifolds $\mathfrak{M}_{1}$ and $\mathfrak{M}_{2}$. These
can depend on different properties that one is interested in comparing, such
as their shape and spectral content. A common measure to compare two curves
$A(t)$ and $B(t)$ belonging to a metric space $S$ is the Fréchet distance
[22],

$\phi(A,B)=\inf_{\alpha,\beta}\max_{u\in[0,1]}\left\\{d(A(\alpha(u)),B(\beta(u)))\right\\},$
(8)

where $d(\cdot,\cdot)$ is the distance function of $S$. In a nutshell, the
above metric looks at all possible reparameterizations of $A$ and $B$ and
takes the one that minimizes their maximum distance. Given
$\boldsymbol{h}({\boldsymbol{q}},\boldsymbol{s})$ and (8), we can define the
following measures

$\boldsymbol{f}(E,\boldsymbol{s})=\phi(\boldsymbol{h}(\mathcal{G}_{\mathfrak{M}_{1}}(E),\boldsymbol{s}),\boldsymbol{h}(\mathcal{G}_{\mathfrak{M}_{2}}(E),\boldsymbol{s})),$
(9)

and

$\boldsymbol{F}(E)=\int_{\mathcal{S}}\boldsymbol{f}(E,\boldsymbol{v})\mathrm{d}\boldsymbol{v},$
(10)

where the Fréchet distance is applied componentwise to the rows of
$\boldsymbol{h}({\boldsymbol{q}},\boldsymbol{s})$. In the following, we refer
to $\boldsymbol{f}(E,\boldsymbol{s})$ and $\boldsymbol{F}(E)$, respectively,
as the modal Fréchet distance and modal integral Fréchet distance. The measure
$\boldsymbol{f}$ allows us to compare two modes for a given value of $E$ and
$\boldsymbol{s}$. On the contrary, $\boldsymbol{F}$ averages this measure over
the entire space of hyper-parameters for
$\boldsymbol{h}({\boldsymbol{q}},\boldsymbol{s})$. It is also worth remarking
that other choices instead of $\phi$ could have been made, see for example
[22].

Consider again the case of a slender soft arm and take
$\boldsymbol{h}({\boldsymbol{q}},s)=\boldsymbol{p}({\boldsymbol{q}},s)$.
Suppose we are interested in quantifying how close two eigenmanifolds are. A
first estimate can be obtained by looking at the tip evolution only, i.e.,
$\boldsymbol{f}(E,L)$. However, $\boldsymbol{f}(E,L)$ provides only a partial
indication of their similarity, because $\boldsymbol{f}(E,L)=\boldsymbol{0}$
does not imply $\boldsymbol{f}(E,s)=0$ for any other point $s\neq L$. To
obtain a more complete measure of similarity, we introduce
$\boldsymbol{F}(E)$. Note that if $\boldsymbol{F}(E)=\boldsymbol{0}$, then
$\boldsymbol{f}(E,\boldsymbol{s})=\boldsymbol{0}$ for all
$\boldsymbol{s}\in\mathcal{S}$, which implies that the NMs produce the same
motion for all the system particles.

To compare the spectra, we can adopt a similar strategy. To this end, we
consider the minimum value of the magnitude-squared coherence defined as
follows

$\gamma(A,B)=\min_{f}\frac{|g_{AB}(f)|^{2}}{g_{AA(f)}g_{BB(f)}},$ (11)

where $g_{AB}$ is the Cross-spectral power density between $A$ and $B$, and
$g_{AA}$ and $g_{BB}$ the auto power spectral density of $A$ and $B$,
respectively. At a given frequency, the coherence measures the extent to which
$B$ can be estimated from $A$ through a linear process. Note that unlike
$\phi(A,B)$, $\gamma(A,B)$ is not a distance function, as it is not symmetric
and lower values indicate less correlation between the two signals, thus worst
similarity. By applying $\gamma$ to the components of $\boldsymbol{h}$ we
define the modal coherence as

$\boldsymbol{g}(E,\boldsymbol{s})=\gamma(\boldsymbol{h}(\mathcal{G}_{\mathfrak{M}_{1}}(E),\boldsymbol{s}),\boldsymbol{h}(\mathcal{G}_{\mathfrak{M}_{2}}(E),\boldsymbol{s})),$
(12)

and the modal integral coherence

$\boldsymbol{G}(E)=\int_{\mathcal{S}}\boldsymbol{g}(E,\boldsymbol{v})\mathrm{d}\boldsymbol{v}.$
(13)

Considerations analogous to those for $\boldsymbol{f}$ and $\boldsymbol{F}$
hold.

The integral measures defined above can also be integrated over the energy
domain to obtain an overall measure of the similarity between the
eigenmanifolds, i.e.,

$\boldsymbol{F}_{E}=\int_{E_{{\boldsymbol{q}}_{eq}}}^{E_{\max}}\boldsymbol{F}(E)\mathrm{d}E,$

and

$\boldsymbol{G}_{E}=\int_{E_{{\boldsymbol{q}}_{eq}}}^{E_{\max}}\boldsymbol{G}(E)\mathrm{d}E,$

where $E_{\max}$ is the maximum energy level considered during analysis. In
future works, we will explore the use of these measure as loss terms in a
learning algorithm for the computation of control-oriented ROMs.

##  V Simulation Results

In this section, we apply the above metrics to compute and compare the
eigenmanifolds of a slender continuum soft robot. The goal is to investigate
the differences between PCC ROMs obtained by increasing the resolution of the
strain and a FE rigid model having a large number of DoFs. For the PCC models,
we enhance the resolution by increasing the number of bodies per unit of
length at rest. The robot has a cylindrical shape with radius
$r=0.02\leavevmode\nobreak\ [$\mathrm{m}$]$, uniform mass density
$\rho=1062\leavevmode\nobreak\ [$\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$]$, and
rest length $L=0.4\leavevmode\nobreak\ [$\mathrm{m}$]$. Furthermore, we
consider a linear elastic material with Young modulus
$Y=0.66\leavevmode\nobreak\ [$\mathrm{MPa}$]$ and Poisson ratio $\nu=0.5$. The
arm is mounted with the base rotated with the gravitational field so that the
straight stress-free configuration is a stable equilibrium, and motion occurs
in the $(x,z)$ plane.

All the bodies of the PCC models share the same physical properties, except
for the rest length which is chosen as $L/n$, being $n$ the discretization
order. The maximum number of bodies we consider is $n_{\max}=5$. For each of
these models, we compute all their $n\leq n_{\max}$ eigenmanifolds, using the
procedure outlined in Section III. The rigid robot consists of ten cylindrical
bodies, each one with length $L/10$ and mass $m=\displaystyle\rho\pi
r^{2}(L/10)\leavevmode\nobreak\ [$\mathrm{kg}$]$. Elasticity is introduced by
assuming linear elastic revolute joints with uniform stiffness coefficients
$k=Y\pi r^{4}10/(4L)\leavevmode\nobreak\
[$\mathrm{N}\text{\,}\mathrm{m}\text{\,}{\mathrm{rad}}^{-1}$]$. Because the
PCC models have at most five bodies, we compute only the first five
eigenmanifolds of the rigid model, ordered in ascending order by oscillation
frequency. The eigenmanifold computation is perfomed with a desired energy
step of $\Delta E=0.05\leavevmode\nobreak\ $\mathrm{[}\mathrm{J}\mathrm{]}$$
until the energy reaches the maximum value $E_{\max}=1\leavevmode\nobreak\
[$\mathrm{J}$]$. In the following, we use the direct kinematics of the
backbone (7) as task space function. The integral measures $F(E)$ and $G(E)$
are approximated through a Gaussian quadrature rule with ten Gaussian points.

Figure 3 shows stroboscopic plots of the robots in the workspace for the first
three NMs at $E_{\max}$.

(a) Mode 1; $t=0T_{1}$

(b) $t=0.2T_{1}$

(c) $t=0.3T_{1}$

(d) $t=0.5T_{1}$

(e) $t=0.7T_{1}$

(f) $t=0.8T_{1}$

(g) $t=0.9T_{1}$

(h) Mode 2; $t=0T_{2}$

(i) $t=0.2T_{2}$

(j) $t=0.3T_{2}$

(k) $t=0.5T_{2}$

(l) $t=0.7T_{2}$

(m) $t=0.8T_{2}$

(n) $t=0.9T_{2}$

(o) Mode 3; $t=0T_{3}$

(p) $t=0.2T_{3}$

(q) $t=0.3T_{3}$

(r) $t=0.5T_{3}$

(s) $t=0.7T_{3}$

(t) $t=0.8T_{3}$

(u) $t=0.9T_{3}$

Figure 3: Stroboscopic plots of the robots in the workspace for the first
(a)–(g), second (h)–(n) and third (o)–(u) NMs at
$E_{\max}=1\leavevmode\nobreak\ [$\mathrm{J}$]$ and normalized time instants
with respect to the mode period. Each color corresponds to a different model.

(a) PCC $1$ DoF

(b) PCC $2$ DoFs

(c) PCC $3$ DoFs

(d) PCC $4$ DoFs

(e) PCC $5$ DoFs

(f) Rigid $10$ DoFs

Figure 4: Energy-frequency relationships for the PCC (a)–(e) and FE (f)
discretization, respectively. The models show a similar behavior only for the
first two modes.

Figure 4 depicts the energy-frequency plot. Only the first two eigenmanifolds
at lower frequency have trajectories with similar periods. These differences
are reasonably expected to show up also in the task space, at least for the
coherence. Figure 5 reports the energy evolution of the modal Fréchet distance
at the end-effector. All the PCC models achieve a good match with the FE
model. Note indeed that the Fréchet distance is always less than
$6\leavevmode\nobreak\ [$\mathrm{cm}$]$. Furthermore, the accuracy improves
significantly as more bodies are considered. For the $x$ component, the
distance between the PCC and FE models decreases with the order of
discretization, as expected. However, the fifth order PCC model does not
always show the best match with the rigid robot along the $z$ direction.
Nevertheless, the modal integral Fréchet distance in Fig. 6 shows that the
five DoFs PCC model better fits the rigid motion overall. These results
demonstrate that energy-frequency relation is not an exhaustive similarity
measure in general. Moving to the coherence, see Fig. 7, the five DoFs model
has good performance only for the first mode. The same conclusions can be
drawn also for the modal integral coherence in Fig. 8. Based on the overall
performance, we conclude that the PCC robot with four DoFs has modal
evolutions that overlap with those of the rigid arm at the tip more closely.

(a) Mode $1$, $f_{x}(E,L)$

(b) Mode $2$, $f_{x}(E,L)$

(c) Mode $3$, $f_{x}(E,L)$

(d) Mode $4$, $f_{x}(E,L)$

(e) Mode $5$, $f_{x}(E,L)$

(f) Mode $1$, $f_{z}(E,L)$

(g) Mode $2$, $f_{z}(E,L)$

(h) Mode $3$, $f_{z}(E,L)$

(i) Mode $4$, $f_{z}(E,L)$

(j) Mode $5$, $f_{z}(E,L)$

Figure 5: Energy evolution of the modal Fréchet distance $\boldsymbol{f}(E,L)$
for the $x$ (a)–(e) and $z$ (f)–(j) tip position.

(a) Mode $1$, $F_{x}(E)$

(b) Mode $2$, $F_{x}(E)$

(c) Mode $3$, $F_{x}(E)$

(d) Mode $4$, $F_{x}(E)$

(e) Mode $5$, $F_{x}(E)$

(f) Mode $1$, $F_{z}(E)$

(g) Mode $2$, $F_{z}(E)$

(h) Mode $3$, $F_{z}(E)$

(i) Mode $4$, $F_{z}(E)$

(j) Mode $5$, $F_{z}(E)$

Figure 6: Energy evolution of the modal integral Fréchet distance
$\boldsymbol{F}(E)$ for the $x$ (a)–(e) and $z$ (f)–(j) backbone position.

(a) Mode $1$, $g_{x}(E,L)$

(b) Mode $2$, $g_{x}(E,L)$

(c) Mode $3$, $g_{x}(E,L)$

(d) Mode $4$, $g_{x}(E,L)$

(e) Mode $5$, $g_{x}(E,L)$

(f) Mode $1$, $g_{x}(E,L)$

(g) Mode $2$, $g_{z}(E,L)$

(h) Mode $3$, $g_{z}(E,L)$

(i) Mode $4$, $g_{z}(E,L)$

(j) Mode $5$, $g_{z}(E,L)$

Figure 7: Energy evolution of the modal coherence $\boldsymbol{g}(E,L)$ for
the $x$ (a)–(e) and $z$ (f)–(j) tip position.

(a) Mode $1$, $G_{x}(E)$

(b) Mode $2$, $G_{x}(E)$

(c) Mode $3$, $G_{x}(E)$

(d) Mode $4$, $G_{x}(E)$

(e) Mode $5$, $G_{x}(E)$

(f) Mode $1$, $G_{x}(E)$

(g) Mode $2$, $G_{z}(E)$

(h) Mode $3$, $G_{z}(E)$

(i) Mode $4$, $G_{z}(E)$

(j) Mode $5$, $G_{z}(E)$

Figure 8: Energy evolution of the modal integral coherence $\boldsymbol{G}(E)$
for the $x$ (a)–(e) and $z$ (f)–(j) backbone position.

##  VI Conclusions and Future Work

This paper presents the first investigation of eigenmanifolds for control-
oriented reduced-order modeling of continuum soft robots. To improve
computational efficiency, we present a new continuation algorithm for mode
computation with an energy-based term. We propose a systematic approach to
compare the eigenmanifolds of dynamic models resulting from different
discretization hypotheses and number of degrees of freedom (DoFs). To this
end, we project the eigenmanifolds onto a task space of fixed dimension and
introduce therein different similarity measures. The proposed method is used
to compare piecewise constant curvature (PCC) models of increasingly finer
discretization with a finite element (FE) model. The results show that the
task space modal evolutions of the PCC and FE models become similar when the
order of the PCC discretization increases.

Future work will explore the use of the proposed similarity measures to obtain
novel control-oriented models with a reduced number of DoFs.

## References

  * [1] Y. Mikhlin and K. V. Avramov, “Nonlinear normal modes of vibrating mechanical systems: 10 years of progress,” _Applied Mechanics Reviews_ , pp. 1–57, 2023.
  * [2] G. Kerschen, M. Peeters, J.-C. Golinval, and C. Stéphan, “Nonlinear modal analysis of a full-scale aircraft,” _Journal of Aircraft_ , vol. 50, no. 5, pp. 1409–1419, 2013.
  * [3] E. Gavassoni, P. B. Gonçalves, and D. de Mesquita Roehl, “Nonlinear vibration modes of an offshore articulated tower,” _Ocean Engineering_ , vol. 109, pp. 226–242, 2015.
  * [4] M. Matheny, L. Villanueva, R. Karabalin, J. E. Sader, and M. Roukes, “Nonlinear mode-coupling in nanomechanical systems,” _Nano Letters_ , vol. 13, no. 4, pp. 1622–1626, 2013.
  * [5] F. Bjelonic, A. Sachtler, A. Albu-Schäffer, and C. Della Santina, “Experimental closed-loop excitation of nonlinear normal modes on an elastic industrial robot,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 1689–1696, 2022.
  * [6] A. Coelho, A. Albu-Schaeffer, A. Sachtler, H. Mishra, D. Bicego, C. Ott, and A. Franchi, “EigenMPC: An eigenmanifold-inspired model-predictive control framework for exciting efficient oscillations in mechanical systems,” in _61st IEEE Conference on Decision and Control_ , 2022, pp. 2437–2442.
  * [7] Y. P. Wotte, S. Dummer, N. Botteghi, C. Brune, S. Stramigioli, and F. Califano, “Discovering efficient periodic behaviors in mechanical systems via neural approximators,” _Optimal Control Applications and Methods_ , 2023.
  * [8] T. Simpson, G. Tsialiamanis, N. Dervilis, K. Worden, and E. Chatzi, “On the use of variational autoencoders for nonlinear modal analysis,” in _Nonlinear Structures & Systems, Volume 1_, M. Brake, L. Renson, R. Kuether, and P. Tiso, Eds., 2023, pp. 297–300.
  * [9] E. Pesheck, C. Pierre, and S. W. Shaw, “Modal reduction of a nonlinear rotating beam through nonlinear normal modes,” _Journal of Vibration and Acoustics_ , vol. 124, no. 2, pp. 229–236, 2002.
  * [10] C. Touzé, M. Amabili, and O. Thomas, “Reduced-order models for large-amplitude vibrations of shells including in-plane inertia,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 197, no. 21-24, pp. 2030–2045, 2008.
  * [11] R. J. Kuether, M. R. Brake, and M. S. Allen, “Evaluating convergence of reduced order models using nonlinear normal modes,” in _Model Validation and Uncertainty Quantification, Volume 3_ , H. Atamturktur, B. Moaveni, C. Papadimitriou, and T. Schoenherr, Eds., 2014, pp. 287–300.
  * [12] R. J. Kuether, B. J. Deaner, J. J. Hollkamp, and M. S. Allen, “Evaluation of geometrically nonlinear reduced-order models with nonlinear normal modes,” _AIAA Journal_ , vol. 53, no. 11, pp. 3273–3285, 2015.
  * [13] C. Della Santina, M. G. Catalano, and A. Bicchi, “Soft robots,” in _Encyclopedia of Robotics_ , M. Ang, O. Khatib, and B. Siciliano, Eds. Springer, 2021.
  * [14] A. Albu-Schäffer and C. Della Santina, “A review on nonlinear modes in conservative mechanical systems,” _Annual Reviews in Control_ , vol. 50, pp. 49–71, 2020.
  * [15] S. H. Sadati, S. E. Naghibi, I. D. Walker, K. Althoefer, and T. Nanayakkara, “Control space reduction and real-time accurate modeling of continuum manipulators using Ritz and Ritz–Galerkin methods,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 1, pp. 328–335, 2017.
  * [16] S. Grazioso, G. Di Gironimo, and B. Siciliano, “A geometrically exact model for soft continuum robots: The finite element deformation space formulation,” _Soft Robotics_ , vol. 6, no. 6, pp. 790–811, 2019.
  * [17] F. Renda, C. Armanini, V. Lebastard, F. Candelier, and F. Boyer, “A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 3, pp. 4006–4013, 2020.
  * [18] B. Caasenbrood, A. Pogromsky, and H. Nijmeijer, “Control-oriented models for hyperelastic soft robots through differential geometry of curves,” _Soft Robotics_ , vol. 10, no. 1, pp. 129–148, 2023.
  * [19] C. Della Santina, C. Duriez, and D. Rus, “Model-based control of soft robots: A survey of the state of the art and open challenges,” _IEEE Control Systems Magazine_ , vol. 43, no. 3, pp. 30–65, 2023.
  * [20] G. Kerschen, M. Peeters, J.-C. Golinval, and A. F. Vakakis, “Nonlinear normal modes, part I: A useful framework for the structural dynamicist,” _Mechanical Systems and Signal Processing_ , vol. 23, no. 1, pp. 170–194, 2009.
  * [21] M. Peeters, R. Viguié, G. Sérandour, G. Kerschen, and J.-C. Golinval, “Nonlinear normal modes, Part II: Toward a practical computation using numerical continuation techniques,” _Mechanical Systems and Signal Processing_ , vol. 23, no. 1, pp. 195–216, 2009.
  * [22] Y. Tao, A. Both, R. I. Silveira, K. Buchin, S. Sijben, R. S. Purves, P. Laube, D. Peng, K. Toohey, and M. Duckham, “A comparative analysis of trajectory similarity measures,” _GIScience & Remote Sensing_, vol. 58, no. 5, pp. 643–669, 2021.


# Determining the viewing angle from TeV light curve of GRB 221009A

Lin Zhou1,Yuan-Chuan Zou1,2
1Department of Astronomy, School of Physics, Huazhong University of Science
and Technology, Wuhan 430074, China
2Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing, 210023,
China
E-mail<EMAIL_ADDRESS>(Y-CZ)

(Accepted XXX. Received YYY; in original form ZZZ; Accepted XXX. Received YYY;
in original form ZZZ)

###### Abstract

Gamma-ray bursts (GRBs) are among the most powerful explosive events in the
universe. LHAASO recently observed the most luminous one: GRB 221009A, and
unveiled its TeV light curve. The light curve exhibits a distinct jet break at
around 670 seconds, enabling the derivation of the viewing angle based on the
smoothness of the jet break. We constructed two models with or without
considering the high-latitude radiation, where the viewing angle was treated
as a free parameter, to fit the TeV light curve. We obtained the viewing
angles being 9.4 $\times 10^{-4}$ radians and 5.9 $\times 10^{-3}$ radians,
respectively. These values closely resemble an on-axis scenario, given the
opening angle is 1.4 $\times 10^{-2}$ radians.

###### keywords:

gamma-ray burst: individual: GRB 221009A - stars: jets

††pubyear: 2015††pagerange: Determining the viewing angle from TeV light curve
of GRB 221009A–Determining the viewing angle from TeV light curve of GRB
221009A

##  1 Introduction

Gamma-ray bursts (GRBs) are one of the most powerful explosive events in the
universe, releasing extremely highly energetic gamma-ray radiation. For a long
time, GRBs have been a fascinating research topic for astronomers because
their origins and nature involve a series of complex physical processes.

On October 9, 2022, the brightest GRB (GRB 221009A) initially triggered at
13:16:59.988 UTC by the Fermi Gamma-Ray Burst Monitor (Fermi-GBM; Meegan et
al. (2009)) in the gamma-ray band (Veres et al., 2022; Lesage et al., 2022).
This event was also detected (Pillera et al., 2022) by the Fermi Large Area
Telescope (Fermi-LAT; Atwood et al. (2009)). The Large High Altitude Air
Shower Observatory (LHAASO; Cao et al. (2019)) also reported the detection of
gamma rays, reaching up to 18 TeV (Huang et al., 2022). Swift Burst Alert
Telescope (Swift-BAT; Barthelmy et al. (2005)) also detected the same event
approximately an hour after the initial trigger. About 143 seconds later, the
Swift X-ray Telescope (Swift-XRT; Burrows et al. (2005)) pinpointed the
target, with the Swift Ultraviolet/Optical Telescope (UVOT; Roming et al.
(2005)), locating it at (RA, DEC) = (288.26452 ∘, 19.77350 ∘) with a 90%
confidence error radius of about 0.61 arcsec (Dichiara et al., 2022). The
event was quickly classified as a gamma-ray burst (Veres et al., 2022), and
according to X-shooter/VLT reports. it occurred at a redshift of 0.151 (Izzo
et al., 2022), confirming it as the brightest event ever recorded by any
gamma-ray burst detector. After two days, reports of a soft X-ray ring due to
dust scattering were released, based on observations from the Swift-XRT
(Tiengo et al., 2022). Multiwavelength observations were reported later(e.g.
Kennea et al., 2022; Bright et al., 2022).

The variation in the light curve and the corresponding physical processes of
GRB 221009A have become a highly discussed research topic. The half-opening
angle and the viewing angle are elements that influence the light curve during
the afterglow decay phase. The half-opening angle $\theta_{\rm jet}$ is the
angle at which the cone-shaped region extends from the central axis of the jet
towards one side. The viewing angle $\theta_{\rm obs}$, also known as the
observer’s angle, represents the angular deviation of the observer from the
cone axis. Both lead to a steepening decay in the late afterglow.

GRB 221009A is believed to have a narrow opening angle and a slightly off-axis
direction, as evidenced by the jet break in the later stages of the afterglow
light curve.

Quite amount of studies have been devoted to the investigation of the jet
opening angles of this Gamma-Ray Burst (e.g. Ren et al., 2023; Williams et
al., 2023; Laskar et al., 2023). However, limited attention has been given to
the consideration of off-axis angles. LHAASO Collaboration et al. (2023)
estimated a half-opening angle $\theta_{\rm jet}$ to be approximately
$0.8^{\circ}$, equivalent to 0.140 radians. This determination was made using
the jet break time and applying their late afterglow model (e.g. Sari et al.,
1999), assuming a viewing angle of zero. O’Connor et al. (2023)’s findings
yielded $\theta_{\rm obs}\lesssim 0.016$ rad and $\theta_{\rm jet}\geq$ 0.4
rad with a structured jet model. This study discussed a situation about a
structured jet with extended wings and determined the viewing angle.

Gill & Granot (2023) built a semi-analytical model to analyze an adiabatic
spherical blast wave and extended the model to angular structured jets. They
found that the viewing angle was comparable to the core angle of the energy
profile, $\theta_{\rm obs}\sim\theta_{\rm c,\epsilon}$, which was determined
to be 0.02 rad. This was employed to explain the high fluence of the gamma-ray
emission and to interpret the multi-waveband light curve.

Polarimetry was used to infer the viewing angle (e.g. Ghisellini & Lazzati,
1999). Negro et al. (2023) took a set of parameters to give a polarization
consistent with the IXPE spectropolarimetric measurement and found it favors a
jet opening angle approximately $1.5^{\circ}$, and a viewing angle about 2/3
of the jet opening angle, which correspond to 0.026 rad and 0.017 rad,
respectively.

In this work, we focus on the TeV light curve around the jet break. We draw
Fig. 1 as a sketch. If the viewing angle is zero, one may expect a sharp jet
break once the emission of the jet edge theoretically reaches the observer. If
the viewing angle is non-zero, the different part of the edge is then
gradually reaches the observer, and consequently the sharp break is smoothed.
By fitting the light curve with a top-hat jet model, we aim to determine the
viewing angle. We present the details of the model and the results in section
2.1 and 2.2. Data selecting is shown in section 2.3 and the method of fitting
is presented in section 2.4. The results are shown in section 2.5. The
conclusions and discussions are shown in section 3.

##  2 Constraining the viewing angle

###  2.1 top-hat model without considering high-latitude emission

After the peak of GRB afterglow light curve, the energy flux follows a power-
law decay with an index $\alpha$, which can be influenced by the visible area
and the viewing angle $\theta_{\rm obs}$. With the viewing angle, the power-
law decay is smoothed. Following the study by LHAASO Collaboration et al.
(2023), which assumes no high-latitude radiation, we constructed a simple
geometry model based on the top-hat model. This model describes the visible
area, which varies with time and the viewing angle.

Figure 1: Illustrating the modification of fitting results for the GRB
afterglow jet break through incorporation of the viewing angle effect. The
green dashed line is the fitting without viewing angle ($\theta_{\rm obs}=0$)
while the red solid line accounts for the viewing angle effect. Green dashed
line represents jet break time.

At a distance $R$, as shown in Fig. 2, with time passing and Lorentz factor
$\Gamma$ decreasing, the visible area ($R\theta$, green part) expands. The
blue circle represents the extent of the jet, with the blue dot in the center
indicating the direction of the jet. The dashed red circle represents the
observable range for the observer as the Lorentz factor decreases. The green
region represents the portion of the jet that actually enters the observer’s
observable range, and the gray area represents the part of the jet temporarily
invisible to the observer. $t_{\rm in}$ and $t_{\rm out}$ represent the
moments when the observing area first and second tangentially touch the edge
of the jet, respectively.

This simple geometry model constructs a simple physical scenario, as
exemplified by Woosley (2001). Because of the relativistic beaming effect, the
radiation outside the angle $1/\Gamma$ is ignored, while the radiation inside
is considered as viewing angle independent. It is denoted as Model 1 or
geometric model. (A more realistic model by considering the precise frame
transformation is discussed in section 2.2, which is denoted as Model 2.) As
the Lorentz factor decreases, the visible area expands, taking on the form of
a circle. For $t<t_{\rm in}$, the visible area is in circular shape, as shown
in the left panel of Fig. 2. During this period, the decay index $\alpha$ of
the GRB afterglow remains constant. At $t=t_{\rm in}$, which indicates that
the edge of the jet first becomes visible to the observer, the afterglow
undergoes a steepening decay. The alteration in the decay rate between $t_{\rm
in}$ and $t_{\rm out}$ is exclusively determined by the viewing angle and
half-opening angle. The latter is proportional to the jet’s visible area
divided by the total visible area. After $t>t_{\rm out}$, indicating the
complete visibility of the entire jet to the observer, the decay rate
stabilizes. The difference comparing to the on-axis ($\theta_{\rm obs}=0$)
case is during $t_{\rm in}<t<t_{\rm out}$. Utilizing this simple geometry
model, we can derive a plausible value for $\theta_{\rm obs}$.

Figure 2: Diagram of the variation in jet area as observed by an observer.
Left: when $t=t_{\rm in}$, the observing area is first tangent to the jet
edge;  Mid: when $t_{\rm in}<t<t_{\rm out}$, the observing circle overlaps
with a part of the jet; Right: when $t=t_{\rm out}$, the observing area covers
the jet and is tangent to the edge for the second time.

As relativistic particles emit radiation, their emissions become concentrated
within a narrow cone due to the relativistic beaming effect. This cone has a
half angle of $1/\Gamma$, where $\Gamma$ represents the Lorentz factor of the
jet. The direction of the cone aligns with the motion of the jet. As a result
of the relativistic beaming effect, the observer can only detect a portion of
the jet that is located within an angular distance smaller than the half
angle. In other words, the observable region is confined to a cone-shaped area
with a finite opening angle. The half-angle of the narrow cone can be
described as

$\theta(t)\propto\frac{1}{\Gamma(t)}=\frac{C_{1}}{\Gamma(t)},$ (1)

where $C_{1}$ is a constant. The Lorentz factor $\Gamma(t)$ obeys:

$\Gamma(t)\propto t^{-\alpha}=C_{2}t^{-\alpha},$ (2)

where $C_{2}$ is another constant determined by the actual value of the
Lorentz factor. The power law index $\alpha$ depends on the profile of the
number density within the surrounding environment. When the circum-burst
number density remains uniformly distributed, the power law index is -3/8
(Sari et al., 1998).

Combing equations (1) and (2) derives

$\theta(t)=\cfrac{C_{1}}{C_{2}}\cdot t^{\alpha}.$ (3)

The parameters $C_{1}$ and $C_{2}$ are combined into the ratio $C_{1}/C_{2}$
in equation (3) to represent the $\theta(t)$. To determine the variation of
$\theta(t)$, we need to consider constants $C_{1}$ and $C_{2}$. In this case,
we assume that $C_{1}=1$ and $C_{2}=\Gamma_{\rm peak}\cdot t_{\rm
peak}^{3/8}$. With equation (3), we can estimate the time boundaries $t_{\rm
in}$ and $t_{\rm out}$. Once the viewing angle $\theta_{\rm obs}$ is provided,
the instances at which the decay rate begins and ceases its steepening is
determined. The time boundaries are determined by

$\displaystyle\theta_{\rm jet}-\theta_{\rm obs}$
$\displaystyle=\frac{1}{\Gamma(t_{\rm in})},$ (4) $\displaystyle\theta_{\rm
jet}+\theta_{\rm obs}$ $\displaystyle=\frac{1}{\Gamma(t_{\rm out})}.$

As the external shock decelerates by the swept-up ambient mass, its Lorentz
factor decreases with time as $\Gamma\propto t^{-3/8}$. It is reasonable to
assume the flux density $F_{\nu}(t)$ follows

$F_{\nu}(t)=S_{\rm ratio}\cdot F_{\rm c}t^{-\alpha_{1}},$ (5)

where $S_{\rm ratio}$ is a function proportional to the visible area of the
jet divided by the total visible area as predicted by theory, $\alpha_{1}$ is
the power law index and $F_{\rm c}$ is an arbitrary constant for fitting. The
$S_{\rm ratio}$ satisfies

$\displaystyle S_{\rm ratio}=\left\\{\begin{array}[]{ll}\cfrac{\pi
R^{2}\theta^{2}}{\pi R^{2}\theta^{2}},&t<t_{\rm in},\\\ \cfrac{S_{\rm
overlap}}{\pi R^{2}\theta^{2}},&t_{\rm in}<t<t_{\rm out},\\\ \cfrac{\pi
R^{2}\theta_{\rm jet}^{2}}{\pi R^{2}\theta^{2}},&t>t_{\rm
out}.\end{array}\right.$ (9)

$S_{\rm overlap}$ represents the area where the visible area predicted by
theory overlaps with the jet, as shown in the middle panel of Fig. 2. With the
law of cosines, $S_{\rm overlap}$ can be expressed by

$\displaystyle S_{\rm overlap}$
$\displaystyle=R^{2}\times\left[2\arccos\left(\cfrac{\theta_{\rm
obs}^{2}+\theta_{\rm jet}^{2}-\theta^{2}}{2\theta_{\rm obs}\theta_{\rm
jet}}\right)\cdot\theta_{\rm jet}^{2}\right.$ (10)
$\displaystyle+2\arccos\left(\cfrac{\theta_{\rm
obs}^{2}+\theta^{2}-\theta_{\rm jet}^{2}}{2\theta_{\rm
obs}\theta}\right)\cdot\theta^{2}$
$\displaystyle\left.-2\sqrt{p(p-\theta)(p-\theta_{\rm obs})(p-\theta_{\rm
jet})}\right],$

where $p=(\theta+\theta_{\rm obs}+\theta_{\rm jet})/2$.

By substituting $S_{\rm overlap}$ into equation (9), we derive an expression
for the energy flux based on the previously assumed variables. To simplify
calculations, we opt for the logarithmic representation,

$\log(F_{\nu}(t))=-\alpha_{1}\log(t)+b+\log(S_{\rm ratio}).$ (11)

This equation represents the simple model with five parameters:
$\\{\theta_{\rm obs},\theta_{\rm jet},C_{1}/C_{2},\alpha_{1},b\\}$, where
$b=\rm log(F_{\rm c})$, $C_{1}$ and $C_{2}$ are degenerated as one parameter.
By combining the TeV light curve data, we estimate the optimal viewing angle.

###  2.2 top-hat model considering high-latitude emission

Though the geometric model in section 2.1 is simple and straightforward, the
high-latitude emission may be considered as a comparing model to show whether
the simplification is reasonable or not.

Given that the findings of the LHAASO Collaboration et al. (2023) were based
on the assumption of excluding high-latitude radiation, the previously used
half-opening angle $\theta_{\rm jet}=0.8^{\circ}$ becomes irrelevant. Instead,
we should treat $\theta_{\rm jet}$ as a free parameter. Our approach now
involves computing the radiation emitted from various segments of the jet and
aggregating them to determine the overall energy flux. This requires a
detailed analysis of the radiation contributions from different parts of the
jet to refine our model and accurately represent the observational data.

The flux from different segments of the jet can be described by the following
formula (Granot et al., 2002):

$f_{\nu}(\theta,t)=a^{3}f_{\nu/a}(0,at),$ (12)

where $a$ satisfies:

$a=\frac{1-\beta}{1-\beta\cos(\theta)},$ (13)

$\theta$ represents the angular distance between a segment of the jet and the
observer, and $\beta$ is velocity.

The flux density of each segment of the jet obeys (e.g. Granot & Sari, 2002):

$f_{\nu}(0,t)\propto t^{-\alpha_{2}}\nu^{-\beta_{2}}$ (14)

where $\alpha_{2}$ and $\beta_{2}$ are parameters that describe how the flux
varies with time and frequency. By summing up the flux density of every
segment, we can rewrite the flux formula as follows:

$F_{\nu}=F_{\nu,\rm const}t^{-\alpha_{2}}\nu^{-\beta_{2}}\iint_{\rm
jet}a^{3+\beta_{2}-\alpha_{2}}{\rm d}S$ (15)

where $F_{\nu,\rm const}$ is a coefficient of the formula.

To compare with the LHAASO observation, which is the flux in the range of 0.3
to 5 TeV, we need to integrate the flux density over frequency to fluence,
i.e., $F=\int_{0.3{\rm TeV}}^{5{\rm TeV}}F_{\nu}{\rm d}(\nu/a)$. Then we have

$F=\int_{0.3{\rm TeV}}^{5{\rm TeV}}F_{\rm\nu,const}\nu^{-\beta_{2}}{\rm
d}\nu\cdot t^{-\alpha_{2}}\iint_{\rm jet}a^{2+\beta_{2}-\alpha_{2}}{\rm d}S,$
(16)

where $F_{\rm const}=\int_{0.3{\rm TeV}}^{5{\rm
TeV}}F_{\rm\nu,const}\nu^{-\beta_{2}}{\rm d}\nu$ is a parameter for fitting.

By summing up all the energy flux from different angles, we can describe the
variation of the flux. This model incorporates several parameters: $\\{F_{\rm
const},C_{1}/C_{2},\alpha_{2},\beta_{2},\theta_{\rm obs},\theta_{\rm jet}\\}$.
When these parameters are provided, the energy flux of the jet is uniquely
determined.

###  2.3 The data

LHAASO is a new generation multi-component facility located in Daocheng,
Sichuan province of China, at an altitude of 4410 meters. It is designed for
the detection of very high energy cosmic rays and gamma rays. The cosmic rays
are measured in the energy range between $10^{12}$ and $10^{18}$ eV, while a
gamma-ray telescope is operated in the energy range between $10^{11}$ and
$10^{15}$ eV. With its capability to measure simultaneously different shower
components (electrons, muons, and Cherenkov light), LHAASO enables the
exploration of the origin, acceleration, and propagation of CR through its
very high-energy spectrum and elemental composition.

GRB 221009A happened to occur within the field of view of LHAASO, at a zenith
angle of $28.1^{\circ}$. During the first 3000 seconds, more than 64,000
photons (above 0.2 TeV) were detected. With a VHE gamma-ray extensive air
shower detector, which consists of three interconnected detectors, the LHAASO
collaboration analyses the light curve of GRB 221009A and tries to find the
best model to explain this extremely intense burst (LHAASO Collaboration et
al., 2023). In their observation, the TeV photon flux started to increase a
few minutes after the GRB trigger and reached its peak approximately 10
seconds later. A decay phase followed, which became more rapid at around 650
seconds after the peak.

They made estimations based on the data collected by LHAASO, which covered an
energy range from 0.3 to 5 TeV. With the data collected, they confirmed that
the power law break was caused by a jet break. With a top-hat model, they
estimated a half-opening angle of approximately $0.8^{\circ}$ (LHAASO
Collaboration et al., 2023).

We use the TeV light curve data of LHAASO Collaboration et al. (2023) for the
analysis. The energy flux light curve of GRB 221009A exhibits four-segment
features: rapid rise, slow rise, slow decay and steep decay. To determine the
viewing angle, it suffices to utilize data exclusively from the portions of
the light curve associated with slow decay and steep decay. To avoid involving
the rising part of the light curve, we extract the relevant data commencing
approximately 276 seconds after the Fermi GBM trigger. This interval also
corresponds to around 50 seconds post the TeV trigger.

After a new dataset is done, we employ the curve_fit python package (Virtanen
et al., 2020) to determine the optimal initial conditions for the MCMC
procedure and start MCMC sampling. Since our focus is the viewing angle
instead of other parameters, our requirements are confined to the parameters
explicitly detailed in Table 1.

Parameter | $\Gamma_{0}$ | $t_{\rm peak}$ | $\theta_{\rm jet}$ | $t_{\rm break}$ | $\alpha_{2}$ | $\alpha_{3}$ | $p$
---|---|---|---|---|---|---|---
value | 560 | 18.0 s | $0.8^{\circ}$ | 670 s | -1.115 | -2.21 | 2.2
Table 1: Parameters we choose from LHAASO Collaboration et al. (2023).
$\Gamma_{0}$ is initial Lorentz factor, $t_{\rm peak}$ is the moment of flux
peak, $\theta_{\rm jet}$ is the half-opening angle, $t_{\rm break}$ is the
moment of the jet break. $\alpha_{2}$ represents the decay index when light
curve started decay, and the $\alpha_{3}$ is the decay index after the jet
break. Notice $\alpha_{2}$ and $\alpha_{3}$ here (taken from LHAASO
Collaboration et al. (2023)) is different from the definition used in our
work. $p$ is power law index of the accelerated electron’s distribution.

###  2.4 MCMC analysis

The two models have several parameters, which can be collected as two
collections: $\Theta_{1}$ and $\Theta_{2}$ as

$\Theta_{1}\equiv\\{\theta_{\rm obs},\theta_{\rm
jet},C_{1}/C_{2},\alpha_{1},b\\},$ (17)

and

$\Theta_{2}\equiv\\{F_{\rm const},C_{1}/C_{2},\alpha_{2},\beta_{2},\theta_{\rm
obs},\theta_{\rm jet}\\}.$ (18)

To fit the light curve of GRB 221009A to our model, we utilize a MCMC
approach, following Ryan et al. (2015). This algorithm generates samples from
the posterior probability distribution $p(\Theta|D)$, which represents the
probability distribution of the parameters $\Theta$ given the light curve
data. By employing this approach, we gain a better understanding of the
overall structure of the posterior distribution and obtain the best estimates
for the parameters involved in our model. Ultimately, this method has the
potential to provide information regarding the uncertainties associated with
the inferred values of individual parameters.

For convenience, we have opted to change some parameters in our model. First,
we introduce a new parameter: q, which represents the ratio of $\theta_{\rm
obs}/\theta_{\rm jet}$. Second, we set $C_{1}$ as 1 and $C_{2}$ as
$\Gamma_{\rm peak}t_{\rm peak}^{3/8}$, so that parameter $C_{1}/C_{2}$ can be
simplified as a single parameter $\Gamma_{\rm peak}$. Third, we simplify
$F_{\rm const}$ to $F_{0}$, which equals log($F_{\rm const}$).

To reduce the models’ dimension, for model 1, we fix three parameters:
$\theta_{\rm jet}$, $t_{\rm peak}$, p, which can be found in Table 1. The
value of $\theta_{\rm jet}$ is set to $0.8^{\circ}$, which is about 0.0140
rad. $t_{\rm peak}$ is set to 18 s so that the parameter $C_{1}/C_{2}$ is
determined as $\Gamma_{\rm peak}$ is provided. And $\beta_{2}$ is set as p-1,
which is about 1.2.

By fixing these parameters, we transform the parameter set from $\Theta$ to
$\Theta_{\rm fit}$, which will be directly utilized in the MCMC routine, which
become

$\Theta_{\rm fit1}\equiv\\{q,\Gamma_{\rm peak},b,\alpha_{1}\\},$ (19)

and

$\Theta_{\rm fit2}\equiv\\{F_{0},\Gamma_{\rm peak},\theta_{\rm
jet},q,\alpha_{2}\\}.$ (20)

The posterior $p(\Theta_{\rm fit}|D)$ is calculated from the likelihood
$p(D|\Theta_{\rm fit})$ and the prior $p(\Theta)$ via Bayes’ Theorem:

$p(\Theta_{\rm fit}|D)\propto p(\Theta_{\rm fit})p(D|\Theta_{\rm fit}).$ (21)

The proportionality constant in equation (21) serves the purpose of
normalization and does not affect the MCMC analysis. When considering the data
points $(t_{i},F_{i})$, which are assumed to be independent of each other and
have Gaussian uncertainties $\sigma_{i}$, we can represent the likelihood as a
conventional $\chi^{2}$ form:

$\displaystyle p(D|\Theta_{\rm fit})$ $\displaystyle\propto\rm
exp(-\frac{1}{2}\chi^{2}),$ (22) $\displaystyle\chi^{2}$
$\displaystyle=\sum_{i}\left(\cfrac{F_{i}-F_{\rm model}(t_{i},\Theta_{\rm
fit})}{\sigma_{i}}\right)^{2},$ (23)

where $F_{\rm model}(t_{i},\Theta_{\rm fit})$ is the flux calculated from our
model at the time $t_{i}$ with parameters $\Theta_{\rm fit}$.

To restrict and determine the parameters of the fit based on prior information
about the data and model, the prior distribution $p(\Theta_{\rm fit})$ is
employed. The prior is assumed to be uniformly distributed within specific
limits as provided in Table 2. With the foundational theory in place, we
proceed to initiate the MCMC analysis.

During the fitting, the high precision and abundance of data points before the
jet break ($\sim$ 670 s) causes high weighting, which makes the data after the
jet break unimportant. Consequently, a model with even no jet break could fit
the data as well. To mitigate this effect, we randomly select 12 data points
from the high-precision part (before 670 s), and keep all the data after the
break, to form a new dataset for fitting. To minimize the uncertainty
associated with random selection, we repeated this data sampling procedure 100
times.

The MCMC analysis is performed using the emcee (Foreman-Mackey et al., 2013)
python package. This algorithm employs an ensemble of “walkers” moving
simultaneously through parameter space instead of the standard single walker
approach. Once the target distribution is defined, an initial state is chosen
within the bounds of the parameters in $\Theta_{\rm fit}$, and the process of
run$\\_$mcmc is initiated. Sampling is performed for 3800 iterations after
burn-in, generating approximately 550,000 samples of $p(\Theta_{\rm fit}|D)$.
To visually present the outcomes of the MCMC sampling process, corner package
(Foreman-Mackey, 2016) is utilized to plot the diagnostics.

###  2.5 The result

We focus on the afterglow data during the slow decay and steep decay phases.
The light curve data is integrated in the energy range 0.3–5 TeV. To ensure
reliable fitting results, we exclude the rising phase by selecting data points
beyond the initial 50 seconds. After this procedure, MCMC simulation is
employed to fit our models.

Figure 3: The fitting result of the geometry model (model 1). The best-fit
values are highlighted in blue. The diagonals along the contours depict the
marginalized probabilities for each parameter. The off-diagonal contours
display the covariances between different parameter pairs. Figure 4:  The
fitting result of the model considering high-latitude radiation (model 2). The
best-fit values are highlighted in blue.

The MCMC analysis results in the generation of contours, as depicted in Fig. 3
and Fig. 4. In these plots, each element corresponds to a histogram
representing the distribution of an individual variable, while the
correlations between variables are displayed within the contours.

It becomes evident that the distribution of all four variables in Fig. 3
closely resembles a Gaussian distribution. However, as shown in Fig. 4, it is
notable that the variables do not follow a Gaussian distribution for model 2.
Still, they exhibit some form of concentration that indicating a set of best
values.

By selecting the MAP (maximum posterior probability), the best-fit values are
determined. These values, along with their uncertainties, are presented in
Table 2.

Figure 5: The fitting result of the light curve is presented with the best
fitting values. The blue data points are the late TeV afterglow light curve of
GRB 221009A, selected after the initial 50 seconds. The red solid curve is the
best fitting result of model 1. The green dashed curve is the best fitting
result of model 2 by fitting the sampled data as described in section 2.4.

Utilizing the determined best-fit values, we depict the fitting results in
Fig. 5. It displays the fitting results of the models and the TeV data after
50 seconds since the start of the TeV emission. The consideration of high
latitude radiation in the model fitting process is found to be less effective
than excluding it, as high latitude radiation itself contributes to a
smoothing effect on the curvature (e.g. Granot et al., 2002). The program
successfully converges to stable parameter values, resulting in relatively
small uncertainties. In our geometry model without considering high-latitude
radiation, the outcome of the best-fit analysis indicates that the ratio q
approximates 0.0673, while the viewing angle $\theta_{\rm obs}$ is estimated
to be around 9.40$\times 10^{-4}$ rad, with an uncertainty of approximately
4.7$\times 10^{-4}$ rad. This result aligns with the assumption posited by
(LHAASO Collaboration et al., 2023), suggesting that the jet under
consideration exhibits a highly focused core and is nearly aligned with our
line of sight. However, upon accounting for high-latitude radiation, the
results obtained through MCMC analysis appear to diverge from the previously
obtained results. The best value of $\theta_{\rm jet}$ and $q$ after
considering high-latitude radiation is 0.0169 rad and 0.3518, respectively.
The value of $\theta_{\rm obs}$ is about 5.96$\times 10^{-3}$ rad, which is
about 6 times greater than the previous results. However, in the context of
comparing this value to the findings of other researchers, both results can be
considered negligible or close to zero. This discrepancy between our models
highlights the significant impact of model choice on the determination of
$\theta_{\rm obs}$.

In Fig. 5, one can see the model 1 (the simplified geometric model) fits the
data even better than the model 2 (the model considering the high-latitude
emission as well). We checked the dashed line by extending it to a much longer
time, and it returns to the steeper slope same as the solid one. This
indicates that the jet break has been smoothed in a much wider time range. It
is reasonable that model 2 has a shallower light curve after the jet break
compared with model 1 because the high-latitude emission can smear out the
over simplified geometric model. However, it is a puzzle that the data matches
the simplified model better (the red solid curves cross most of the data after
break, while the green dashed curve does not). A structured jet cannot relieve
this issue, as a top-hat jet is the sharpest structure. It might be a jet
break overlapping with spectral frequency crossing, which makes the light
curve steeping. However, the spectral index after the break is not softened as
argued by LHAASO Collaboration et al. (2023).

Model 1 |  Median |  Bound |  1$\sigma$ uncertainty |  3$\sigma$ uncertainty
---|---|---|---|---
q |  0.0673 |  [0, 1] |  (0.020, 0.127) |  (0.000, 0.255)
$\Gamma_{\rm peak}$ |  261.8 |  (70,560) |  (255.5,270.1) |  (245.9,304.5)
b |  -3.340 |  (-4,-2.5) |  (-3.363,-3.317) |  (-3.406,-3.278)
$\alpha_{1}$ |  1.109 |  (1.0,1.5) |  (1.098, 1.120) |  (1.077, 1.138)
Model 2 |  Median |  Bound |  1$\sigma$ Quantity |  3$\sigma$ Quantity
q |  0.3518 |  [0, 1] |  (0.1041,0.7275) |  (0.0017,0.9961)
$\Gamma_{\rm peak}$ |  291.6 |  (70,560) |  (183.6,453.7) |  (111.7,557.8)
$\theta_{\rm jet}$ |  0.0169 |  (0.005,0.05) |  (0.0088,0.0334) |  (0.0051,0.0496)
$F_{0}$ |  1.6241 |  (30,35) |  (1.2184,2.0098) |  (1.0047,2.3414)
$\alpha_{2}$ |  1.6757 |  (1.11,2) |  (1.5549,1.7816) |  (1.3558,1.8891)
Table 2: Best-fit values of $\Theta_{\rm fit1}$ and $\Theta_{\rm fit2}$ are
listed. Bounds and their uncertainties are also shown. Model 1 only applies
simple geometry, which is described in section 2.1. Model 2 is to consider the
high-latitude radiation, which is described in section 2.2.

##  3 conclusion and discussion

We developed two models to determine the viewing angle of GRB 221009A by
fitting its TeV light curve measured by LHAASO, which exhibited a distinct jet
break feature in the TeV band. Both models treat the jet as having a top-hat
angular distribution and assume the Lorentz factor decays as $\Gamma\propto
t^{-3/8}$. Through the application of MCMC sampling, we derived optimal
viewing angles of 9.4$\times 10^{-4}$ rad and 5.9$\times 10^{-3}$ rad. Both
viewing angles are sufficiently small to be regarded as approaching zero, thus
providing robust evidence supporting the on-axis assumption.

Our method has certain limitations. Firstly, we did not consider the
structured jet model (either power law (Dai & Gou, 2001; Rossi et al., 2002),
or Gaussian (Zhang et al., 2004)). Given the complex nature of the
environments in the TeV band, which require further discussion, it was
challenging to construct a reliable model for fitting. The second issue is the
determination of $\Gamma_{\rm peak}$. $\Gamma_{\rm peak}$ is assumed to be
derived from the initial Lorentz factor $\Gamma_{0}$, a value already
established by LHAASO Collaboration et al. (2023). However, the complex
parameters of the environment hindered our ability to obtain a reliable
$\Gamma_{\rm peak}$. Therefore, we opted to abandon the computation of an
exact peak Lorentz factor in favor of fitting. In this context, the evolution
of the jet within its environment is not accounted for. The third limitation
in our study is that we solely focused on the TeV energy flux observed by
LHAASO, neglecting multi-wavelength fitting. While this approach imposes
certain restrictions on the scope of our study, it effectively shielded us
from encountering challenging scenarios, such as data absence during jet
breaks (e.g. Kann et al., 2023). As a result, we were able to draw a clear
conclusion that the ejecta emitting TeV photons possesses a very small off-
axis angle, essentially treatable as zero.

waveband | $\theta_{\rm jet}$ | $\theta_{\rm obs}$ |
---|---|---|---
radio to X-rays |  $\geq$ 0.4 |  $\lesssim$ 0.016 | O’Connor et al. (2023)
radio to X-rays |  $\approx$ 0.021 |  $\sim$ 0.02 | Gill & Granot (2023)
2-8 keV X-ray | 0.0262 | 0.0175 | Negro et al. (2023)
radio to TeV | 0.01 |  | Sato et al. (2023)
| 0.1 |  |
X-ray and optical | 0.0245 |  | Ren et al. (2023)
0.3-10 keV X-ray | 0.0349 |  | Williams et al. (2023)
X-ray and optical |  $<$0.02 |  | Levan et al. (2023)
Radio to GeV | 0.0286 |  | Laskar et al. (2023)
X-ray to $\gamma$-ray | 0.0124 |  | An et al. (2023)
X-ray to $\gamma$-ray | 0.0127 |  | Zheng et al. (2024)
TeV | 0.0105 |  | Zhang et al. (2024)
X-ray and TeV | 0.0012 |  | Derishev & Piran (2024)
radio to TeV | 0.0234 |  | Ren et al. (2024)
| 0.0051 |  |
TeV | 0.0140 | 0.0009 | Model 1 (this work)
| 0.0169 | 0.0060 | Model 2 (this work)
Table 3: References related to jet opening angles and off-axis angles are
listed. Blank spaces indicate cases where observing angles were not discussed.
All these results are calibrated in units of radians.

Noticing different authors have obtained different jet opening angle and
viewing angle for this GRB, we have listed these angles from other literature
in Table 3 for a comparison. It is apparent that our results are significantly
smaller compared to those reported by other researchers. This is very likely
because of the data was used in different bands. For example, the “jet break”
is at around 670 s for TeV band (LHAASO Collaboration et al., 2023), while the
“break” is 0.255 day for the ultraviolet band (Williams et al., 2023), 0.98
day for the optical band (Kann et al., 2023). There is also a trend that
opening angle is smaller, constrained by the lower bands (such as optical band
and radio band). It might be that the jet opening angle is expanding while it
penetrates in the environment medium, because of lateral expansion (Rhoads,
1999). It is also might be caused by the structured jet, e.g., two-component
jet. The different viewing angles may correspond to the two components of the
jet (Sato et al., 2023), i.e., the smaller opening angle component (with
higher Lorentz factor) corresponds to the TeV emission, while the wider one
corresponds to the lower band emission.

As illustrated in Fig. 5, the rate of decay of jet bending after introducing
the high-latitude radiation is noticeably slower compared to the scenario
where high-latitude radiation is not considered. This suggests that the model
incorporating high-latitude radiation may not be comprehensive, and there may
be other influences leading to variations in the rate of jet breaking. We
speculate that the jet emitting TeV photons may also have its own structure,
or there are some spectral evolution overlapping at the jet breaking.

## Acknowledgements

We thank the helpful discussions with A. M. Chen, Weihua Lei, Kai Wang,
DuanYuan Gao, and the hospitality of Yao’an station of Purple Mountain
Observatory. The language was refined by ChatGPT. This work is supported in
part by the National SKA Program of China (2022SKA0130100), and by the science
research grants from the China Manned Space Project with No. CMS-
CSST-2021-B11. The simulation was completed on the HPC Platform of Huazhong
University of Science and Technology.

## Data Availability

The data used are publicly available at LHAASO Collaboration et al. (2023)
(Figs. 3 and 4), and can be downloaded from
https://www.nhepsdc.cn/resource/astro/lhaaso/paper.Science2023.adg9328/.

## References

  * An et al. (2023) An Z.-H., et al., 2023, arXiv e-prints, p. arXiv:2303.01203
  * Atwood et al. (2009) Atwood W. B., et al., 2009, ApJ, 697, 1071
  * Barthelmy et al. (2005) Barthelmy S. D., et al., 2005, Space Sci. Rev., 120, 143
  * Bright et al. (2022) Bright J., Rhodes L., Fender R., Farah W., Pollak A., Siemion A., 2022, GRB Coordinates Network, 32653, 1
  * Burrows et al. (2005) Burrows D. N., et al., 2005, Space Sci. Rev., 120, 165
  * Cao et al. (2019) Cao Z., et al., 2019, arXiv e-prints, p. arXiv:1905.02773
  * Dai & Gou (2001) Dai Z. G., Gou L. J., 2001, ApJ, 552, 72
  * Derishev & Piran (2024) Derishev E., Piran T., 2024, MNRAS,
  * Dichiara et al. (2022) Dichiara S., Gropp J. D., Kennea J. A., Kuin N. P. M., Lien A. Y., Marshall F. E., Tohuvavohu A., Williams M. A., 2022, The Astronomer’s Telegram, 15650, 1
  * Foreman-Mackey (2016) Foreman-Mackey D., 2016, The Journal of Open Source Software, 1, 24
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Ghisellini & Lazzati (1999) Ghisellini G., Lazzati D., 1999, MNRAS, 309, L7
  * Gill & Granot (2023) Gill R., Granot J., 2023, MNRAS, 524, L78
  * Granot & Sari (2002) Granot J., Sari R., 2002, ApJ, 568, 820
  * Granot et al. (2002) Granot J., Panaitescu A., Kumar P., Woosley S. E., 2002, ApJ, 570, L61
  * Huang et al. (2022) Huang Y., Hu S., Chen S., Zha M., Liu C., Yao Z., Cao Z., Experiment T. L., 2022, GRB Coordinates Network, 32677, 1
  * Izzo et al. (2022) Izzo L., et al., 2022, GRB Coordinates Network, 32765, 1
  * Kann et al. (2023) Kann D. A., et al., 2023, ApJ, 948, L12
  * Kennea et al. (2022) Kennea J. A., et al., 2022, GRB Coordinates Network, 33002, 1
  * LHAASO Collaboration et al. (2023) LHAASO Collaboration et al., 2023, Science, 380, 1390
  * Laskar et al. (2023) Laskar T., et al., 2023, ApJ, 946, L23
  * Lesage et al. (2022) Lesage S., Veres P., Roberts O. J., Burns E., Bissaldi E., Fermi GBM Team 2022, GRB Coordinates Network, 32642, 1
  * Levan et al. (2023) Levan A. J., et al., 2023, ApJ, 946, L28
  * Meegan et al. (2009) Meegan C., et al., 2009, ApJ, 702, 791
  * Negro et al. (2023) Negro M., et al., 2023, ApJ, 946, L21
  * O’Connor et al. (2023) O’Connor B., et al., 2023, Science Advances, 9, eadi1405
  * Pillera et al. (2022) Pillera R., Bissaldi E., Omodei N., La Mura G., Longo F., Fermi-LAT team 2022, GRB Coordinates Network, 32658, 1
  * Ren et al. (2023) Ren J., Wang Y., Zhang L.-L., Dai Z.-G., 2023, ApJ, 947, 53
  * Ren et al. (2024) Ren J., Wang Y., Dai Z.-G., 2024, ApJ, 962, 115
  * Rhoads (1999) Rhoads J. E., 1999, ApJ, 525, 737
  * Roming et al. (2005) Roming P. W. A., et al., 2005, Space Sci. Rev., 120, 95
  * Rossi et al. (2002) Rossi E., Lazzati D., Rees M. J., 2002, MNRAS, 332, 945
  * Ryan et al. (2015) Ryan G., van Eerten H., MacFadyen A., Zhang B.-B., 2015, ApJ, 799, 3
  * Sari et al. (1998) Sari R., Piran T., Narayan R., 1998, ApJ, 497, L17
  * Sari et al. (1999) Sari R., Piran T., Halpern J. P., 1999, ApJ, 519, L17
  * Sato et al. (2023) Sato Y., Murase K., Ohira Y., Yamazaki R., 2023, MNRAS, 522, L56
  * Tiengo et al. (2022) Tiengo A., Pintore F., Mereghetti S., Salvaterra R., a larger Collaboration 2022, GRB Coordinates Network, 32680, 1
  * Veres et al. (2022) Veres P., Burns E., Bissaldi E., Lesage S., Roberts O., Fermi GBM Team 2022, GRB Coordinates Network, 32636, 1
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Williams et al. (2023) Williams M. A., et al., 2023, ApJ, 946, L24
  * Woosley (2001) Woosley S., 2001, Nature, 414, 853
  * Zhang et al. (2004) Zhang B., Dai X., Lloyd-Ronning N. M., Mészáros P., 2004, ApJ, 601, L119
  * Zhang et al. (2024) Zhang B., Wang X.-Y., Zheng J.-H., 2024, Journal of High Energy Astrophysics, 41, 42
  * Zheng et al. (2024) Zheng C., et al., 2024, ApJ, 962, L2