# Motion Planning under Uncertainty: Integrating Learning-Based Multi-Modal
Predictors into Branch Model Predictive Control

Mohamed-Khalil Bouzidi1,3, Bojan Derajic 2,3, Daniel Goehring1, Joerg
Reichardt3 1 Free Universitity of Berlin, Germany {firstname.lastname@fu-
berlin.de}2 Technical Universitity of Berlin, Germany {firstname.lastname@tu-
berlin.de}3 Continental<EMAIL_ADDRESS>work is
funded by the German Federal Ministry forEconomic Affairs and Climate Action
within the project ”nxtAIM”.

###### Abstract

In complex traffic environments, autonomous vehicles face multi-modal
uncertainty about other agents’ future behavior. To address this, recent
advancements in learning-based motion predictors output multi-modal
predictions. We present our novel framework that leverages Branch Model
Predictive Control(BMPC) to account for these predictions. The framework
includes an online scenario-selection process guided by topology and collision
risk criteria. This efficiently selects a minimal set of predictions,
rendering the BMPC real-time capable. Additionally, we introduce an adaptive
decision postponing strategy that delays the planner’s commitment to a single
scenario until the uncertainty is resolved. Our comprehensive evaluations in
traffic intersection and random highway merging scenarios demonstrate enhanced
comfort and safety through our method.

##  I Introduction

Autonomous vehicles (AVs) must anticipate and predict the behavior of other
traffic participants (TPs) to plan safe and comfortable trajectories. In
dynamic environments, however, predictions of how a traffic scene will evolve
are quite uncertain over the entire planning horizon. This uncertainty arises
from the unknown TP intentions but also from other factors such as sensor
noise and occlusion. In addition, the probability distributions for the
development of a traffic situation are often multi-modal, i.e., at the end of
the planning horizon, several quite different environmental states may be
similarly probable (s. Fig 1). Consequently, planning algorithms must find a
plan that may not be optimal for any single future scenario, but rather
represents a compromise against all possibilities of the future.

Recently, various approaches are proposed to address such cases. In [1, 2]
uncertainty-aware planners are presented where the problem is formulated as a
Partially Observable Markov Decision Process (POMDP). These approaches suffer
from the fact that solving POMPD with increasing time horizon and action space
is no longer real-time capable. For this reason, the action space and time
steps are only coarsely sampled, causing suboptimal solutions hindering their
suitability for dealing with highly dynamic environments. On the other hand,
there are some variations of Model Predictive Control (MPC), which also take
into account the uncertainty.

Figure 1: An illustration of the challenge of planning in traffic scenes with
uncertain multi-modal predictions. To address this, Branch MPC optimizes over
a scenario tree.

One approach is robust MPC [3, 4]. These methods consider all possible cases
or the worst case to ensures collision avoidance in all possible scenarios.
However, these approaches lead to conservative, overcautious behaviors,
possibly preventing the AV from moving forward at all. To reduce conservatism,
stochastic MPCs [5, 6] (or variants such as scenario-based MPC [7, 8]) are
developed to exclude very unlikely predictions and guarantee constraints
satisfaction with a certain probability. Still, this can be very conservative,
especially with multi-modal uncertainty. The main source of conservatism is
that the planner ignores future information gathering about the scene that the
AV can use, i.e., it does not take into account feedback in its predictions.
For instance, in the scene depicted in Fig. 1 an experienced driver would be
aware that in some timesteps, he will get more information about which mode
will occur (i.e., whether the TP will turn right or cross the street) and use
it for his motion plan. In contrast, the mentioned robust or stochastic MPC
plan to comply with both modes for the entire planning horizon. This concept
of future feedback is included in the so-called Branch MPC which can plan with
different predictions of the traffic by allowing different control
inputs/strategies for different scenario evolutions.

Branch MPC for Motion planning and Control task has already been implemented
in [9, 10, 11, 12, 13, 14]. However, the integration of TP’s predictions in
these works were restricted to rule-based predictors[9, 10] or model-based
predictors [12, 13]. These prediction methods lack accuracy and flexibility
due to rigidity or incompleteness of the underlying model or rulesets (e.g.,
the proposed rule-based approaches can only take one TP into account),
limiting the prediction quality in scenes with complex interaction of
heterogeneous agents. In contrast, learning-based methods excel due to their
ability to handle complex and uncertain nature of other TP’s behavior
extracting patterns from high-dimensional data and providing multi-modal
predictions. Consequently, recent works are predominantly focused on learning-
based predictors [15, 16, 17]. However, integrating these predictors poses a
challenge for which viable solutions have yet to be proposed. Branch MPC
struggles with scalability regarding the number of scenarios which does not
allow to use the raw output of these predictor. Therefore, we propose a novel
approach to select only relevant predictions for the scenario tree based on
topology concepts and collision risk. Additionally, we further leverage the
learning-based output to design a decision postponing strategy by adaptively
setting the branching time (s. fig 1). This parameter significantly impacts
planner performance but is typically treated as a tuning parameter selected
offline in previous works. The main contributions of our work are therefore
summarized as:

  * •

An uncertainty-aware Motion Planning framework integrating learning-based
predictor into Branch MPC

  * •

A scenario selection strategy based on topology concepts and collision risk

  * •

An adaptive decision postponing strategy to adapt the branching time based on
the predictions

The remainder of the paper is organized as follows: We introduce the employed
multi-modal Predictor and the Baseline MPC in Section II. In Section III, we
augment the Baseline MPC to the Branch MPC Framework. Section IV presents our
strategy to select and integrate the relevant predictions of the TPs into the
Branch MPC. Section V presents our adaptive decision postponing method. The
evaluation is provided in Section V. Section VI concludes the paper.

##  II Preliminaries

Consider a motion planning problem of an AV in an arbitrary traffic scene with
$O$ TPs $o\in\\{0,..,O-1\\}$ (such as pedestrians, cars, bicycles). The AV
lacks knowledge about the intentions and future trajectory of TPs. It relies
solely on current and past states obtained through its environment perception
and tracking module. Additionally, a map $M_{in}$ containing polylines with
the respective road attributes of each point (such as location and lane type)
is given. A reference path with road boundaries is then extracted from
$M_{in}$ possibly by a high-level route planner. This reference path denoted
as $\mathcal{P}_{ref}:[0,\theta_{max}]\to\mathbb{R}^{4}\times[0,2\pi]$ maps
the arclength $\theta$ to a tuple consisting of the reference pose and
distances to the left and right road boundaries
$\theta\mapsto(x_{ref}(\theta),y_{ref}(\theta),d_{lb}(\theta),d_{rb}(\theta),\psi_{ref}(\theta))$
where $\theta_{max}$ is the maximum arclength.

###  II-A Learning-based Multi-Modal Motion Predictor

To effectively operate in a scene with TPs, the AV must anticipate their
behavior using a prediction module. Due to the TP’s unknown intentions and
range of potential choices, this needs to be inferred probabilistically.
Consequently, most SotA learning-based predictors provide multi-modal
distributed predictions for different possible behaviors of TPs. Catalyzed by
the availability of several large-scale datasets of traffic scenarios, a wide
range of these approaches has been developed in the recent literature. One of
these predictors is Motion Transformer (MTR) [15].We employ the MTR, noting
that any other predictor capable of multi-modal predictions could serve as an
alternative. The MTR is built using a transformer-based context encoder and
motion decoder. The inputs to the model are the history of position, heading
angle, and velocity of all surrounding TPs over a specified time frame.
Additionally, it incorporates the road map $M_{in}$ including the lane segment
attributes to provide scene context. Based on these inputs, the MTR generates
multiple predictions regarding how TPs could potentially interact. The MTR’s
output is structured as a tuple of Gaussian Mixture Models (GMMs) for each
object’s future state
$\mathcal{N}^{o}_{k,m}(\boldsymbol{\mu}^{o}_{k,m},\boldsymbol{\Sigma}^{o}_{k,m},\pi^{o}_{m})$
at each time step within the prediction horizon $N$. Each component
$m\in\mathcal{M}$ of this mixture represents a predicted mode. Hence, the
predicted distribution of each TP for each timestep $k$ can be expressed as:

$p^{o}_{k}=\sum_{m\in\mathcal{M}}\pi^{o}_{m}\mathcal{N}^{o}_{k,m}(\boldsymbol{\mu}^{o}_{k,m},\boldsymbol{\Sigma}^{o}_{k,m})$
(1)

Conversely, the predicted trajectories of each object $o$ are extracted from
the GMM along with the respective probability $\pi^{o}_{m}$ of each mode $m$.
These consist of the predicted states for each timestep $k$, denoted as
$\boldsymbol{o}_{k,m}=[x^{o}_{k,m},y^{o}_{k,m},\psi^{o}_{k,m},v^{o}_{i,m}]^{\top}$,
including position, heading, and velocity.

###  II-B Nominal Model Predictive Contouring Control

Model Predictive Contouring Control (MPCC) as introduced in [18] does not take
into account multi-modal predictions of the TPs and only plans with one
(usually with the highest probability $\pi^{o}_{m}$) prediction. As it serves
as our baseline that we will extend into Branch MPCC in the next section, we
now briefly recall its approach.

For the state space model
$\dot{\boldsymbol{z}}(t)=f(\boldsymbol{z}(t),\boldsymbol{u}(t))$ we use the
kinematic bicycle model:

$\dot{\boldsymbol{z}}=\left[v\cos(\psi),\ v\sin(\psi),\
v\frac{\tan(\psi)}{l},\ a,\ j,\ \dot{\delta},\dot{\theta}\right]^{\top}$ (2)

where $\boldsymbol{z}=[x,y,\psi,v,a,\delta,\theta]^{\top}$ is the state vector
and $\boldsymbol{u}=[j,\dot{\delta},\dot{\theta}]^{\top}$ the control input
vector. The acceleration, steering angle, jerk, steering angle rate, virtual
speed projected on the path, and wheelbase are denoted as
$a,\delta,j,\dot{\delta},\dot{\theta},l$, respectively. Variables with a
subscript $k$ denote the timestep $k$ of the discretized variables defined in
eq. (2).

The MPCC maximizes path progress while minimizing lag error $\hat{e}^{l}_{k}$
and contouring error $\hat{e}^{c}_{k}$, approximated by:

$\displaystyle\left[\begin{array}[]{cr}\hat{e}^{c}_{k}\\\
\hat{e}^{l}_{k}\end{array}\right]=\left[\begin{array}[]{cr}\sin(\psi_{ref}\left(\theta_{k}\right))&-\cos(\psi_{ref}\left(\theta_{k}\right))\\\
-\cos(\psi_{ref}\left(\theta_{k}\right))&-\sin(\psi_{ref}\left(\theta_{k}\right))\end{array}\right]\Delta\boldsymbol{p}_{ref}$

where $\Delta\boldsymbol{p}_{ref}=[x-x_{ref}\left(\theta_{k}\right),\
y-y_{ref}\left(\theta_{k}\right)]^{\top}$. This allows us to define the
running cost:

$J_{k}=\left[\begin{array}[]{cr}\hat{e}^{c}_{k}\\\
\hat{e}^{l}_{k}\end{array}\right]^{\top}Q\left[\begin{array}[]{cr}\hat{e}^{c}_{k}\\\
\hat{e}^{l}_{k}\end{array}\right]-q_{v}v_{k}^{p}+\boldsymbol{u}_{k}^{\top}R\boldsymbol{u}_{k}+J^{p}_{k}$
(3)

where $Q,q_{v},R$ are the respective weights. Additionally, we include the
cost term $J^{p}_{k}$, employing the Potential Field method [19], to address
the uni-modal uncertainty (i.e. of one mode) inherent in the predicted TP
trajectories.

$\begin{split}\@add@centering\centering
J^{p}_{k}&=q_{ob}\sum_{i=0}^{O-1}\cdot\exp\left(-\left(\frac{\Delta
x^{o}_{k}}{l^{o}}\right)^{2}-\left(\frac{\Delta
y^{o}_{k}}{w^{o}}\right)^{2}\right)\\\
&+q_{lm}\sum_{l=0}^{L-1}\exp\left(-\left(\frac{d_{lm}^{l}-\hat{e}^{c}_{k}\left(\theta_{k}\right)}{\sigma}\right)^{2}\right)\end{split}$
(4)

where $\Delta x^{i}_{k},\Delta y^{i}_{k}$ are the distances to the respective
obstacle, $d_{lm}^{l}$ is the signed distance from the reference path to the
respective $L$ lane marker, $\sigma$ a scaling factor and $l^{i},w^{i}$ a
conservative estimation of the length and width of the obstacle considering
the uni-modal uncertainty and $q_{ob},q_{lm}$ are the respective weights. This
supplements the existing hard constraints for obstacle avoidance and road
boundaries.

For feasible vehicle trajectories, we add the actuator constraints
$j_{k}\in[j_{\text{min}},j_{\text{max}}]$,
$\dot{\delta}_{k}\in[\dot{\delta}_{\text{min}},\dot{\delta}_{\text{max}}]$ and
$\delta_{k}\in[\delta_{\text{min}},\delta_{\text{max}}]$ and limit
longitudinal and lateral acceleration based on [20] .

Figure 2: Our learning-aided Motion Planning and Control Framework using a
Multi-modal Predictor, Scenario Selection, Adaptive Decision Postponing and
Branch MPCC

##  III Branch Model Predictive Contouring Control

To account for multi-modal predictions, we augment the MPCC.

The Branch MPCC (BMPCC) constructs a scenario tree to represent the
uncertainty in the scene. Each branch of the tree corresponds to a distinct
prediction (mode) of the traffic scene evolution, i.e., we consider only a
single branching (cf. Fig. 1). This single branch is denoted as a scenario
consisting of the $N$ nodes from the root node to the leaf node. Each scenario
$s\in\mathcal{S}$ has its own set of states
$\boldsymbol{z}_{k,s}\in\mathbb{R}^{n_{z}}$ (i.e. the respective nodes) and
control inputs $\boldsymbol{u}_{k,s}\in\mathbb{R}^{n_{u}}$, where
$n_{z},n_{u},|S|$ are the number of states, control inputs and scenarios,
respectively. We defer our discussion of which subset $\mathcal{S}$ of modes
from $\mathcal{M}$ to include in the scenario tree to sec. IV. We define the
set $I$ to encapsulate all elements $k,s$ within the scenario tree. This gives
us a discrete-time formulation of the dynamic system using eq. (2) :

$\boldsymbol{z}_{k+1,s}=f(\boldsymbol{z}_{k,s},\boldsymbol{u}_{k,s})$ (5)

In formulating the optimization problem, we employ the running cost and
constraint from the nominal MPCC as defined in sec. II-B.

$\displaystyle\underset{\boldsymbol{z}_{k,s},\boldsymbol{u}_{k,s}\forall(k,s)\in
I}{\min}\
\sum_{s\in\mathcal{S}}\omega_{s}\sum_{k=0}^{N-1}J_{k}(\boldsymbol{z}_{k,s},\boldsymbol{u}_{k,s})$
(6a) $\displaystyle s.t.\
\boldsymbol{z}_{k+1,s}=f(\boldsymbol{z}_{k,s},\boldsymbol{u}_{k,s})\
\forall(k,s)\in I$ (6b) $\displaystyle\boldsymbol{z}_{0,s}=\boldsymbol{z}(0)\
\forall s\in\mathcal{S}$ (6c)
$\displaystyle\boldsymbol{z}_{k,s}\in\mathcal{Z}(\boldsymbol{o}_{k,s}),\
\boldsymbol{u}_{k,s}\in\mathcal{U}\ \forall(k,s)\in I$ (6d)
$\displaystyle\boldsymbol{u}_{k,i}=\boldsymbol{u}_{k,j}\
\forall(i,j)\in\mathcal{S}\ \forall\ k\in[0,b]$ (6e)

where $\mathcal{Z}$ is a set of state constraints, including the obstacle
constraints for each scenario, and $\mathcal{U}$ is the set of control input
constraints. The process for selecting scenarios and defining the cost weight
$\omega_{s}$ is detailed in Section IV.

Each branch’s states and control inputs adapt to the selected prediction
(mode) of the obstacles $o$ to satisfy the constraints and optimize the cost
of the respective scenario, only. However, the BMPCC includes the constraints
in eq. (6e) which are called non-anticipatory

######  Remark 1.

The non-anticipatory constraints force the control inputs of all scenarios to
equal each other in the first timesteps $[t_{0},t_{b}]$, i.e., the so-called
branching time. These are introduced to model that the controller cannot
_anticipate_ in the beginning which of the scenarios will occur.

This setup, termed passive information gathering [14], models the intuition
that more information on scene evolution will be available in future. Up to
time $t_{b}$, the AV follows a consensus plan and postpones commitment to a
single mode. After the branching time $t_{b}$ these constraints in eq. (6e)
are no longer active. We detail this postponing in sec. V.

##  IV Scenario Selection

This section describes the construction of the scenario tree. To achieve this,
all $|\mathcal{M}|$ predictions outputted by the neural network in sec. II-A
are considered. The objective is to limit the number of considered scenarios
$\mathcal{S}$ by a small subset of the predicted modes,
$\mathcal{S}\subseteq\mathcal{M}$, $|\mathcal{S}|\leq|\mathcal{M}|$, to make
the BMPCC tractable for real-time applications while still not excluding
scenarios that may restrict the robustness of the controller. For that,
initially, the predicted trajectories that have similar effects on the planner
are clustered in sec. IV-A. Then, the most risk-relevant prediction from each
of these clusters is selected and added as a scenario in the scenario tree in
sec. IV-B. To determine these similarities and relevance of predictions for
the AV, i.e., for both of our methods presented subsequently, the future
behavior of the ego vehicle needs to be known. Since the current motion plan
$\boldsymbol{\tau}^{*}_{t}$ is not available in advance, we make the following
assumption:

######  Assumption 1.

The optimal motion plan calculated in the last planning cycle
$\boldsymbol{\tau}^{*}_{t-1}$ is sufficiently close to the current optimal
motion plan $\boldsymbol{\tau}^{*}_{t}$ to judge the topological relationships
of the agents.

Figure 3: Illustration of the Scenario Selection procedure. Predictions
enclosed in grey ellipses belong to same cluster. Red ellipses mark scenarios
selected from each cluster.

###  IV-A Topology-based Clustering

For clustering the predictions, we examine the topology of individual
predicted trajectories. The underlying idea is based on the insight that
certain obstacle behaviors, even if not identical, lead to similar optimal AV
behavior, while others require fundamentally different responses. An example
of that is provided in Fig. 3 where the ego vehicle must either brake or
accelerate at the intersection to maintain a safe distance from the TP
depending on the TP´s behavior111In this example, stopping avoids collision in
both cases but is far from being optimal due to excess of conservatism.. Our
methodology relies on grouping predictions of TPs according to this criterion,
considering exactly one scenario per cluster $\mathcal{C}_{s}$ where
$\bigcup_{s\in\mathcal{S}}\mathcal{C}_{s}=\mathcal{M}$ i.e. leading to a
partition of the set $\mathcal{M}$. To find these clusters, we employ uniform
visibility deformation (UVD) to assess the similarity of predicted
trajectories.

######  Definition 1.

Two trajectories
$\boldsymbol{\tau}_{1}(\alpha)$,$\boldsymbol{\tau}_{2}(\alpha)$ parameterized
by $\alpha\in[0,1]$ and satisfying
$\boldsymbol{\tau}_{1}(0)=\boldsymbol{\tau}_{2}(0)$,
$\boldsymbol{\tau}_{1}(1)=\boldsymbol{\tau}_{2}(1)$ belong to the same
_uniform visibility deformation_ , if for all $\alpha$, the line connecting
$\boldsymbol{\tau}_{1}(\alpha)$,$\boldsymbol{\tau}_{2}(\alpha)$ is collision-
free.[4]

UVD is typically utilized to compare ego trajectories for global trajectory
planning, aiming to efficiently identify multiple topologically distinct
trajectories[4, 21]. We use this method not to compare ego trajectories but to
cluster TP trajectories (i.e., the ego serve as obstacle). For this purpose,
the final position condition
($\boldsymbol{\tau}_{1}(1)=\boldsymbol{\tau}_{2}(1)$) is irrelevant to us, so
we disregard it.

The precise procedure for clustering entails checking whether two predicted
trajectories belong to the same group by uniformly discretizing them and
selecting equidistant timesteps as the parametrization. Subsequently, we
connect corresponding points on the two trajectories, i.e., points at the same
timestep, and determine if these lines intersect with the ego vehicle motion
plan $\boldsymbol{\tau_{t-1}}$ (s. Fig. 2). Thus, they do not belong to the
same cluster if at least one connected line intersects with the ego. In the
case of multiple TPs, if it is found that two trajectories are not
topologically equivalent for at least one TP, the respective predictions of
all TPs are assigned to different clusters.

###  IV-B Ranking and Selection based on Collision Risk

The aim of building the scenario tree is to map various reactions of the AV
that cover possible behaviors of the TPs for safe and comfortable driving.
Therefore, we add one prediction from each of the formed clusters to the
scenario tree. The most relevant prediction that best represents the entire
cluster must be selected. Representation here means that if the selected
scenario is taken into account in the planning, i.e. to avoid collisions, the
collision with the other predictions in the respective cluster are implicitly
also likely avoided (s. Fig. 3). For this process, we introduce the criterion
of collision risk. We utilize the work of [22], who presents an efficient way
to calculate the collision event probability between two oriented rectangles
analytically.

######  Definition 2.

The probability density of a collision event between a TP
$\boldsymbol{o_{k,m}}$ and the ego $\boldsymbol{z_{k}}$ at a certain timestep
$t_{k}$ is called Collision Event Probability (CEP) density
$\frac{dC^{o}_{m}(t_{k})}{dt}$. The CEP $C^{o}_{m}(t_{i},t_{j})$, the
probability of a collision during a period of time $[t_{i},t_{j}]$, i.e.,
between the ego plan and TP prediction, is the integral of the CEP density
over time.

Leveraging the procedure detailed in [22], we estimate
$\frac{dC^{o}_{m}(t_{k})}{dt}$ based on a first-order Taylor approximation.
The formula depends on the TP’s pose, velocity, and the respective Gaussian
uncertainty at timestep $t_{k}$. This we can extract from the output of our
prediction in sec. II-A. Furthermore, it depends on the ego plan, i.e., the
pose and velocity of the AV at timestep $t_{k}$. The trajectory of the ego
vehicle is considered deterministic where we rely on Assumption 1.

To obtain the CEP, we accumulate the timesteps over the prediction horizon
with the sample time $\Delta t$:

$C^{o}_{m}(t_{0},t_{N-1})=\sum_{k=0}^{N}\frac{dC^{o}_{m}(t_{k})}{dt}\Delta t$
(7)

We calculate this collision risk for each prediction of each TP. This
collision risk is utilized to define a metric for ranking and selecting the
predictions. An exemplary decision function for this purpose may be defined as
follows:

$D^{o}_{m}=C^{o}_{m}(t_{0},t_{N})+\lambda\pi^{o}_{m}$ (8)

The choice of the decision function is motivated to avoid excluding very
likely predictions (i.e. high occurrence probability $\pi^{o}_{m}$) that have
a low collision probability with the AV. If one would only consider
trajectories with a high collision risk and exclude the very likely
trajectory, one would adopt an excessively conservative approach. Accordingly,
the tuning factor $\lambda$ may be adapted as trade-off between conservatism
and effectiveness. Based on this decision metric, for each cluster
$\mathcal{C}_{s}$, the mode with the maximum value is introduced as scenario
$s$ into the scenario tree.

$s=\underset{m\in\mathcal{C}_{s}}{\text{argmax}}\sum_{o=0}^{O-1}D_{m}^{o}$ (9)

Once the scenario tree is constructed, the weighting $\omega_{s}$ for each
scenario in the cost function in eq. (6a) can be determined. This weighting is
derived through the sum of the probabilities of the individual predictions per
cluster.

##  V Adaptive Decision Postponing

This section outlines how the branching time, introduced in Remark 1, is
determined within our framework. The underlying assumption of the BMPCC is
that after the branching time is reached the uncertainty (which of the
predictions is the correct one) resolves, and the AV knows which option in the
scenario tree it should take. This uncertainty depends on various factors of
the traffic scene and cannot be predetermined offline independently of the
current situation. Thus, this necessitates to adapt the branching time online
based on the observed situation and depending on when it is expected to have
sufficient information for a decision. A wrong selection of the branching time
can harm the AV performance. If set too short, decisions may be made
prematurely, potentially resulting in the selection of the wrong scenario due
to the remaining uncertainty. A longer branching time postpones the decision,
aiming to accommodate all scenarios, but this can result in overly
conservative behavior where an opportunity may be missed, e.g., to merge into
a promising gap.

To determine the branching time, we estimate after how many timesteps the AV
can confidently assign the TPs state measurement to one of the predicted
trajectories. For this purpose, we introduce the Bhattacharyya distance [23].

######  Definition 3.

The Bhattacharyya Distance $B$ is a measure to quantify the similarity between
two probability distributions $p(x)$ and $q(x)$,
$B=-\ln(\int(\sqrt{p(x)q(x)}dx))$.

We leverage this measure to assess for every timestep in the prediction
horizon how much two Gaussian distributions
$\mathcal{N}^{o}_{k,i},\mathcal{N}^{o}_{k,j},\ i,j\in\mathcal{S}$ from the GMM
in eq. (1) overlap (s. Fig. 4). For two multivariate Gaussian distributions,
the formula simplifies to:

$\begin{split}B(o,k,i,j)&=\frac{1}{8}\left(\boldsymbol{\mu}^{o}_{k,i}-\boldsymbol{\mu}^{o}_{k,j}\right)^{T}\left(\boldsymbol{\Sigma}^{o}_{k}\right)^{-1}\left(\boldsymbol{\mu}^{o}_{k,i}-\boldsymbol{\mu}^{o}_{k,j}\right)\\\
&\quad+\frac{1}{2}\ln\left(\frac{\operatorname{det}\boldsymbol{\Sigma}^{o}_{k}}{\sqrt{\operatorname{det}\boldsymbol{\Sigma}^{o}_{k,i}\operatorname{det}\boldsymbol{\Sigma}^{o}_{k,j}}}\right)\end{split}$
(10)

where
$\boldsymbol{\Sigma}^{o}_{k}=\frac{\boldsymbol{\Sigma}^{o}_{k,i}+\boldsymbol{\Sigma}^{o}_{k,j}}{2}$.
If a predefined threshold is exceeded $B(o,k,i,j)\geq B_{th}$ for a certain
timestep we set the branching time at that particular time $t_{b}$. If more
than one TP is in the scene, we use the previously calculated collision risk
to threshold initially which TPs are relevant for decision postponing. After
that, if there are more than two TPs left, we calculate the time $t^{o}_{b}$
where $B(o,k,i,j)\geq B_{th}$ for each of these TPs and simply select the
maximum branching time $t_{b}=\text{max}_{o}t_{b}^{o}$. Similarly, if
$|\mathcal{S}|>2$, we calculate the distinguishability for all combinations
and select the maximum branching time.

Figure 4: Procedure of the online estimation of the branching time based on
the Bhattacharya Distance

##  VI Performance Evaluation

Our evaluation is divided into qualitative and quantitative results.
Qualitatively, we analyze our methods in a traffic intersection. For
quantitative results, we use Monte Carlo analysis of random highway merging
scenarios. We implement our framework in CasADi using the IPOPT solver for the
BMPCC. We use a prediction horizon of $N=40$ and a sample time of $\Delta
t=100\ ms$ for both our framework and the baselines. We limit the number of
scenarios selected by our method to $S_{max}=2$ and show that it can perform
better than using BMPCC with up to 4 scenarios. If our approach outputs more
than $S_{max}$ clusters $\mathcal{C}_{s}$, we select the clusters that have
the highest relevance based on our metric in eq.(8).

Our first Baseline is the conventional MPCC (_CMPCC_) introduced in sec. II-B.
Second, we use Scenario-based MPCC (_SCMPCC_), which takes the 5 TP
predictions with highest probability $\pi^{m}_{o}$, i.e., only excluding very
unlikely predictions. The main difference to the BMPCC is that SCMPCC does not
consider future feedback. In other words, the SCMPCC equals the BMPCC only if
the branching time is set identical to the prediction horizon ($t_{b}=t_{N}$).
The following variants are BMPCCs where we ablate different parts of our
method. _BMPCCnoSS2_ , _BMPCCnoSS3_ and _BMPCCnoSS4_ leaves out the scenario
selection and decision postponing part and only add the $S_{max}=2$,
$S_{max}=3$, and $S_{max}=4$ scenarios with highest probability to the
scenario tree. Then, in BMPCCCnoDP, we employ the scenario selection with
$S_{max}=2$ but leave out the adaptive decision postponing. For all four of
these variants, the branching time is set to the very first timestep as done
in [12, 10, 11].

TABLE I: Comparison using Monte Carlo analysis of Highway Merging Scenarios | Merging Execution | Convergence Quality
---|---|---
|  Success |  Aborted |  Collision |  Mean Cost |  Mean Solving time (ms)
CMPCC (1 Scen.) |  89% |  5% |  6% |  1397.6 |  22
SCMPCC |  86% |  14% |  0% |  1226.5 |  73
BMPCC no Scenario Selection (2 Scen.) |  89 % |  8 % |  3 % |  1039.9 |  50
BMPCC no Scenario Selection (3 Scen.) |  89 % |  9 % |  2 % |  883.4 |  92
BMPCC no Scenario Selection (4 Scen.) |  90 % |  9 % |  1 % |  1074.0 |  167
BMPCC no adapt. Dec. Postponing (2 Scen.) |  89 % |  9 % |  2 % |  817.8 |  64 + 9 (BMPCC+ rest)
Our Framework (2 Scen.) |  94 % |  5% |  1% |  575.5 |  54 + 10 (BMPCC+ rest)

###  VI-A Traffic Intersection

Figure 5: Prediction of TP (blue car) and the plan of the AV (yellow car) at a
certain timestep. Predictor outputs 6 predictions (shown in different colors)
from which 2 most-relevant are selected. Corresponding trajectories are shown
in corresponding colors. Right: Planned velocity profile.

In this traffic scenario, both our AV and a TP approach an intersection where
the TP has the right of way. Our predictor predicts multiple trajectories for
the TP that either turns right before the intersection or crosses it (s. Fig.
5). Accordingly, as depicted in Fig. 5, our approach estimates that
uncertainty resolves after $t_{b}=0.8s$ and only slightly brakes to
accommodate both possible TP behaviors. We run the scenario twice, once where
the TP turns right before the intersection and once where it crosses. In Fig.
6, the performance of our method is compared to the other methods. Costs are
calculated for the entire scenario duration, following the same criteria
(comfort, reference tracking, etc.) as defined in eq. (3). When the TP crosses
the intersection, CMPCC and BMPCCnoSS2 perform poorly because the most
probable predictions indicate a right turn in the beginning. In the case of
CMPCC, emergency braking is even required to avoid a collision. Even when the
most-relevant scenarios are selected with our method, BMPCCnoDP performs
poorly because it always expects uncertainty to resolve in the next time step,
which is not the case. As a result, it starts braking later than our
framework, leading to significantly less comfortable behavior. SCMPCC shows
the best behavior for this scenario because it plans for the entire prediction
horizon for the crossing scenario. However, in the scenario where the TP turns
right, SCMPCC exhibits the worst behavior because it plans for the worst-case
scenario (TP crosses intersection) for a long time, i.e., too conservatively.
The other methods all exhibit satisfactory behavior for this scene. In
conclusion, our method succeeds in both cases, while the others only perform
well in one of them.

Figure 6: Qualitative comparison of the performance of our approach and the
Baselines in the traffic intersection scene

###  VI-B Highway Merging

We conduct a Monte Carlo analysis, sampling 100 random highway merging
scenarios (s. Tab. I). Fig. 6 depicts an example of the evolution of a single
merging scene using our method. To simulate TP behavior, we utilize the
Intelligent Driver Model(IDM), with IDM parameters randomly sampled. For
method comparison, we assess the rate of merging successes, occurrences of
getting stuck on the merging lane (aborted), and collisions. Additionally, we
evaluate computation time (including time for scenario selection and adaptive
decision postponing) and, again, the costs. Results indicate, as expected,
that SCMPC is very safe but tends to be overly conservative. Our approach
achieves a balance between robustness and performance. The impact of scenario
selection is best observed in the comparison between BMPCCnoDP and BMPCCnoSS
(with 2, 3, and 4 scenarios). With scenario selection, BMPCCnoDP partly
outperforms BMPCC with 3 or 4 scenarios, with significantly reducing
computation time. Furthermore, the distinct effect of adaptive decision
postponing is evident in the comparison between our framework and BMPCCnoDP.

Figure 7:  AV behavior at three different timesteps. Initially, the AV does
not commit to merging in front or behind the first TP. As uncertainty reduces,
the branches converge.

##  VII Conclusion

We proposed a framework for motion planning and control that addresses multi-
modal uncertainty of traffic participants. By integrating a learning-based
multi-modal motion predictor into a Branch MPCC we improve AV safety without
excessive conservatism. However, this enhancement required designing
intermediary modules. Specifically, we developed a scenario selection method
leveraging topology and risk metrics as selection criteria of multi-modal
predictions. By reducing the number of considered predictions, we rendered
Branch MPC tractable for real-time applications. Our findings demonstrate that
few predictions are sufficient if the selected predictions cover this
uncertainty well enough. Furthermore, we introduced an adaptive decision-
postponing approach to estimate when sufficient future information allows to
commit to one of the predictions, i.e., when branching in the Branch MPCC
should begin. Our evaluation has shown that this also has a major impact on
the performance.

## References

  * [1] C. Hubmann, J. Schulz, M. Becker, D. Althoff, and C. Stiller, “Automated Driving in Uncertain Environments: Planning With Interaction and Uncertain Maneuver Prediction,” _IEEE Transactions on Intelligent Vehicles_ , vol. 3, no. 1, pp. 5–17, Mar. 2018.
  * [2] A. Somani, N. Ye, D. Hsu, and W. S. Lee, “DESPOT: Online POMDP Planning with Regularization,” in _Advances in Neural Information Processing Systems_ , vol. 26. Curran Associates, Inc., 2013.
  * [3] G. Schildbach, M. Soppert, and F. Borrelli, “A collision avoidance system at intersections using Robust Model Predictive Control,” in _2016 IEEE Intelligent Vehicles Symposium (IV)_ , Jun. 2016, pp. 233–238.
  * [4] S. Dixit, U. Montanaro, M. Dianati, D. Oxtoby, T. Mizutani, A. Mouzakitis, and S. Fallah, “Trajectory planning for autonomous high-speed overtaking in structured environments using robust mpc,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 21, no. 6, pp. 2310–2323, 2020.
  * [5] T. Benciolini, T. Brudigam, and M. Leibold, “Multistage Stochastic Model Predictive Control for Urban Automated Driving,” in _2021 IEEE International Intelligent Transportation Systems Conference (ITSC)_. Indianapolis, IN, USA: IEEE, Sep. 2021, pp. 417–423.
  * [6] T. Brudigam, M. Olbrich, M. Leibold, and D. Wollherr, “Combining Stochastic and Scenario Model Predictive Control to Handle Target Vehicle Uncertainty in an Autonomous Driving Highway Scenario,” in _2018 21st International Conference on Intelligent Transportation Systems (ITSC)_. Maui, HI: IEEE, Nov. 2018, pp. 1317–1324.
  * [7] R. Kensbock, M. Nezami, and G. Schildbach, “Scenario-Based Decision-Making, Planning and Control for Interaction-Aware Autonomous Driving on Highways,” in _2023 IEEE Intelligent Vehicles Symposium (IV)_ , Jun. 2023, pp. 1–6, iSSN: 2642-7214.
  * [8] O. de Groot, B. Brito, L. Ferranti, D. Gavrila, and J. Alonso-Mora, “Scenario-Based Trajectory Optimization in Uncertain Dynamic Environments,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 3, pp. 5389–5396, 2021.
  * [9] R. Oliveira, S. H. Nair, and B. Wahlberg, “Interaction and decision making-aware motion planning using branch model predictive control,” in _2023 IEEE Intelligent Vehicles Symposium (IV)_ , 2023.
  * [10] Y. Chen, U. Rosolia, W. Ubellacker, N. Csomay-Shanklin, and A. D. Ames, “Interactive multi-modal motion planning with branch model predictive control,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 5365–5372, 2022.
  * [11] J. P. Alsterda, M. Brown, and J. C. Gerdes, “Contingency Model Predictive Control for Automated Vehicles,” in _2019 American Control Conference (ACC)_. Philadelphia, PA, USA: IEEE, Jul. 2019, p. 717.
  * [12] V. Fors, B. Olofsson, and E. Frisk, “Resilient Branching MPC for Multi-Vehicle Traffic Scenarios Using Adversarial Disturbance Sequences,” _IEEE Transactions on Intelligent Vehicles_ , vol. 7, no. 4, Dec. 2022.
  * [13] C. H. Ulfsjoo and D. Axehill, “On Integrating POMDP and Scenario MPC for Planning under Uncertainty – with Applications to Highway Driving,” in _2022 IEEE Intelligent Vehicles Symposium (IV)_. Aachen, Germany: IEEE, Jun. 2022, pp. 1152–1160.
  * [14] Ö. Ş. Taş, P. H. Brusius, and C. Stiller, “Decision-theoretic MPC: Motion Planning with Weighted Maneuver Preferences Under Uncertainty,” Oct. 2023, arXiv:2310.17963 [cs, math] version: 1.
  * [15] S. Shi, L. Jiang, D. Dai, and B. Schiele, “Motion transformer with global intention localization and local movement refinement,” _Advances in Neural Information Processing Systems_ , vol. 35, pp. 6531–6543, 2022.
  * [16] N. Nayakanti, R. Al-Rfou, A. Zhou, K. Goel, K. S. Refaat, and B. Sapp, “Wayformer: Motion Forecasting via Simple & Efficient Attention Networks,” _arXiv preprint arXiv:2207.05844_ , 2022.
  * [17] Z. Zhou, J. Wang, Y. Li, and Y. Huang, “Query-centric trajectory prediction,” in _2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2023, pp. 17 863–17 873.
  * [18] B. Brito, B. Floor, L. Ferranti, and J. Alonso-Mora, “Model Predictive Contouring Control for Collision Avoidance in Unstructured Dynamic Environments,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 4, pp. 4459–4466, 2019.
  * [19] F. Siebenrock, M. Günther, and S. Hohmann, “LTV-MPC Based Trajectory Planning Considering Uncertain Object Prediction Through Adaptive Potential Fields,” in _2020 IEEE Conference on Control Technology and Applications (CCTA)_. IEEE, 2020, pp. 666–672.
  * [20] F. Altché, P. Polack, and A. de La Fortelle, “High-speed trajectory planning for autonomous vehicles using a simple dynamic model,” in _2017 IEEE 20th International Conference on Intelligent Transportation Systems (ITSC)_ , 2017, pp. 1–7.
  * [21] O. de Groot, L. Ferranti, D. Gavrila, and J. Alonso–Mora, “Globally guided trajectory planning in dynamic environments,” in _2023 IEEE International Conference on Robotics and Automation (ICRA)_ , 2023, pp. 10 118–10 124.
  * [22] A. Philipp and D. Goehring, “Analytic Collision Risk Calculation for Autonomous Vehicle Navigation,” in _2019 International Conference on Robotics and Automation (ICRA)_. Montreal, QC, Canada: IEEE, May 2019, pp. 1744–1750.
  * [23] A. Bhattacharyya, “On a measure of divergence between two multinomial populations,” _Sankhyā: The Indian Journal of Statistics (1933-1960)_ , vol. 7, no. 4, pp. 401–406, 1946.


11institutetext: Instituto de Astrofísica de Andalucía-CSIC, Glorieta de la
Astronomía s/n, 18008 Granada, Spain. 22institutetext: Max-Planck-Institut für
Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany. 33institutetext:
INAF Istituto di Radioastronomia, via Gobetti 101, 40129 Bologna, Italy.
44institutetext: Foundation for Research and Technology - Hellas, IESL &
Institute of Astrophysics, Voutes, 7110 Heraklion, Greece. 55institutetext:
Department of Physics, University of Crete, GR-70013 Heraklion, Greece.
66institutetext: Institute for Astrophysical Research, Boston University, 725
Commonwealth Avenue, Boston, MA 02215. 77institutetext: Astronomical
Institute, St. Petersburg State University, Universitetskij Pr. 28,
Petrodvorets, St. Petersburg 198504, Russia. 88institutetext: Korea Astronomy
and Space Science Institute, Daedeok-daero 776, Yuseong-gu, Daejeon 34055,
Republic of Korea. 99institutetext: Santa Cruz Institute for Particle Physics
and Department of Physics, University of California, Santa Cruz, CA 95064,
USA. 1010institutetext: Physics Department, McGill University, Montreal, QC
H3A 2T8, Canada. 1111institutetext: Institut de Radioastronomie
Millimétriquea, Avenida Divina Pastora 7, Local 20, E-18012 Granada, Spain.

# New jet feature in the parsec-scale jet of the blazar OJ 287 connected to
the 2017 teraelectronvolt flaring activity

R. Lico, Email<EMAIL_ADDRESS>C. Casadio 445522 S. G. Jorstad 6677 J. L.
Gómez 11 A. P. Marscher 66 E. Traianou 1122 J.-Y. Kim 8822 G.-Y. Zhao 11 A.
Fuentes 11 I. Cho 11 T. P. Krichbaum 22 O. Hervet 99 S. O’Brien 1010 B.
Boccardi 22 I. Myserlis 111122 I. Agudo 11 A. Alberdi 11 Z. R. Weaver 66 J. A.
Zensus 22

(Received …; accepted …)

###### Abstract

Context. In February 2017 the blazar OJ 287, one of the best super-massive
binary-black-hole-system candidates, was detected for the first time at very
high energies (VHEs; E$>100$ GeV) with the ground-based $\gamma$-ray
observatory VERITAS.

Aims. Very high energy $\gamma$ rays are thought to be produced in the near
vicinity of the central engine in active galactic nuclei. For this reason, and
with the main goal of providing useful information for the characterization of
the physical mechanisms connected with the observed teraelectronvolt flaring
event, we investigate the parsec-scale source properties by means of high-
resolution very long baseline interferometry observations.

Methods. We use 86 GHz Global Millimeter-VLBI Array (GMVA) observations from
2015 to 2017 and combine them with additional multiwavelength radio
observations at different frequencies from other monitoring programs. We
investigate the source structure by modeling the brightness distribution with
two-dimensional Gaussian components in the visibility plane.

Results. In the GMVA epoch following the source VHE activity, we find a new
jet feature (labeled K) at $\sim 0.2$ mas from the core region and located in
between two quasi-stationary components (labeled S1 and S2). Multiple periods
of enhanced activity are detected at different radio frequencies before and
during the VHE flaring state.

Conclusions. Based on the findings of this work, we identify as a possible
trigger for the VHE flaring emission during the early months of 2017 the
passage of a new jet feature through a recollimation shock (represented by the
model-fit component S1) in a region of the jet located at a de-projected
distance of $\sim 10$ pc from the radio core.

###### Key Words.:

Galaxies: active – BL Lacertae objects: OJ 287 – Galaxies: jets

##  1 Introduction

Super-massive binary black-hole (SMBBH) systems are commonly believed to form
during hierarchical galaxy formation as an outcome of merger events (e.g.,
Begelman et al. 1980). The blazar OJ 287 ($z=0.306$) is currently considered
one of the best SMBBH system candidates, supported by a number of lines of
observational evidence and theoretical arguments (Lehto & Valtonen 1996;
Valtonen et al. 2008). While the vast majority of blazar objects show erratic
variability across the electromagnetic spectrum, from the optical light curve
of OJ 287 two quasi-periodic variability patterns have been revealed, on
timescales of $\sim 12$ yr and $\sim 60$ yr, respectively (Valtonen et al.
2006). Within the context of a gravitationally bound system of two black
holes, the $\sim 12$ yr and $\sim 60$ yr optical outbursts are connected with
the orbital period and the periastron advance timescales, respectively.
Additional scenarios, not necessarily requiring the presence of a secondary
black hole, are also being investigated with the goal of explaining the
optical emission variability (e.g., Britzen et al. 2018; Liska et al. 2018).

Figure 1: Natural weighted total intensity images for the five 86 GHz GMVA
observing epochs from 2015 (left) to 2017 (right). The MJD epoch is indicated
in the top-left corner in each image. The restoring beam size, shown in the
bottom-right corner in each image, is 0.29 mas $\times$ 0.06 mas (at -12.8∘),
0.25 mas $\times$ 0.04 mas (at -10.1∘), 0.23 mas $\times$ 0.06 mas (at
-11.5∘), 0.21 mas $\times$ 0.05 mas (at -11.2∘), and 0.24 mas $\times$ 0.04
mas (at -6.5∘) from left to right, respectively. The overlaid lowest total
intensity contour is at $1\%$ of the map peak (0.54, 1.11, 2.30, 2.53, and
0.68 Jy/beam from left to right, respectively), with the contours that follow
a factor of two higher. The image rms noise from left to right is 0.54, 1.03,
2.71, 4.15, and 3.47 mJy/beam, respectively. The color scale represents the
total intensity. The overlaid red circles represent the model-fit components.

Since February 5, 2017, OJ 287 has aroused considerable interest for another
reason: the Very Energetic Radiation Imaging Telescope Array System (VERITAS)
collaboration reported the first ever very high-energy (VHE; E$>100$ GeV)
$\gamma$-ray emission detection of OJ 287 (Mukherjee & VERITAS Collaboration
2017). Different theoretical models related to the mass accretion, the
relativistic plasma outflows, and the particle pair-production processes have
been considered for explaining the production of such VHE emission in active
galactic nuclei (AGNs; e.g., Rieger & Levinson 2018). We note that blazars
dominate the population of the VHE $\gamma$-ray-detected AGNs; this is mainly
because of the strong Doppler boosting effects produced in their relativistic
jets, which are closely aligned to our line-of-sight. Currently, $93\%$ of the
AGNs listed in the online TeVCat111http://tevcat2.uchicago.edu/ catalog are
blazars (74 out of 80 objects). However, the OJ 287-like blazars, the so-
called low-frequency peaked blazars, with their synchrotron emission component
peaking at frequencies $<10^{14}$ Hz, represent less than $15\%$ of the whole
TeVCat blazar sample.

In the present work we investigate the source parsec-scale properties by means
of high-resolution very long baseline interferometry (VLBI) observations, with
the main goal of helping to characterize the physical mechanisms connected
with the reported teraelectronvolt (TeV) emission activity. We make use of
five 86 GHz Global Millimeter-VLBI Array (GMVA) observations, covering the
period 2015-2017, which allow us to access the innermost jet region at a
resolution of $\sim 50$ $\mu$as ($\sim 0.2$ pc in linear scale), complemented
by additional multifrequency radio observations from other monitoring
programs.

The paper is laid out as follows: observations and methods in Sect. 2; main
results in Sect. 3; and general discussion and concluding remarks in Sect. 4.
Throughout the paper we use a $\Lambda$ cold dark matter cosmology with
$H_{0}=67.4$ km s-1 Mpc-1, $\Omega_{m}=0.31$, and $\Omega_{\Lambda}=0.69$
(Planck Collaboration et al. 2020). At a redshift of $z=0.306$ (Nilsson et al.
2010), 1 mas corresponds to $\sim 4.7$ pc.

##  2 Observations and data analysis

In this work we analyze five 86 GHz GMVA observations of OJ 287 performed from
May 16, 2015, to September 29, 2017 – modified Julian days (MJDs) 57158,
57529, 57661, 57843, and 58025 – in total intensity emission. These
observations are part of a monitoring
program222http://www.bu.edu/blazars/vlbi3mm/ (PIs: A. P. Marscher and T. P.
Krichbaum) of a sample of $\gamma$-ray-bright blazars and radio galaxies
selected among the targets of the VLBA-BU-BLAZAR
program333https://www.bu.edu/blazars/VLBAproject.html.

Figure 2: Time evolution of the position and flux density for each model-fit
component. Left image: Time evolution of the distance from the core
(represented by the horizontal solid black line) of each model-fit component.
The light and dark blue circles represent components S1 and S2, respectively.
The green squares and red crosses represent components E and K, respectively.
The vertical dashed line represents the TeV detection (MJD 57785). Right
image: Light curves for each model-fit component. The core component is
represented by the empty black circles, and all the other symbols are as in
the left image. The nominal total intensity emission uncertainties were
calculated by considering a calibration error of $\sim 20\%$ of the flux
density and a statistical error equal to three times the map rms noise, but no
amplitude scaling factors are provided in this work (see Sect. 2 for more
details).

The GMVA stations that took part in the different observations presented in
this work are: the Robert C. Byrd Green Bank Telescope; eight stations of the
Very Long Baseline Array (VLBA): BR, FD, KP, LA, MK, NL, OV, and PT; the
European stations Effelsberg, Yebes, Onsala, Metsäovi, and Pico Veleta; and
the Korean VLBI Network array. For the a priori phase and amplitude
calibration, we used the software package Astronomical Image Processing System
(Greisen 2003), while for the final images we used the CLEAN and self-
calibration procedures in the DIFMAP software package (Shepherd 1997). More
details about the calibration process of millimeter-VLBI data reduction are
available in, for example, Jorstad et al. (2017), Casadio et al. (2019), and
Kim et al. (2019).

In each observing epoch we modeled the source brightness distribution by means
of two-dimensional Gaussian circular components, using the model-fitting
routine in the DIFMAP package. The cross-identification of the model-fit
components across the different epochs is based on their position, size, and
flux density. The uncertainties in the position of the model-fit components
were estimated by using the ratio of the size of each component to the signal-
to-noise ratio (as in Lico et al. 2012).

We note that, for several known reasons, significant flux losses could affect
the GMVA observations (e.g., Casadio et al. 2019). The amplitude scaling
factors used to correct the GMVA flux density values are usually on the order
of 2.0-2.5 and are extrapolated by using close-in-time lower-frequency
observations (e.g., Koyama et al. 2016; Kim et al. 2019). Since OJ 287 is an
extremely variable source, showing significant radio and optical emission
variations on timescales of a few days (e.g., Liu et al. 2017), it is
important to use observations that are as close as possible in time to the
GMVA observing epochs. Although several monitoring programs provide us with
both single-dish and VLBI observations at different frequencies (as shown and
discussed in Sect.3.2), it was not possible to get observations that were
sufficiently contemporaneous to obtain reliable amplitude scaling factors, and
we thus do not provide them here. The main effect of the lack of the scaling
factors is that we cannot quantify the GMVA flux density evolution across the
different epochs. We overcame this by using observations at multiple radio
frequencies and resolutions during the same observing period (Sect. 3.2).
Furthermore, the scaling factors affect neither the source structure nor the
relative flux density among the different model-fit components within the same
observing epoch.

##  3 Results

###  3.1 86 GHz GMVA images and model fits

In Fig. 1 we report the five 86 GHz GMVA total intensity images. The source
shows a jet structure in the northwest direction, extending for $\sim 0.4$ mas
($\sim 2.0$ pc in linear scale). More details about the image beam and total
intensity contours can be found in the figure caption. As is apparent from a
visual inspection of the images, the source jet orientation varies across the
different epochs. This is a well-known and well-investigated feature of OJ287
(e.g., Cohen 2017; Britzen et al. 2018; Gómez et al. 2022). A combined
multifrequency ridgeline and position angle time evolution analysis of the OJ
287 jet, which includes all the GMVA epochs reported in this Letter (except
for the last epoch, MJD 58025), was recently presented in Dey et al. (2021).

The red circles in each sub-image represent the model-fit components. A
summary of the model-fit parameters is reported in Table 1. The overall source
brightness distribution is well represented by four circular model-fit
components, labeled Core, S1, S2, and E. In agreement with previous works
(e.g., Agudo et al. 2012; Hodgson et al. 2017; Jorstad et al. 2017), we
identified the most upstream component, which is unresolved, bright, and
variable, as the core. The distance from the core of each model-fit component
(under the assumption that the core is stationary) is reported in the left
panel of Fig 2. Components S1 and S2 (light and dark blue circles in Fig 2,
respectively) are consistent with being stationary during the whole observing
period. We note that at lower frequencies, because of the lower resolution,
only a single quasi-stationary component has been detected in the images of
the OJ 287 jet since 2005, at an average distance of $0.2$ mas from the core,
as reported in several OJ 287 kinematics studies (e.g., Agudo et al. 2012;
Jorstad et al. 2017; Sasada et al. 2018; Weaver et al. 2021). The outermost
jet component, labeled E, represents the more extended emission.

A new model-fit component, labeled K (red cross in Fig 2), was detected in the
region between S1 and S2 on MJD 57843 (March 31, 2017), about two months after
the reported TeV detection. Although the identification of K is not obvious
and straightforward, we identified the different model-fit components on MJD
57843 according to their distance from the core and flux density across the
different epochs.

###  3.2 Light curves

The light curves of the five GMVA epochs for the different model-fit
components are reported in the right panel of Fig 2. We find that on MJD 57661
(September 30, 2016) the overall source emission was dominated by the total
intensity emission from component S1, while on MJD 57843 (March 31, 2017) it
was dominated by the emission from the new component, K. During the last
observing epoch, MJD 58025, component S2 became brighter than the core and S1.

To investigate the overall radio emission activity during the period
2015-2017, we produced the light curves at different radio frequencies from
different monitoring programs (Fig. 3). We used: 15 GHz single-dish
observations from the Owens Valley Radio Observatory (OVRO) blazar monitoring
program (Richards et al. 2011); 86 GHz single-dish observations from the
Polarimetric Monitoring of AGN at Millimetre Wavelengths (POLAMI)
program444http://polami.iaa.es/ (Agudo et al. 2018a, b; Thum et al. 2018); 15
GHz VLBA observations from the Monitoring of Jets in Active galactic nuclei
with VLBA Experiments (MOJAVE)
program555https://www.physics.purdue.edu/astro/MOJAVE/; and 43 GHz VLBA
observations from the VLBA-BU-BLAZAR
program666https://www.bu.edu/blazars/VLBAproject.html. As shown in Fig. 3, two
periods of enhanced radio emission activity can be identified before (gray
shaded area) and during (red shaded area) the VHE flaring activity, and an
increasing trend is also observed toward the end of the observing window
(green shaded area) at all the different radio frequencies. We note that we
have not quantified the 86 GHz GMVA total intensity emission variability
across the different epochs because we do not provide the scaling factors (see
Sect. 2 for more details), and this is why in Fig. 3 the total GMVA flux
density (dark blue circles) is overall underestimated.

Figure 3: Multiwavelength light curves: 15 GHz single-dish OVRO (black
circles), 86 GHz single-dish POLAMI program (red circles), 86 GHz GMVA (dark
blue circles), 15 GHz VLBA MOJAVE program (light blue circles), and 43 GHz
VLBA-BU-BLAZAR program (green circles) observations from 2015-2017. The
vertical dashed line represents the time of the TeV detection (MJD 57785). The
shaded gray, red, and green areas represent periods of enhanced emission
activity (see Sects. 3.2 and 4 for more details).

##  4 Discussion and conclusion

During the period February 1-4, 2017 (MJD 57785-57788), the blazar OJ287 was
observed with VERITAS for a total exposure time of 13 hours, resulting in the
first ever VHE $\gamma$-ray emission detection of the source, at $>5.7\sigma$
significance level. The source VHE detection was found to be temporally
coincident with increased activity at X-ray and $\gamma$-ray energies
(Mukherjee & VERITAS Collaboration 2017), as well as at centimeter and
millimeter wavelengths, as shown in Fig.3. This event triggered several
multiwavelength observing efforts and follow-ups with the primary purpose of
characterizing the underlying physical mechanisms (e.g., O’Brien 2017; Gupta
et al. 2019).

In the present work we have investigated the innermost regions of the OJ 287
jet, where the VHE $\gamma$ rays are thought to be produced, by means of high-
resolution millimeter-VLBI observations with the GMVA. One of the main
findings is that during the MJD 57158 (March 31, 2017) 86 GHz GMVA epoch we
detected a new model-fit component (K), in the region between the quasi-
stationary components S1 and S2, that dominates the source total intensity
emission. This is the closest-in-time GMVA epoch to the VHE event detection,
separated by about two months (58 days). During the same period, an enhanced
overall activity is detected at different radio frequencies (shaded red area
in Fig. 3). These two findings, a high radio emission state and the detection
of a new model-fit component, could be related to the VHE event and the
passage (ejection) of K through (from) the S1 quasi-stationary jet component.
We note that S1, as investigated and suggested in several works (e.g., Agudo
et al. 2011; Sasada et al. 2018; Weaver et al. 2021), is considered a
recollimation shock. The passage of a new component through a recollimation
shock associated with an enhanced activity at high energies is a quite common
event that has been observed in several AGNs (e.g., Marscher et al. 2010;
Jorstad & Marscher 2016; Casadio et al. 2019).

The identification of K on MJD 57843 (March 31, 2017) is driven and supported
by the following considerations. On MJD 57661 (September 30, 2016), when the
source underwent another episode of enhanced radio emission activity (shaded
gray area in Fig. 3), the overall GMVA total intensity emission was dominated
by component S1, found to be more than two times brighter than the core and S2
(right image in Fig.2). We interpret this finding as an indication that during
that epoch the new component, K, is blended with and passing through S1 and
that it emerges from the other side in the early months of 2017. Within this
scenario, numerical simulations suggest that fluctuations in the position of a
recollimation shock undergoing the crossing of a moving disturbance are
expected (e.g., Gómez et al. 1997). We note that as K exits the S1 region and
approaches S2, as predicted by simulations, both S1 and S2 show fluctuations
in their positions with respect to the core (as is apparent in the left image
of Fig.2). In more detail, from MJD 57661 (September 30, 2016) to MJD 57843
(March 31, 2017) component S2 moves away from the core from $0.19\pm 0.01$ mas
to $0.29\pm 0.01$ mas, while S1 moves toward the core region from $0.11\pm
0.01$ mas to $0.08\pm 0.01$ mas.

If we assume that component K was blended together with S1 (on average at
$\sim 0.1$ mas from the core) during the VHE outburst (MJD 57785) and then
observed at a distance of 0.16 mas during the GMVA epoch MJD 57843, we obtain
a propagation speed of $0.32$ mas/year, corresponding to $\sim$4.8 c. We note
that more than half of the OJ 287 knots found at 43 GHz, from the VLBA-BU-
BLAZAR monitoring program, show a proper motion between 0.3 and 0.5 mas/year
(e.g., Jorstad et al. 2017).

The main findings, to be taken into account when drawing a plausible physical
scenario, can be summarized as follows: (1) detection of a new jet feature
between S1 and S2 in the GMVA image during the closest GMVA epoch to the TeV
detection; (2) high radio emission state during the TeV flaring activity
period from both single-dish and VLBI observations; (3) source total intensity
emission dominated by S1 in September 2016 (prior to the TeV detection); (4)
fluctuations in the positions of S1 and S2 with respect to the core as the new
component, K, leaves the S1 region and approaches S2. Based on all these
results, we conclude that the February 2017 TeV flaring activity in OJ 287
could have been triggered by the passage of K through a recollimation shock
represented by the quasi-stationary model-fit component S1. We note that, for
OJ 287, both the core and S1 regions were found to be suitable sites for
particle acceleration associated with moving shocks and responsible for high
energy flares (e.g., Agudo et al. 2011; Hodgson et al. 2017; Sasada et al.
2018). In the specific case of the 2017 OJ 287 TeV detection, we conclude that
the region where the VHE emission originates coincides with S1, located at
$\sim 0.1$ mas from the core, corresponding to a de-projected distance of
$\sim 10.0$ pc (by assuming a viewing angle of $2.6^{\circ}$ as reported in
Jorstad et al. 2017).

Some additional remarks about the last GMVA observing epoch, MJD 58025
(September 29, 2017), further support the above proposed scenario. On one
hand, during this last GMVA epoch component K is not detected and component E
has a position and relative flux density consistent with the previous epochs.
On the other hand, component S2 has a higher flux density than the core and
S1, unlike during the previous epochs after 2015. These findings may indicate
that during this epoch K had reached and was passing though component S2,
which produced an increase in its flux density. According to the speed that we
estimated for K ($0.32$ mas/year), it is expected to have reached the region
where component S2 is found by MJD 58025. A detailed analysis of the OJ 287
jet kinematic properties will be presented in a dedicated and forthcoming
paper based on observations from the 43 GHz VLBA-BU-BLAZAR monitoring (Weaver
et al. 2021). We also note that after the overall flux density decrease
following the TeV flaring activity, the light curves at different radio
frequencies reported in Fig. 3 (shaded green area) show a new increasing trend
toward the end of 2017.

Table 1: Summary of the model-fit parameters.

Obs. date | MJD | Comp. ID |  $S$aa$a$Flux density ($S$) in Jy; |  $\sigma_{S}$bb$b$uncertainties on $S$; |  rcc$c$$r$ and $PA$ represent the distance and position angle (from north to east) of the component in polar coordinates with respect to the core; |  $\sigma_{r}$dd$d$uncertainties on $r$; |  PAcc$c$$r$ and $PA$ represent the distance and position angle (from north to east) of the component in polar coordinates with respect to the core; |  FWHMee$e$radius of the circular Gaussian components.
---|---|---|---|---|---|---|---|---
yy/mm/dd |  |  | (Jy) | (Jy) | (mas) | (mas) | (deg) | (mas)
2015/05/16 | 57158 | Core | 0.34 | 0.07 | … | … | ….. | 0.006
|  | S1 | 0.32 | 0.06 | 0.07 | 0.01 | -22.2 | 0.049
|  | S2 | 0.23 | 0.05 | 0.21 | 0.01 | -46.2 | 0.051
|  | E | 0.17 | 0.03 | 0.39 | 0.01 | -40.9 | 0.122
2016/05/21 | 57529 | Core | 1.08 | 0.22 | … | … | ….. | 0.017
|  | S1 | 0.53 | 0.11 | 0.10 | 0.01 | -32.2 | 0.033
|  | S2 | 0.06 | 0.01 | 0.15 | 0.01 | -68.4 | 0.004
2016/09/30 | 57661 | Core | 0.78 | 0.16 | … | … | ….. | 0.006
|  | S1 | 2.00 | 0.40 | 0.11 | 0.01 | -23.2 | 0.030
|  | S2 | 0.60 | 0.12 | 0.19 | 0.01 | -42.8 | 0.090
|  | E | 0.04 | 0.01 | 0.36 | 0.05 | -83.5 | 0.095
2017/03/31 | 57843 | Core | 0.83 | 0.17 | … | … | … | 0.016
|  | S1 | 0.57 | 0.12 | 0.08 | 0.01 | -41.8 | 0.005
|  | S2 | 0.25 | 0.05 | 0.29 | 0.01 | -47.0 | 0.038
|  | K | 1.39 | 0.28 | 0.16 | 0.01 | -40.6 | 0.038
|  | E | 0.06 | 0.02 | 0.38 | 0.05 | -50.2 | 0.096
2017/09/29 | 58025 | Core | 0.58 | 0.12 | … | … | … | 0.041
|  | S1 | 0.81 | 0.16 | 0.09 | 0.01 | -22.6 | 0.067
|  | S2 | 1.08 | 0.22 | 0.25 | 0.01 | -24.3 | 0.068
|  | E | 0.06 | 0.02 | 0.42 | 0.04 | -39.6 | 0.040
777

In conclusion, with the current work we (1) identify as a possible trigger of
the 2017 TeV enhanced activity the passage of a new jet feature through a
recollimation shock and (2) determine that the jet region where this occurs is
located at a de-projected distance of $\sim 10$ pc from the core. These
findings add a new important piece of information to a wider ongoing
multiwavelength analysis that will help to constrain emission models and
characterize details about the mechanism(s) producing the observed VHE
emission (VERITAS collaboration et al. in prep.). A detailed multiband optical
photometric and polarimetric analysis of OJ 287, covering the period from
September 2016 to December 2017, was presented in Gupta et al. (2019), and the
86 GHz polarization analysis with GMVA+ALMA observations from 2017 will be
presented in a forthcoming publication by Zhao et al. Furthermore, the results
of this Letter are meant to complement and be part of a series of other
ongoing works that feature the analysis of OJ 287 at a finer resolution by
means of 230 GHz Event Horizon Telescope observations (Gómez et al. in prep)
and 22 GHz space-VLBI RadioAstron observations (Lico et al. in prep).

###### Acknowledgements.

This research has made use of data obtained with the Global Millimeter VLBI
Array (GMVA), which consists of telescopes operated by the MPIfR, IRAM,
Onsala, Metsahovi, Yebes, the Korean VLBI Network, the Greenland Telescope,
the Green Bank Observatory and the Very Long Baseline Array (VLBA). The VLBA
and the GBT are a facility of the National Science Foundation operated under
cooperative agreement by Associated Universities, Inc. The data were
correlated at the correlator of the MPIfR in Bonn, Germany. RL, JLG, GYZ, AF,
TT, IC, IA and AA acknowledge financial support from the State Agency for
Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa”
award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709), from the
Spanish Ministerio de Economía y Competitividad, and Ministerio de Ciencia e
Innovación (grants AYA2016-80889-P, PID2019-108995GB-C21,
PID2019-107847RB-C44), the Consejería de Economía, Conocimiento, Empresas y
Universidad of the Junta de Andalucía (grant P18-FR-1769), the Consejo
Superior de Investigaciones Científicas (grant 2019AEP112). This study makes
use of 43 GHz VLBA data from the VLBA-BU Blazar Monitoring Program (BEAM-ME
and VLBA-BU-BLAZAR;http://www.bu.edu/blazars/VLBAproject.html), funded by NASA
through the Fermi Guest Investigator Program. The research at Boston
University was supported in part by NASA Fermi Guest Investigator program
grant 80NSSC20K1567. The VLBA is an instrument of the National Radio Astronomy
Observatory. The National Radio Astronomy Observatory is a facility of the
National Science Foundation operated by Associated Universities, Inc. The
POLAMI observations were carried out at the IRAM 30m Telescope. IRAM is
supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This research
has made use of data from the MOJAVE database that is maintained by the MOJAVE
team (Lister et al. 2018). This research has made use of data from the OVRO
40-m monitoring program (Richards et al. 2011), supported by private funding
from the California Insitute of Technology and the Max Planck Institute for
Radio Astronomy, and by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and
NSF grants AST-0808050 and AST- 1109911.

## References

  * Agudo et al. (2011) Agudo, I., Jorstad, S. G., Marscher, A. P., et al. 2011, ApJ, 726, L13.
  * Agudo et al. (2012) Agudo, I., Marscher, A. P., Jorstad, S. G., et al. 2012, ApJ, 747, 63.
  * Agudo et al. (2018a) Agudo, I., Thum, C., Ramakrishnan, V., et al. 2018, MNRAS, 473, 1850.
  * Agudo et al. (2018b) Agudo, I., Thum, C., Molina, S. N., et al. 2018, MNRAS, 474, 1427.
  * Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307.
  * Britzen et al. (2018) Britzen, S., Fendt, C., Witzel, G., et al. 2018, MNRAS, 478, 3199 .
  * Casadio et al. (2019) Casadio, C., Marscher, A. P., Jorstad, S. G., et al. 2019, A&A, 622, A158.
  * Cohen (2017) Cohen, M. 2017, Galaxies, 5, 12. doi:10.3390/galaxies5010012
  * Dey et al. (2021) Dey, L., Valtonen, M. J., Gopakumar, A., et al. 2021, MNRAS, 503, 4400.
  * Gómez et al. (1997) Gómez, J. L., Martí, J. M., Marscher, A. P., et al. 1997, ApJ, 482, L33.
  * Gómez et al. (2022) Gómez, J. L., Traianou, E., Krichbaum, T. P., et al. 2022, ApJ, 924, 122
  * Greisen (2003) Greisen, E. W. 2003, ASSL, 285, 109
  * Gupta et al. (2019) Gupta, A. C., Gaur, H., Wiita, P. J., et al. 2019, ApJ, 157, 95.
  * Hodgson et al. (2017) Hodgson, J. A., Krichbaum, T. P., Marscher, A. P., et al. 2017, A&A, 597, A80.
  * Jorstad & Marscher (2016) Jorstad, S. & Marscher, A. 2016, Galaxies, 4, 47.
  * Jorstad et al. (2017) Jorstad, S. G., Marscher, A. P., Morozova, D. A., et al. 2017, ApJ, 846, 98.
  * Kim et al. (2019) Kim, J.-Y., Krichbaum, T. P., Marscher, A. P., et al. 2019, A&A, 622, A196.
  * Koyama et al. (2016) Koyama, S., Kino, M., Giroletti, M., et al. 2016, A&A, 586, A113.
  * Lehto & Valtonen (1996) Lehto, H. J., & Valtonen, M. J. 1996, ApJ, 460, 207.
  * Lico et al. (2012) Lico, R., Giroletti, M., Orienti, M., et al. 2012, A&A, 545, A117.
  * Liska et al. (2018) Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2018, MNRAS, 474, L81.
  * Lister et al. (2018) Lister, M. L., Aller, M. F., Aller, H. D., et al. 2018, ApJ, 234, 12.
  * Liu et al. (2017) Liu, X., Yang, P. P., Liu, J., et al. 2017, MNRAS, 469, 2457.
  * Marscher et al. (2010) Marscher, A. P., Jorstad, S. G., Larionov, V. M., et al. 2010, ApJ, 710, L126.
  * Mukherjee & VERITAS Collaboration (2017) Mukherjee, R. & VERITAS Collaboration 2017, The Astronomer’s Telegram, 10051 (2017).
  * Nilsson et al. (2010) Nilsson, K., Takalo, L. O., Lehto, H. J., et al. 2010, A&A, 516, A60
  * O’Brien (2017) O’Brien, S. 2017, arXiv:1708.02160
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6.
  * Richards et al. (2011) Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJ, 194, 29.
  * Rieger & Levinson (2018) Rieger, F. & Levinson, A. 2018, Galaxies, 6, 116.
  * Sasada et al. (2018) Sasada, M., Jorstad, S., Marscher, A. P., et al. 2018, ApJ, 864, 67.
  * Shepherd (1997) Shepherd, M. C. 1997, ASPC, 125, 77.
  * Thum et al. (2018) Thum, C., Agudo, I., Molina, S. N., et al. 2018, MNRAS, 473, 2506.
  * Valtonen et al. (2006) Valtonen, M. J., Lehto, H. J., Sillanpää, A., et al. 2006, ApJ, 646, 36.
  * Valtonen et al. (2008) Valtonen, M. J., Lehto, H. J., Nilsson, K., et al. 2008, Nature, 452, 851.
  * Weaver et al. (2021) Weaver, Z. R., Jorstad, S. G., et al. 2021, ApJ, submitted.