# Light Travel Time Effects in Blazar Flares

Justin D. Finke1,∗

###### Abstract

I present a model for light travel time effects for emission from a plasma
blob in a blazar jet. This calculation could be incorporated into more complex
models with particle acceleration and radiation mechanisms, but as presented
here it is a agnostic as to these mechanisms. This model includes light travel
time effects for an expanding or contracting blob. As an example, this model
is applied to a flare observed by VERITAS and MAGIC from Mrk 421 in 2013; and
to a flare observed by the Fermi Large Area Telescope from 3C 454.3 in 2010.

##  1 Keywords:

BL Lacertae object: general, BL Lacertae objects: individual (Mrk 421),
quasars: general, quasars: individual (3C 454.3), gamma-rays: general, gamma-
rays: galaxies, time-domain astronomy

##  2 Introduction

Blazars are active galactic nuclei (AGN) with relativistic jets moving close
to our line of sight. They are routinely detected at all wavelengths, from
radio to $\gamma$-rays, and often extremely variable at all wavelengths as
well. Consequently, a complete understanding of blazars requires a deep
understanding of the time domain. The Fermi Large Area Telescope (LAT) is
extremely useful for this, since it monitors the entire sky in $\approx 30$
MeV to 100 GeV $\gamma$-rays every 3 hours. Observations by the LAT can be
supplemented with observations by many other observatories, for instance by
Swift in space in the optical and X-rays; and on the ground by MAGIC,
H.E.S.S., and VERITAS at very-high energy (VHE) $\gamma$-rays, Very Large
Baseline Array, TANAMI, Effelsberg, Owens Valley Radio Observatory in the
radio, and SMARTS and the GASP-WEBT consortium in the optical (among others).
The loss of Fermi will be a major blow to the study of blazars, especially
considering there is no planned observatory that can replace its capabilities.

Blazar variability is often analyzed in terms of flares, explained by particle
acceleration and radiative cooling in a homogeneous one-zone nonthermal plasma
“blob”. Extensive time-dependent modeling analysis on flares has been
performed (e.g., Boettcher et al., 1997; Kirk et al., 1998; Chiaberge and
Ghisellini, 1999; Dermer and Schlickeiser, 2002; Joshi and Böttcher, 2011;
Zacharias and Schlickeiser, 2013; Dotson et al., 2015; Zacharias, 2023).
Particle acceleration could be from diffusive shock or magnetic reconnection.
Leptonic radiation processes include synchrotron, typically at low energies,
and Compton scattering of a number of radiation fields at high energies.
Hadronic processes may also be included (e.g., Petropoulou and Mastichiadis,
2012; Diltz et al., 2015). Light travel time effects are sometimes, though not
always, ignored, implicitly assuming that the acceleration and/or radiative
timescales are much larger than the light crossing timescale.

Observationally, flares are often characterized in terms of an exponential
rise and decay (e.g., Nalewajko, 2013; Meyer et al., 2019; Roy et al., 2019;
Bhatta et al., 2023). Among other properties, the flares are often analyzed in
terms of symmetry, with symmetric flares having the same rising and decaying
timescales. Symmetric flares are often taken as an indication that the flare
variability is dominated by light travel time rather than particle
acceleration or cooling. Most blazar flares, but not all, are consistent with
being symmetric within the uncertainties.

Here I explore a simple model for emission from a plasma blob in a blazar jet
that takes into account light travel time effects. The simple model is
agnostic as to the particle acceleration and radiation mechanisms, assuming
only that the acceleration and radiative timescales are much less than the
light crossing timescale. In Section 3 I explore a simple model from a non-
expanding blob, which will result in a symmetric light curve. This should be a
good approximation for the majority of blazar flares that show no evidence for
asymmetry. Next in Section 4 I explore a model for a blob that can have a size
that changes with time. This could reproduce asymmetric flares, while still
being agnostic as to the exact acceleration and radiation mechanisms. Finally
I conclude with a discussion in Section 5.

##  3 Flare From a Constant Size Blob

###  3.1 Formalism

Consider a homogeneous blob. If the emission of the entire blob is changing
simultaneously, the observer will see the portion of the blob closer to the
observer before the portion of the blob that is farther away. This effect was
explored by Chiaberge and Ghisellini (1999) and Joshi and Böttcher (2011), and
described by Zacharias and Schlickeiser (2013) for a spherical blob with
constant size (see also Finke and Becker, 2014, 2015). A blob in a blazar jet
will be moving at a high relativistic speed, but I use all quantities of the
blob in the observer’s frame, so no relativistic transformations are
necessary.

I take the emitting blob to be centered on $x=R$ with two-dimensional cross-
sectional surface described by $(R-x)^{g}+y^{g}=R^{g}$, where $(x,y)$ are the
standard Cartesian coordinates. The observer is in the $-x$ direction. See
Figure 1 for an illustration of the geometry. The three-dimensional surface is
created by rotating this surface around the $x$-axis. The blob will then have
a volume given by

$\displaystyle V$
$\displaystyle=2\pi\int_{0}^{R}dxR^{2}\left[1-\left(\frac{x}{R}\right)\right]^{2/g}$
$\displaystyle=2\pi R^{3}\int^{1}_{0}d\tilde{x}\ (1-\tilde{x}^{g})^{2/g}\ .$
(1)

Here $g=2$ corresponds to the case of a sphere. For a constant size blob with
this geometry, and with intrinsic time series light curve $F(t)$ starting at
$t=0$, the observed flux as a function of the observer’s time $t_{\rm obs}$ is

$\displaystyle F_{\rm obs}(t_{\rm obs})=\frac{\pi c}{V}\int^{\min(2R/c,t_{\rm
obs})}_{0}dtF(t_{\rm obs}-t)\left[R^{g}-|R-tc|^{g}\right]^{2/g}\ .$ (2)

For the spherical case, $g=2$, Equation (2) gives the result from Zacharias
and Schlickeiser (2013) and Finke and Becker (2015),

$\displaystyle F_{\rm obs}(t_{\rm obs})=\frac{3c}{R}\int^{\min(2R/c,t_{\rm
obs})}_{0}dtF(t_{\rm
obs}-t)\left[\frac{tc}{2R}-\left(\frac{tc}{2R}\right)^{2}\right]\ .$ (3)

For the shortest possible flare,

$\displaystyle F(t)=F_{0}\delta(t-t_{0})\ ,$ (4)

where $F_{0}$ is the fluence of the flare and $\delta(x)$ is the usual Dirac
delta function, one gets from Equation (2)

$\displaystyle F(t_{\rm obs})$ $\displaystyle=F_{0}\frac{\pi
c}{V}\left[R^{g}-|R-(t_{\rm obs}-t_{0})c|^{g}\right]^{2/g}$
$\displaystyle\times H\left(t_{\rm obs}-t_{0};0,\frac{2R}{c}\right)\ .$ (5)

Here I use the step function

$H(x;a,b)\equiv\left\\{\begin{array}[]{ll}1&a<x<b\\\
0&\mathrm{otherwise}\end{array}\right.\ .$ (6)

Figure 1:  An illustration of the geometry of the emitting blob (red) in the
jet with the observer labeled.

The solution using Equation (3.1) has four free parameters: $F_{0}$, $R$, $g$,
and $t_{0}$. This model will create a symmetric flare, with the rising and
falling timescales being equal. If light travel time is indeed the only source
of variability in a flare, the shape of the light curve should be the same as
all energies (wavelengths, frequencies), although the normalization would be
different. So if one had a light curve flux at energy $E_{1}$, $F(E_{1},t)$,
the flux at $E_{2}$ would be $F(E_{2},t)=r_{1}F(E_{1},t)$, where $r_{1}$ is
another free parameter that is constant with time;
$F(E_{3},t)=r_{2}F(E_{1},t)$ where $r_{2}$ is another free parameter; and so
on. Each additional light curve adds another free parameter.

In Figure 2 I plot this model for different values of the parameter $g$. This
parameter can control the general shape of the flare. For lower values of $g$,
it becomes more peaked.

Figure 2:  Example of the non-changing blob size model. The parameters are
$F_{0}=10^{-6}\ \mathrm{erg}\ \mathrm{cm}^{-2}$, $R=10^{15}\ \mathrm{cm}$,
$t_{0}=0$, and various values of $g$ as shown in the Figure.

.

This model can be a more physically-motivated way of computing the size scale
from a flare, although it should be clear that the size scale found here ($R$)
is in the observer frame, not the frame co-moving with the jet, as is often
desired.

As described in Section 2, a large number of blazar flares do indeed appear to
be symmetric, so this simple model should be able to provide a good fit to a
large number of blazar flare light curves. I next apply this model to two
bright, well-observed $\gamma$-ray flares taken from the literature.

###  3.2 2013 April flare from Mrk 421

A highly detailed multiwavelength campaign on the high synchrotron peaked BL
Lac object Mrk 421 took place in 2013 April, as reported by Acciari et al.
(2020). During this time the source was very bright, and exhibited several
extremely bright flares. Here I apply the simple light travel time model for a
non-expanding blob to the flare on 2013 April 15 (MJD 56397) as observed in
the $\gamma$-rays by MAGIC and VERITAS. This bright flare was also captured in
similar detail in the X-rays by NuSTAR, but I do not explore that here.

I do a Markov chain Monte Carlo Foreman-Mackey et al. (2013) with the model
described above in Section 3.1 to the 0.2-0.4 TeV, 0.4-0.8 TeV, and $>0.8$ TeV
light curve for the 2013 April 15 flare from Mrk 421. The data are taken from
Acciari et al. (2020). The result can be seen in the left side of Figure 3 and
model parameters can be found in Table 1. The model appears to be a reasonably
good fit to the data, although it misses some of the points at early and late
times. The light curves at the different energy ranges do look very similar,
indicating that the light travel time may dominate the variability for this
flare.

The fit deviates substantially from a spherical blob, with the fit resulting
in $g\approx 0.60$, versus $g=2$ for a spherical blob. There is no real reason
to think that the shape of the blob is spherical, this is often used just for
simplicity. The resulting radius of the blob, $R\approx 3\times 10^{15}\
\mathrm{cm}$, is consistent with previous blazar modeling.

Figure 3:  MAGIC and VERITAS light curves for Mrk 421 in 2013 April. The green
symbols indicate the 0.2-0.4 TeV band; the red indicates the 0.4-0.8 TeV band;
and the blue indicates the $>0.8$ TeV band. The shaded regions show the 68%
confidence intervals from the model MCMC result. Left: constant size model.
Right: changing size model.

###  3.3 2010 November flare from 3C 454.3

In 2010 November 3C 454.3 exhibited a bright $\gamma$-ray outburst detected by
the Fermi-LAT that was the brightest $\gamma$-ray flare from a blazar up until
that point (Abdo et al., 2011). I also apply the non-expanding blob model to
the 0.1-1.0 and $>1.0$ GeV light curve of the brightest flare from this
outburst. The result can be seen in the left side of Figure 4 and the
parameters can be seen in Table 2. As can be seen, the model does not provide
a particularly good fit to the data. Further, the shape of the light curve at
the two different energy ranges is quite different, indicating that the simple
light travel time model may not be a good approximation for this burst. It is
interesting to note that the fit results in $g\approx 2$, consistent with a
spherical geometry for the emitting region.

Figure 4:  Fermi-LAT light curves for 3C 454.3 in 2010 November and model
results. The green symbols indicate the 0.1-1.0 GeV band; the red indicates
the $>1.0$ GeV band. The shaded regions show the 68% confidence intervals from
the model MCMC result. Left: constant size model. Right: changing size model.

##  4 Flare From a Blob That is Changing Size

###  4.1 Formalism

Here I generalize the light travel time effect for an expanding or contracting
axisymmetric blob, which, to the best of my knowledge, has not been explored
before. In this case, $R\rightarrow R(t)$. The solution for a constant size
blob $\delta$ function, Equation (3.1) above, can be used as a Green’s
function for the expanding blob case (letting $F_{0}\rightarrow 1$):

$\displaystyle G(t_{\rm obs},t_{0})$ $\displaystyle=\frac{\pi
c}{V(t_{0})}\left[R(t_{0})^{g}-|R(t_{0})-(t_{\rm
obs}-t_{0})c|^{g}\right]^{2/g}$ $\displaystyle\times H\left(t_{\rm
obs}-t_{0};0,\frac{2R(t_{0})}{c}\right)\ .$ (7)

The solution for the expanding/contracting blob will then be

$\displaystyle F_{\rm obs}(t_{\rm obs})=\int^{t_{\rm
obs}}_{0}dt_{0}F(t_{0})G(t_{\rm obs},t_{0})\ ,$ (8)

where again $F(t_{0})$ is the intrinsic light curve of the blob. Putting the
Green’s function, Equation (4.1) in Equation (8), one gets

$\displaystyle F_{\rm obs}(t_{\rm obs})=\pi c\int^{t_{\rm
obs}}_{\max(0,t_{\min})}dt_{0}\frac{F(t_{0})}{V(t_{0})}\left[R(t_{0})^{g}-|R(t_{0})-(t_{\rm
obs}-t_{0})c|^{g}\right]^{2/g}$ (9)

where $t_{\min}$ is defined by $t_{\rm obs}<t_{\min}+2R(t_{\min})/c$, based on
the step function in Equation (3.1). It can be determined once $R(t)$ is
specified. Equation (9) can be re-written with the substitution $t=t_{\rm
obs}-t_{0}$, leading to

$\displaystyle F_{\rm obs}(t_{\rm obs})=\pi c\int^{\min(t_{\max},t_{\rm
obs})}_{0}dt\frac{F(t_{\rm obs}-t)}{V(t_{\rm obs}-t)}\left[R(t_{\rm
obs}-t)^{g}-|R(t_{\rm obs}-t)-tc|^{g}\right]^{2/g}\ .$ (10)

Here $t_{\rm max}$ is defined by $t_{\rm max}<R(t_{\max})/c$, and must be
determined once $R(t)$ is specified. For a blob that is neither expanding nor
contracting, i.e., $R(t)$ is constant, Equation (10) reduces to Equation (2).

I take the blob length scale to expand as

$\displaystyle R(t_{0})=R_{0}\left(1+\frac{t_{0}-t_{\min}}{T}\right)^{a}\ ,$
(11)

which has free parameters $R_{0}$, $T$, and $a$. For this parameterization for
$R(t_{0})$, in Equation (9) $t_{\min}=t_{\rm obs}-2R_{0}/c$. In Equation (10),
$t_{\max}$ does not have a closed-form solution and so must be solved for
numerically.

For a $\delta$-function intrinsic flux, as used in the constant size case, the
expanding blob would simplify to the constant size case. This is because only
the flux at one infinitesimally small instant will be seen by the observer, so
any change in the blob’s size will not make a difference. So I take the
intrinsic light curve to be a Gaussian,

$\displaystyle
F(t_{0})=\frac{F_{0}}{\sqrt{2\pi}\sigma_{t}}\exp\left(\frac{-(t_{0}-t_{\rm
min})^{2}}{2\sigma_{t}^{2}}\right)\ .$ (12)

This model then has free parameters $F_{0}$, $R_{0}$, $g$, $t_{\min}$,
$\sigma_{t}$, $a$, $T$, along with any ratios for light curves at other
wavelengths ($r_{2}$, $r_{3}$). For $\sigma_{t}\rightarrow 0$, it is well
known that a Gaussian reduces to a Dirac $\delta$ function; so Equation (12)
reduces to Equation (4). This means constant size model is nested with the
changing size model, with the changing size model having 3 additional free
parameters.

An example of light curves with this model are shown in Figure 5 for different
values of the parameter $a$. Changes in $a$ can account for flare asymmetry,
allowing this model to explain a wide variety of flares with an emitting blob
that changes in size.

Figure 5:  Example of the changing blob size model. The parameters are
$F_{0}=10^{-6}\ \mathrm{erg}\ \mathrm{cm}^{-2}$, $\sigma_{t}=10^{3}\
\mathrm{s}$, $R_{0}=10^{15}\ \mathrm{cm}$, $T=10^{4}\ \mathrm{s}$, $g=2$,
$t_{\min}=0$, and various values of $a$ as shown in the Figure.

.

One prediction of this model is that the light curves for a source at
different wavelengths, produced by the same emitting region, should have the
same overall shape. This means if one observes simultaneous light curves at
different wavelengths, and their overall spaces are not consistent with one
another, this model alone is not sufficient. The light curves may differ to to
different radiative properties or mechanisms at the different wavelengths, or
contamination by other emitting regions.

###  4.2 2013 April flare from Mrk 421

I apply the changing size model described above to the 2013 April flare from
Mrk 421. The result can be seen in the right side of Figure 3 and Table 1. The
model does appear to provide a bit better fit. Since $a>0$, this indicates the
blob is expanding rather than contracting. A likelihood ratio test indicates
that the constant size blob model is rejected in favor of the changing blob
model at $4.8\sigma$ significance. Since the light curves at different
energies are so similar, this is a reasonably strong indication that the
variability is dominated by the light travel time for this flare.

###  4.3 2010 November flare from 3C 454.3

A also apply the changing size model to the 2010 November flare from 3C 454.3.
The model can be seen overplotted with the data in Figure 4, and the resulting
parameters can be seen in Table 2. The model does not appear to provide a
significantly better fit, and indeed looks very similar to the non-expanding
blob case. This is confirmed with the likelihood ratio test, which indicates
that the changing blob model is preferred at $10^{-7}\sigma$ over the non-
expanding model, which any reasonable researcher would interpret as the more
complicated model is not significantly preferred over the simpler one. The
reason for this seems to be that, unlike the flare in Mrk 421 explored above,
the light curve shapes for the different Fermi-LAT energy bins are quite
different (indeed, as noted by Abdo et al., 2011). This is an indication the
flare shape is not dominated by light travel time effects, and that a more
complicated model is needed to explain it.

##  5 Discussion

I have described a simple model for light travel time effects in blazar
flares, where the emitting region could be changing in size (either expanding
or contracting). The change in size of the region can lead to an asymmetry,
i.e., different rising and decay timescales. The constant size model uses an
instantaneous turning on and off of emission, i.e., a Dirac $\delta$ function
for the intrinsic emission. The expanding blob model assumes a narrow Gaussian
for the intrinsic emission. Although unrealistic, these models should be good
approximations for flares where the particle acceleration and energy loss
timescales are much less than the light travel timescale for the “blob”. For
many flares, including asymmetric flares, light travel time effects alone
should be able to account for the flare. Observing at multiple wavelengths is
a good way to confirm or rule out this model. If the flares have the same
shape at different wavelengths, can be explained by this model, this is a good
indication that the variability is dominated by light travel time. However, if
the light curve shapes during the flare are different at different
wavelengths, this could be an indication that other processes are important.
For instance, it could be an indication that the particle acceleration and
energy loss timescales are not much smaller than the light travel timescale,
and cannot be neglected. Although here one must be careful. For instance, in
comparing optical emission and $\gamma$-ray emission in FSRQs, one must be
careful to take into account the emission from the “blue bump”, i.e., from the
accretion disk.

This story is complicated by recent observations by the Imaging X-ray
Polarimetry Explorer (IXPE). Simultaneous observations of Mrk 421 (Di Gesu et
al., 2023) and PG 1553+113 (Middei et al., 2023) reveal optical polarizations
that are significantly different than X-ray polarizations measured by IXPE.
Mrk 421 and PG 1553+113 are both high synchrotron-peaked blazars, so their
optical and X-ray emission are both thought to be from synchrotron. One
possible explanation is that the electrons producing the optical synchrotron
can travel a larger distance from the acceleration site than the electron
producing the X-ray synchrotron, since the lower energy, optical-emitting
electrons will have a longer energy loss timescale than the higher energy,
X-ray emitting electrons (Zhang et al., 2024). In this case, one would indeed
expect the optical emission region to be larger than the X-ray emission
region, and thus the optical would have a larger variability timescale, and
likely a time delay, relative to the X-ray emission. In any case, light travel
time effects, such as those described here, would need to be taken into
account in accurate time-dependent modeling calculations.

As two examples, I applied this model to the 2013 April flare from the BL Lac
Mrk 421, as observed by MAGIC and VERITAS; and to the 2010 November flare from
the FSRQ 3C 454.3, as observed by the Fermi-LAT. In the case of the Mrk 421,
the expanding blob model provides a good fit to the data, and is preferred
over the constant size model with a significance of 4.8$\sigma$. This flare
was extremely bright, and the data are quite good. A more complicated would
probably not provide a better description of the data. For the flare from 3C
454.3, neither model provides a particularly good reproduction of the data,
and the changing size model is not significantly preferred over the constant
size model.

The models here could be usefully applied to studies of light curves making up
multiple flares (e.g., Meyer et al., 2019; Roy et al., 2019; Bhatta et al.,
2023). Since most flares are symmetric, the constant size blob model should be
sufficient to explain most flares. The parameter $g$ can have a substantial
impact on the shape of the flare, and varying it could probably describe most
flares (Figure 2). The fits can be used to at least put an upper limit on the
size scale of the flare (the parameter $R$). Care should be taken, since this
size scale is in the observer’s frame. In the frame co-moving with the blob,
the size scale could be larger by a factor of the Doppler factor. Also, the
data from most flares are not as good as the data for the flares explored
here. So a lack of improvement in the fit for most flares could be due to the
large uncertainty in the data, and may not be due to the flare being dominated
by light travel time effects.

Although these models do not contain any information about particle
acceleration or radiation mechanisms, incorporating light travel time, as
discussed here, can be incorporated into such models. To some extend light
travel time effects have been incorporated into more complicated models
already (e.g., Chiaberge and Ghisellini, 1999; Joshi and Böttcher, 2011;
Zacharias and Schlickeiser, 2013). I will incorporate the formalism described
here into more complicated models in future work.

Particle-in-cell simulations indicate that magnetic reconnection in blazar
jets could create “plasmoids”, i.e., magnetized, nonthermal plasma of various
sizes with acceleated, radiating particles (e.g., Petropoulou et al., 2018).
Christie et al. (2019) performed time-dependent modeling of emission from
plasmoids created by magnetic reconnection in a blazar jet. They included
light travel time effects between the different plasmoids, but did not include
the light travel time effects within the plasmoid. The light travel time
effect described here could be integrated into their model to take into
account intra-plasmoid light travel time, including for a plasmoid that is
changing in size. If the plasmoids are not changing in size rapidly during a
flare, the flares could be more symmetric than previously assumed, affecting
$\gamma$-ray observations and the inferred flare timescales (Meyer et al.,
2021).

Table 1: Model parameters for fits to 2013 April flare from Mrk 421.  Parameter | Constant Size Blob | Changing Size Blob
---|---|---
$t_{0}$ or $t_{\min}$ [MJD] | $56395.0^{+0.4}_{-0.3}$ | $56395.7\pm 0.3$
$\log_{10}F_{0}$ [erg cm-2] | $-4.59^{+0.04}_{-0.08}$ | $-4.42^{+0.29}_{-0.17}$
$\sigma_{t}$ [days] | - | $1.77^{+0.31}_{-0.26}$
$\log_{10}R_{0}$ [cm] | $15.46^{+0.10}_{-0.19}$ | $15.55\pm 0.09$
$T$ [days] | - | $3.09\pm 0.45$
$a$ | - | $2.7^{+2.1}_{-1.8}$
$g$ | $0.60^{+0.08}_{-0.03}$ | $0.49\pm 0.04$
$r_{1}$ | $0.51\pm 0.01$ | $0.51\pm 0.02$
$r_{2}$ | $0.36\pm 0.01$ | $0.36\pm 0.01$
$\chi^{2}$/dof | 56/24 | 24/21
Table 2: Model parameters for fits to 2010 November flare from 3C 454.3.  Parameter | Constant Size Blob | Changing Size Blob
---|---|---
$t_{0}$ or $t_{\min}$ [MJD] | $55518.17^{+0.03}_{-0.01}$ | $55518.17^{+0.03}_{-0.01}$
$\log_{10}F_{0}$ [erg cm-2] | $1.13^{+0.01}_{-0.01}$ | $1.13^{+0.01}_{-0.04}$
$\log_{10}\sigma_{t}$ [s] | - | $1.97^{+0.90}_{-0.68}$
$\log_{10}R_{0}$ [cm] | $15.68^{+0.01}_{-0.01}$ | $15.68^{+0.01}_{-0.01}$
$T$ [days] | - | $7.8^{+2.5}_{-2.7}$
$a$ | - | $-0.66^{+2.78}_{-3.51}$
$g$ | $2.07\pm 0.11$ | $2.08^{+0.11}_{-0.10}$
$r_{1}$ | $0.057\pm 0.02$ | $0.057\pm 0.03$
$\chi^{2}$/dof | 87/19 | 86/16

## Conflict of Interest Statement

The author declares that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

J.D.F. performed all of the work presented in this manuscript.

## Funding

I am supported by NASA through the Fermi GI program and contract S-15633Y; and
by the Office of Naval Research. This work made use of a grant of computing
time from the DoD’s High Performance Computing Modernization Program.

## Acknowledgements

I would like to thank the anonymous referees for helpful comments, and R.
Sambruna and F. M. Civano for the invitation to contribute to this collection.

## References

  * Abdo et al. (2011) Abdo, A. A., Ackermann, M., Ajello, M., Allafort, A., Baldini, L., Ballet, J., et al. (2011).  Fermi Gamma-ray Space Telescope Observations of the Gamma-ray Outburst from 3C454.3 in November 2010.  _ApJ_ 733, L26.  10.1088/2041-8205/733/2/L26
  * Acciari et al. (2020) Acciari, V. A., Ansoldi, S., Antonelli, L. A., Arbet Engels, A., Baack, D., Babić, A., et al. (2020).  Unraveling the Complex Behavior of Mrk 421 with Simultaneous X-Ray and VHE Observations during an Extreme Flaring Activity in 2013 April.  _ApJS_ 248, 29.  10.3847/1538-4365/ab89b5
  * Bhatta et al. (2023) Bhatta, G., Zola, S., Drozdz, M., Reichart, D., Haislip, J., Kouprianov, V., et al. (2023).  Catching profound optical flares in blazars.  _MNRAS_ 520, 2633–2643.  10.1093/mnras/stad280
  * Boettcher et al. (1997) Boettcher, M., Mause, H., and Schlickeiser, R. (1997).  $\gamma$-ray emission and spectral evolution of pair plasmas in AGN jets. I. General theory and a prediction for the GeV - TeV emission from ultrarelativistic jets.  _A &A_ 324, 395–409.  10.48550/arXiv.astro-ph/9604003
  * Chiaberge and Ghisellini (1999) Chiaberge, M. and Ghisellini, G. (1999).  Rapid variability in the synchrotron self-Compton model for blazars.  _MNRAS_ 306, 551–560.  10.1046/j.1365-8711.1999.02538.x
  * Christie et al. (2019) Christie, I. M., Petropoulou, M., Sironi, L., and Giannios, D. (2019).  Radiative signatures of plasmoid-dominated reconnection in blazar jets.  _MNRAS_ 482, 65–82.  10.1093/mnras/sty2636
  * Dermer and Schlickeiser (2002) Dermer, C. D. and Schlickeiser, R. (2002).  Transformation Properties of External Radiation Fields, Energy-Loss Rates and Scattered Spectra, and a Model for Blazar Variability.  _ApJ_ 575, 667–686.  10.1086/341431
  * Di Gesu et al. (2023) Di Gesu, L., Marshall, H. L., Ehlert, S. R., Kim, D. E., Donnarumma, I., Tavecchio, F., et al. (2023).  Discovery of X-ray polarization angle rotation in the jet from blazar Mrk 421.  _Nature Astronomy_ 7, 1245–1258.  10.1038/s41550-023-02032-7
  * Diltz et al. (2015) Diltz, C., Böttcher, M., and Fossati, G. (2015).  Time Dependent Hadronic Modeling of Flat Spectrum Radio Quasars.  _ApJ_ 802, 133.  10.1088/0004-637X/802/2/133
  * Dotson et al. (2015) Dotson, A., Georganopoulos, M., Meyer, E. T., and McCann, K. (2015).  On the Location of the 2009 GeV Flares of Blazar PKS 1510-089.  _ApJ_ 809, 164.  10.1088/0004-637X/809/2/164
  * Finke and Becker (2014) Finke, J. D. and Becker, P. A. (2014).  Fourier Analysis of Blazar Variability.  _ApJ_ 791, 21.  10.1088/0004-637X/791/1/21
  * Finke and Becker (2015) Finke, J. D. and Becker, P. A. (2015).  Fourier Analysis of Blazar Variability: Klein-Nishina Effects and the Jet Scattering Environment.  _ApJ_ 809, 85.  10.1088/0004-637X/809/1/85
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., and Goodman, J. (2013).  emcee: The MCMC Hammer.  _PASP_ 125, 306.  10.1086/670067
  * Joshi and Böttcher (2011) Joshi, M. and Böttcher, M. (2011).  Time-dependent Radiation Transfer in the Internal Shock Model Scenario for Blazar Jets.  _ApJ_ 727, 21.  10.1088/0004-637X/727/1/21
  * Kirk et al. (1998) Kirk, J. G., Rieger, F. M., and Mastichiadis, A. (1998).  Particle acceleration and synchrotron emission in blazar jets.  _A &A_ 333, 452–458.  10.48550/arXiv.astro-ph/9801265
  * Meyer et al. (2021) Meyer, M., Petropoulou, M., and Christie, I. M. (2021).  The Observability of Plasmoid-powered $\gamma$-Ray Flares with the Fermi Large Area Telescope.  _ApJ_ 912, 40.  10.3847/1538-4357/abedab
  * Meyer et al. (2019) Meyer, M., Scargle, J. D., and Blandford, R. D. (2019).  Characterizing the Gamma-Ray Variability of the Brightest Flat Spectrum Radio Quasars Observed with the Fermi LAT.  _ApJ_ 877, 39.  10.3847/1538-4357/ab1651
  * Middei et al. (2023) Middei, R., Perri, M., Puccetti, S., Liodakis, I., Di Gesu, L., Marscher, A. P., et al. (2023).  IXPE and Multiwavelength Observations of Blazar PG 1553+113 Reveal an Orphan Optical Polarization Swing.  _ApJ_ 953, L28.  10.3847/2041-8213/acec3e
  * Nalewajko (2013) Nalewajko, K. (2013).  The brightest gamma-ray flares of blazars.  _MNRAS_ 430, 1324–1333.  10.1093/mnras/sts711
  * Petropoulou et al. (2018) Petropoulou, M., Christie, I. M., Sironi, L., and Giannios, D. (2018).  Plasmoid statistics in relativistic magnetic reconnection.  _MNRAS_ 475, 3797–3812.  10.1093/mnras/sty033
  * Petropoulou and Mastichiadis (2012) Petropoulou, M. and Mastichiadis, A. (2012).  Temporal signatures of leptohadronic feedback mechanisms in compact sources.  _MNRAS_ 421, 2325–2341.  10.1111/j.1365-2966.2012.20460.x
  * Roy et al. (2019) Roy, N., Chatterjee, R., Joshi, M., and Ghosh, A. (2019).  Probing the jets of blazars using the temporal symmetry of their multiwavelength outbursts.  _MNRAS_ 482, 743–757.  10.1093/mnras/sty2748
  * Zacharias (2023) Zacharias, M. (2023).  Exploring the evolution of the particle distribution and the cascade in a moving, expanding emission region in blazar jets.  _A &A_ 669, A151.  10.1051/0004-6361/202244683
  * Zacharias and Schlickeiser (2013) Zacharias, M. and Schlickeiser, R. (2013).  Synchrotron Lightcurves of Blazars in a Time-dependent Synchrotron-self Compton Cooling Scenario.  _ApJ_ 777, 109.  10.1088/0004-637X/777/2/109
  * Zhang et al. (2024) Zhang, H., Böttcher, M., and Liodakis, I. (2024).  Revisiting High-Energy Polarization from Leptonic and Hadronic Blazar Scenarios.  _arXiv e-prints_ , arXiv:2404.1247510.48550/arXiv.2404.12475


A New Optimal Subpattern Assignment (OSPA) Metric  for Multi-target Filtering

This paper proposes and evaluates a new metric.
This metric will overcome a limitation of the Optimal Subpattern
Assignment (OSPA) metric mentioned by Schuhmacher et al.: the
OSPA distance between two sets of points is insensitive to the
the case where one is empty. This proposed metric called Complete
OSPA (COSPA), retains all the advantages of the OSPA metric for
evaluating the performance of multiple target filtering algorithms
while also allowing separate control over the threshold of physical
distance errors and cardinality errors.

Keywords: Optimal Subpattern Assignment (OSPA), Generalized OSPA, Complete OSPA metric (COSPA metric), Multi-target Filtering (MTF).

§ INTRODUCTION

In the area of multi-target filtering (MTF), a rigorous and robust metric plays an important role for performance evaluation. The application of metric can be found in target tracking such as . A key technical issue underpinning such an evaluation concerns the measurement of `distance' between two sets. Clearly a standard distance measure on Euclidean space cannot be used. Hoffman and Mahler [1] studied this problem and found that the Hausdorff distance is relatively insensitive to errors in the number of targets (which is an important issue in MTF) and proposed a new metric to overcome this shortcoming. In $2008$, Schuhmacher et al. [2] built on this work and proposed a metric called optimal sub-pattern assignment (OSPA) which incorporated the spatial distance and the cardinality distance. Spatial distances between pairs of state vectors across two sets are cut off at $c$, a parameter, and are weighted equally while each extra element in a bigger set will be penalized as if there was a distance error of $c$. However, when one set is empty, the OSPA distance takes on the value of $c$ regardless of the cardinality of the other set. This was pointed out in [2] as a minor inconvenience. However, as discussed in Rahmathullah et al. <cit.>, in reality the OSPA metric is not a desirable tool for evaluating MTF algorithms. They pointed out that the MTF community prefers to understand the missing target and false target performance beyond the cardinality mismatch. Hence Rahmathullah et al. [3] tried to overcome these limitations by proposing a new metric called Generalized OSPA (GOSPA). This metric was derived by removing the normalization of the OSPA metric and multiplying parameter $1/\alpha$ to the cardinality error with the optimal choice of $\alpha=2$. However, by not normalizing, the GOSPA metric will generally grow with the size of the sets.

The normalization of the OSPA metric scales the sum of all distances between the two finite sets to be within the interval $[0,c]$. In contrast, the GOSPA metric is exactly the sum of all distances between the finite sets.

Furthermore, by weighting the cardinality penalty less than the spatial distance cut-off $c$, the GOSPA metric gives unexpected results with most scenarios when the sets are not empty (see example in Figure <ref>). Another issue is that the OSPA metric is insensitive to the distance  cut-off and cardinality penalty.
The cause is that the penalty for each unassignable point and the distance cut-off between two assigned points in the OSPA metric are the same, $c$.

By analysing the advantages and disadvantages of OSPA together with the requirements of the MTF community as discussed in [3], we propose a new metric, namely a Complete OSPA (COSPA). This metric overcomes the shortcoming of the OSPA metric when one of the two finite sets is empty. Also, by choosing the cut-off value to be less  than or equal to the penalty for each unassignable point, together with considering the sensitivity to the empty set, the COSPA metric overcomes the limitations of the OSPA metric while still retaining all the advantages of the OSPA metric for other cases when two finite sets are non-empty. The earlier result is published in [4].

The structure of this paper is as follows. Section <ref> provides some definitions and background. Section <ref> summarizes the OSPA metric and discusses its limitations while
Section <ref> summarizes and discusses the GOSPA metric. The main part of this paper is
Section <ref>, where the COSPA metric is developed and the analysis of OSPA, GOSPA and COSPA is conducted,  with illustrative scenarios discussed in Section <ref>.
Section <ref> presents numerical studies to demonstrate the usefulness of the proposed metric and compares this with the other two metrics mentioned above.
Concluding remarks and suggestions for future work are given in Section <ref>.

§ BACKGROUND

Some background and definitions in this section are used throughout the paper.
Let $\mathcal{X}\subseteq\mathds{R}^{n_x}$
where $n_x$ is the natural number. The distance $\bar{d}^{(a)}
(x,y)$ [An arbitrary function with values in $[0,a]$ could be used.
In particular, $\bar{d}^{(a)}
(x_i,y_{\pi(i)})=t(d(x,y))$ where $t:[0,\infty)\to [0,a]$ is any transformation that is nondecreasing, subadditive and $t(u)=0$ if and only if $u=0$.] for $x,y\in\mathcal{X}$ and $a>0$ is cut off by $a$ if it is equal to or larger than $a$. i.e.

\begin{align}\label{E:dinstance_c}
\bar{d}^{(a)}
(x,y)=\min (a,d(x,y))
\end{align}

where $d(x,y)$ is typically $\|x-y\|_{p'}$ where $1\le p'\le \infty$.

A finite set is a collection of a finite number of elements. If $A$ is a finite set, then $A$ has finite number of elements, i.e. $|A|<\infty$. A finite set $B$ is defined to be bigger than a finite set $A$ if $B$ has more elements than $A$, i.e. $|B|>|A|$.

Let $A$ and $B$ be two finite sets. Denote $\Pi(A,B)$ as a set of all one-to-one functions between $A$ and $B$, i.e. if $|A|\le|B|$ then $\pi\in\Pi(A,B)$ is a one-to-one (injective) function from $A$ to $B$ otherwise
from $B$ to $A$. Mathematically, $\Pi(A,B)$ is

\begin{align}\label{E:Permutation}
\hspace{-0.2cm}\left\{\!\!
\begin{array}{ll}
\{\pi\text{ from }A\text{ to }B:\pi \text{ is injective}\}
,\!\!\!&\hbox{if $|A|\le|B|$;} \\
\{\pi\text{ from }B\text{ to } A:\pi \text{ is injective}\},&\hbox{otherwise.}
\end{array}
\right.\!\!
\end{align}

§ OSPA

We summarize the OSPA metric in Section <ref> and discuss its limitations in Section <ref>.

§.§ Summary

The OSPA distance, the total error, is the sum of a localization error and a cardinality error. The localization error is the smallest sum of distances between all combinations of elements of two finite sets $X$ and $Y$. The cardinality error is given by the mismatch between the number of elements in $X$ and the number in $Y$ scaled using $c$.
The per target error obtained by normalizing total error by the largest cardinality of the two given sets is a proper metric.

Let $X$ and $Y$ be two finite sets. For order parameter $p$ ($p\geq 1$), the cut-off parameter $c$ ($c>0$), the OSPA metric $\bar{d}^{(c)}_p(X,Y)=\bar{d}^{(c)}_p(Y,X)$ is defined as.

* if $X=Y=\emptyset$: $\bar{d}^{(c)}_p(X,Y)=0.$

* otherwise if $p<\infty$, $\bar{d}^{(c)}_p(X,Y)$ is for $|X|\le|Y|$

\begin{align}\label{E:OSPA}
%	&\bar{d}^{(c)}_p(X,Y)\nonumber\\=&
\hspace{-0.7cm}
\left[\frac{1}{|Y|}\left(\!\displaystyle\min_{\pi\in\Pi(X,Y)}\!\sum_{x\in X}\!\bar{d}^{(c)}
\end{align}

* if $p=\infty$, $\bar{d}^{(c)}_p(X,Y)$ is

\begin{align}\label{E:OSPA_infty}
%	\!\!\!\!\!\! \bar{d}^{(c)}_p(X,Y)=	\!
\left\{
\begin{array}{ll}
\displaystyle\min_{\pi\in\Pi(X,Y)}\max_{x\in X}\bar{d}^{(c)}
(x,\pi(x)),&\hbox{if $|X|=|Y|$;} \\
c&\hbox{if $|X|\ne |Y|$;}
\end{array}
\right.
\end{align}

$\Pi(X,Y)$ is
in (<ref>); and the choice of $c$ and the steps required to calculate the OSPA metric can be found in [2].

The OSPA distance is interpreted as a $p-$th order per-target error, comprised of a $p-$th order per-target localization error and a $p-$th order per-target cardinality error.
By (<ref>) for $1\le p<\infty$ , the localization and cardinality errors are given respectively as follows for $|X|\le|Y|$

\begin{align*}%\max(m,n)-\min(m,n)
\bar{e}^{(c)}_{p,\text{loc}}(X,Y)&=\bar{e}^{(c)}_{p,\text{loc}}(Y,X)\\&=\left(\frac{1}{|Y|}\min_{\pi\in\Pi(X,Y)}\sum_{x\in X}\!\bar{d}^{(c)}
\bar{e}^{(c)}_{p,\text{card}}(X,Y)&=\bar{e}^{(c)}_{p,\text{card}}(Y,X)=c\left(\frac{|Y|-|X|}{|Y|}\right)^{\frac{1}{p}}\!\!\!\!\!\!.
\end{align*}

The decomposition of the OSPA metric into separate components is usually not necessary for performance evaluation, but may provide valuable additional information.

§.§ Discussion of the OSPA metric

By Definition <ref>, if $X=\emptyset$ and $Y\ne\emptyset$, $\bar{d}^{(c)}_p(X,Y)=c$ regardless of how many elements there are in set $Y$ (see example in Figure <ref> and the discussion in Section <ref>).

Furthermore, the OSPA metric gives the same result for the distance between two sets $X$ and $Y$ where $|Y|=|X|-1$, and between two sets $X$ and $Y\cup\{z\}$ where $d(x,z)\ge c$ for $z\notin Y$ and $\forall x\in X$ (see example in Figure <ref>). The GOSPA metric [3] overcomes these limitations but it has shortcomings compared with OSPA. The next section summarizes and discusses the GOSPA metric.

§ GENERALIZED OSPA (GOSPA)

The GOSPA metric [3] is summarized in Subsection <ref>. The GOSPA metric is an unnormalized version of OSPA with some scaling of the cardinality parameter.
GOSPA actually helps to overcome the limitation of the  OSPA metric when one of the two sets is empty. However, the lack of normalization causes the GOSPA metric to increase when the sizes of the sets increase.
These limitations are discussed in Subsection <ref>.

§.§ Summary

Let $X$ and $Y$ be two finite sets.
For $p\geq 1$, $c>0,$ $0<\alpha\le 2$, the GOSPA metric,  $\bar{d}^{(c,\alpha)}_p(X,Y)=\bar{d}^{(c,\alpha)}_p(Y,X)$, is defined as follows for $|X|\le|Y|$,

\begin{eqnarray}\label{E:GOSPA}
\left(\displaystyle\min_{\pi\in\Pi(X,Y)}\sum_{x\in X}\!\bar{d}^{(c)}
(x,\pi(x))^p+\frac{c^p}{\alpha }(|Y|-|X|)\right)^{\frac{1}{p}}.&
\end{eqnarray}

Let $X$ and $Y$ be two finite sets.
For $p\geq 1$, $c>0$, $\alpha=2$, the GOSPA metric,  $\bar{d}^{(c,2)}_p(X,Y)=\bar{d}^{(c,2)}_p(Y,X)$, can be expressed as an optimization over assignment sets as for $|X|\le|Y|$,

\begin{align}\label{E:PropGOSPA}
&\Bigg(\displaystyle\min_{\pi\in\Pi(X,Y)}\Bigg[\sum_{x\in X}\bar{\delta}\big(\bar{d}^{(c)}
(x,\pi(x))^p\nonumber\\&+\frac{c^p}{2}\bigg(|X|+|Y|-2\sum_{x\in X}\bar{\delta}\big(\bar{d}^{(c)}
\end{align}

where $\delta(k,h)$ is a Kronecker delta (i.e. $\delta(k,h)$ is $1$ if $k=h$ or $0$ if $k\ne h$) and $\bar{\delta}(k,h)=1-\delta(k,h)$.

§.§ Discussion of the GOSPA metric

The GOSPA metric overcomes
the limitations of the OSPA metric by removing the normalization of the OSPA metric and multiplying the second term of (<ref>) with $1/\alpha$ where $0<\alpha\le 2$. However, these alterations to the OSPA metric have an unexpected side-effect: GOSPA distance between two finite sets in general increases if one or both of the two sets increases in cardinality (details in Section <ref>).
Furthermore, the cut-off $c$, being larger than the penalty given to each extra vector  $c/\sqrt[p]{\alpha}$ if $\alpha\ge 1$, makes the GOSPA metric even less useful (see Figure <ref>).  Definition <ref> shows that it is obviously that $\bar{d}^{(c,2)}_p(X,Y)\le \bar{d}^{(c,\alpha)}_p(X,Y)$ for all
$0<\alpha\le 2$. Hence, it is not clear what an optimization over assignments sets mean in the Proposition <ref> but incidentally, Proposition <ref> is just a way to group the mismatched targets as the second term in (<ref>) when $\alpha=2$ (the mismatched targets are two targets from two different sets that are further than or equal to $c$). If $\alpha\ne 2$, we cannot rewrite (<ref>) in the form of (<ref>).

In order to overcome the limitations of the OSPA metric discussed in Section <ref>, we propose a new metric, namely Complete OSPA (COSPA), in the next section. This metric not only overcomes the above limitations but also retains the advantages of the OSPA metric for evaluating cases where the OSPA metric gives reliable solutions.

§ COMPLETE OPTIMAL SUBPATTERN ASSIGNMENT (COSPA) METRIC

The Optimal Subpattern Assignment (OSPA) metric is very popular in multi-target filtering because it can evaluate the localization error and cardinality error between two sets of vectors. However, as discussed in Section <ref>
the OSPA metric between an empty set and a non-empty set is the same regardless of how many elements the non-empty set has. It is also indistinguishable between the cut-off distance and the cardinality penalty
and discussed in Section <ref>. Hence, in this section, we will derive a new metric which overcomes these shortcomings and still retains all the beneficial properties of the OSPA metric. The new metric, which is called Complete OSPA (COSPA), cuts off the distance between two vectors at
$c$ ($c>0$), penalizes each extra of points in a bigger set by $\dot{c}$ ($\dot{c}\ge c$) and penalized the empty set error by $\xi$ ($0\le\xi\le 1$). In the rest of this section, we show that COSPA will overcome the shortcomings of the current metrics OSPA and COSPA.

Let $X$ and $Y$ be two finite sets. For order parameter $p$ ($1\le p\le\infty$), cut-off parameter $c$ ($0<c$) and cardinality penalty parameter $\dot{c}$ ($\dot{c}\ge c$)  and empty-set parameter $0\le\xi\le 1$, the COSPA metric $d^{(\dot{c},\xi)}_{c,p}(Y,X)=d^{(\dot{c},\xi)}_{c,p}(X,Y)$ is defined as follows.

* If $X=Y=\emptyset$: $d^{(\dot{c},\xi)}_{c,p}(X,Y)=0.$

* Otherwise, for $p<\infty$, $d^{(\dot{c},\xi)}_{c,p}(X,Y)$ is for $|X|\le|Y|$

\begin{eqnarray}
\Bigg(\displaystyle
\dfrac{1}{|Y|}\bigg[\min_{\pi\in\Pi(X,Y)}\sum_{x\in X}\Big(c^p\delta\big(\bar{d}^{(c)}
\label{E:DiffCOSPA1}\\
& +\xi\dfrac{\delta(X,\emptyset)\bar{\delta}(Y,\emptyset)}{|Y|}\displaystyle \dot{c}^p(|Y|-1)\Bigg)^{\frac{1}{p}}\label{E:COSPA2}
\end{eqnarray}

* if $p=\infty$, $d^{(\dot{c},\xi)}_{c,\infty}(X,Y)$ is

\begin{align}\label{E:COSPA_infty}
\left\{
\begin{array}{ll}
\displaystyle\min_{\pi\in\Pi(X,Y)}\max_{x\in X}\bar{d}^{(c)}
(x,\pi(x)) &\!\!\! \hbox{if $|X|=|Y|$} \\
\dot{c}&  \hbox{if $|X|\ne |Y|$}
\end{array}
\right.
\end{align}

The function $d^{(\dot{c},\xi)}_{c,p}(\cdot,\cdot)$ is called the COSPA metric of order $p$ with cut-off $c$, cardinality penalty $\dot{c}$ and empty-set error $\xi$.

COSPA is proved as a metric in [4] (note that $\dot{c}$ and $c$ are swapped in [4]). If we choose $\dot{c}=c$ and $\xi=0$, COSPA is exactly OSPA. The term $\dot{c}$, which is the penalty for each cardinality error and differs from the cut off distance $c$, exists to overcome the limitation of OSPA for scenarios that distinguish between the cut-off distance from a cardinality error and are similar to the example shown in Figure <ref>.
In this example the OSPA distance between the two sets in Figure (<ref>) is the same as the OSPA distance between the two sets in Figure (<ref>). In general, the OSPA distances for two pairs of finite sets are the same when the cardinality of a set of the first pair is smaller than the other set of that pair and the second pair is the same as the first pair with the smaller set now having an extra element whose distance to elements in the other set of that pair is larger than $c$.

Note that the cut-off parameter $c$ is smaller than cardinality penalty $\dot{c}$ if the outline distance between two vectors is penalized less than each mismatched number of elements between two sets of vectors. If $\dot{c}=c$, Definition <ref> (COSPA) is only different from Definition <ref> (OSPA) [2] by the term in (<ref>). The term in (<ref>)  only exists to take into account the case when one of the two sets is empty. When both finite sets are non-empty the term (<ref>) is zero, so COSPA is exactly the same as OSPA [2] if $\dot{c}=c$.

If we do not wish to distinguish what vectors in $X$ are very far from their images in $Y$,
(x,\pi(x))^p$ for
Alternatively, $d^{(\dot{c},\xi)}_{c,p}(X,Y)$ in (<ref>) can be written for simplicity as follows if $|X|\le|Y|$ and $Y\ne\emptyset$.

\begin{eqnarray}
&\bigg[\displaystyle\frac{1}{|Y|}\bigg(\min_{\pi\in\Pi(X,Y)}\!\sum_{x\in X}\!\bar{d}^{(c)}
\!(x,\pi(x))^p+\dot{c}^p(|Y|-|X|)\!\bigg)\label{E:AltCOSPA1}\\
& %\hspace{.15cm}
+\xi\dfrac{\delta(X,\emptyset)\bar{\delta}(Y,\emptyset)}{|Y|}\displaystyle \dot{c}^p(|Y|-1)\bigg]^{\frac{1}{p}}\label{E:AltCOSPA2}
\end{eqnarray}

From Definition <ref>, Definition <ref> and (<ref>), the following are true

* If $X\ne\emptyset$ and $Y\ne\emptyset$, (<ref>) shows that the COSPA distance between two non-empty finite sets is smaller than or equal to $\dot{c}$. In this case, the COSPA distance (<ref>) has the same forms as the OSPA distance (<ref>) because the term in (<ref>) does not exist. Furthermore, if $\dot{c}=c$, then the COSPA distance is the same as the OSPA distance, i.e. $d^{(\dot{c},\xi)}_{c,p}(X,Y)=\bar{d}^{(c)}_p(X,Y).$

* If $X=\emptyset$ or $Y=\emptyset$ but not both,

\begin{align*}
\end{align*}

By (<ref>),
the COSPA assignment is the OSPA assignment. The COSPA assignment between $X$ and $Y$, $\pi^*$, is

\begin{align}\label{E:COSPA_Ass}
\begin{array}{ll}
\displaystyle\arg\min_{\pi\in\Pi(X,Y)}\sum_{x\in X}\bar{d}^{(c)}
\!(x,\pi(x))^p
,\!\!\!&\hbox{if  $|X|\le|Y|$;} \\
\displaystyle\arg\min_{\pi\in\Pi(X,Y)}\sum_{y\in Y}\bar{d}^{(c)}
\!(y,\pi(y))^p,&\hbox{otherwise.}
\end{array}
\end{align}

Assume that $X$ is the set of truth targets[In Remark <ref>, the term `target' and 'vector' are used interchangeably.] and $Y$ is the set of estimated targets. Definition <ref> can be interpreted as follows.

* The term
in (<ref>) is actually the sum of distances between vectors in $X$ and their images in $Y$ by the one to one function $\pi^*$ if each of these distances is smaller than $c$. Each pair $(x,\pi^*(x))\in X\times Y$ is a pair of correctly associated targets (vectors) in $X$ and $Y$ if their distance is smaller than $c$. Note that $\pi^*\in\Pi(X,Y)$ is defined in (<ref>).

* The term in (<ref>) is actually the sum of $|\gamma|$ distances between $|\gamma|$ vectors in $X$ and the $|\gamma|$  vectors in $Y$ by the one to one function $\pi^*$ when each of these distances is larger than or equal to $c$ where

\begin{align*}%\label{E:dinstance_c}
\gamma=\{(x,\pi^*(x))\in X\times Y:\bar{d}^{(c)}
\end{align*}

Each point in the map $\gamma$ is a pair of incorrectly associated targets if their distance is equal to or larger than $c$.
Alternatively, the target in $X$
is called a missing target and the
corresponding estimated target in $Y$ via $\pi^*$
is called a false target.

* The term in (<ref>) is actually the cardinality error between $X$ and $Y$. If $|Y|>|X|$, the targets in $Y$ which do not associate with any targets in $X$ via $\pi^*$ are called false targets.

* The term in (<ref>) only exists if the smaller set is empty and the bigger set is not empty.

The COSPA distance comprises of 3 components: COSPA localization error, COSPA outline error (i.e. the sum of all distances that are larger than or equal to $c$) and COSPA cardinality error. Here, the OSPA localization error is the sum of the COSPA localization error and the COSPA outline error. The COSPA cardinality error is the same as the OSPA cardinality error if both sets are not empty. If either $X$ or $Y$ is empty, the COSPA cardinality error is the COSPA distance (Remark <ref>.2) while the OSPA cardinality error is the OSPA distance $c$. Hence, similar to the OSPA metric, the COSPA metric is interpreted as a $p-$th order per-target error.
By (<ref>) for $1\le p<\infty$, COSPA localization (<ref>), COSPA outline (<ref>) and cardinality errors (<ref>) are given respectively as follows for $|X|\le|Y|$

\begin{align}
\bigg(\frac{1}{|Y|}\sum_{x\in X}\bar{\delta}\big(\bar{d}^{(c)}
&\bar{e}^{(c)}_{p,\text{out}}(X,Y)=\bar{e}^{(c)}_{p,\text{out}}(Y,X)\nonumber\\=&c\left(\frac{\sum_{x\in X}\delta\big(\bar{d}^{(c)}
\end{align}

§ ANALYSIS OF COSPA, GOSPA AND OSPA METRICS

In  this section, we show the solution the COSPA metric offers to the limitations of OSPA and GOSPA mentioned in Sections <ref>  and <ref> via simple scenarios. Furthermore, we also analyze the scenarios by comparing these solutions. Without loss of generality $\xi$ in Definition <ref> is chosen as $1$ (i.e. $\xi=1$) and the order parameter $p=1$ for the three metrics.

§.§ Effect of Cardinality Zero

As discussed in Section <ref>, the OSPA metric gives the same result when one of the two arguments is empty. In Figure <ref> the OSPA distance between the non-empty set $Y$ and $\emptyset$
is the same as the OSPA distance between the non-empty set $Z$ and $\emptyset$. The GOSPA and COSPA metrics give smaller values for  $(\emptyset, Y)$ than $(\emptyset, Z)$. As a result, a set with two vectors $Y$ is closer to $\emptyset$ than a set with three vectors $Z$. This is a natural (meaningful) physical interpretation. Detailed explanations are given in Table <ref>.

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0.2) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=3pt] (x1);

(lx1)at (0,-.2)$y_1$;
[circle,draw=black, fill=black, inner sep=0pt,minimum size=3pt](x2)[right =1.7cm of x1] ;
(x3)[right =3.2cm of x1] ;

(lx3)at (3.35,-.2);
(lx2)at (1.85,-.2)$y_2$;

(lx3.east |-x1.north)+(.05,+.2) node(cor1);
(lx1.west|-lx1.south)+(-.05,-.1) node(cor2);
(cor1) rectangle  (cor2);

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0.2) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=3pt] (x1);

(lx1)at (0,-.2)$z_1$;
[circle,draw=black, fill=black, inner sep=0pt,minimum size=3pt](x2)[right =1.7cm of x1] ;
[circle,draw=black, fill=black, inner sep=0pt,minimum size=3pt](x3)[right =3.2cm of x1] ;

(lx3)at (3.3,-.2)$z_3$;
(lx2)at (1.8,-.2)$z_2$;

(lx3.east |-x1.north)+(.05,+.2) node(cor1);
(lx1.west|-lx1.south)+(-.05,-.1) node(cor2);
(cor1) rectangle  (cor2);

The OSPA metric $\bar{d}^{(c)}_p(\emptyset,Y)=\frac{2c}{2}=c= \frac{3c}{3}=\bar{d}^{(c)}_p(\emptyset,Z)$. OSPA assigns the same value to the distance between any non-empty set and an empty set.

Analysis of the three metrics for Figure <ref>

Metric   Fig.<ref> ($Y$)   Fig. <ref> ($Z$)   Which one is closer to $\emptyset$?
OSPA    $c$     $c$    They are the same
GOSPA   $c\frac{2}{\alpha}$    $c\frac{3}{\alpha}$   $Y$
COSPA   $\dot{c}\frac{3}{2}$    $\dot{c}\frac{5}{3}$   $Y$
Intuition[`Intuition' usually corresponds to a surveillance operator's likely needs, i.e. the cardinality penalty is larger than or equal to the distance cut off due to $\dot{c}\ge c$.]            $Y$

The COSPA metric is different from the other two metrics when comparing one algorithm producing an empty set and another algorithm giving a non-empty set as an output. Compared with the non-empty set ground truth, COSPA gives an algorithm with the empty set a bigger value than an algorithm with the non-empty set. The distance between two empty sets is zero. It is clear that the three metrics will give the algorithm with the empty set a smaller distance when the ground truth is empty. Figure <ref> and Table <ref> give an example when the ground truth is not empty. GOSPA  concludes that the empty set $Z$ is the closer estimate to $X$ than the non empty set $Y$ if $\frac{c}{\alpha}\le\eta$ and $\alpha>1$. These results depend on the choice of cut-off $c$. If $c\le\eta$, OSPA concludes that $Y$ is as good as $Z$ for estimating $X$, otherwise $Y$ is the better estimate of $X$ than $Z$. The COSPA metric concludes that $Y$ is the better estimate of $X$ than $Z$ no matter what values are chosen for the cut-off $c$ and cardinality  error penalty $\dot{c}$ (see more detail in Table <ref>). This is a
natural (meaningful) physical interpretation.

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x1);
(0,-.3) node $x_1$;

(0.8,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x2) ;
(0.8,-.3)node $x_2$;

(1.6,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xi) ;
(1.6,-.3) node $x_i$;

(2.4,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xn-1) ;

(3.2,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xn) ;
(3.3,-.3) node $x_m$;

(xn.east |-x1.north)+(.35,+.25) node(cor1);
(x1.west|-x1.south)+(-.35,-.5) node(cor2);
(cor1) rectangle  (cor2);

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x1);
(0,-.3) node $x_1$;

(0.8,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x2) ;
(0.8,-.3)node $x_2$;

(1.6,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xi) ;
(1.6,-.3) node $x_i$;

(2.4,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xn-1) ;

(3.2,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xn) ;
(3.3,-.3) node $x_m$;

(0,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y1);
(0,1.2) node $y_1$;

(.8,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y2);
(.8,1.2) node $y_2$;

(1.6,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](yi);
(1.6,1.2) node $y_i$;

(2.4,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](yn-1) ;

(3.2,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](yn);
(3.3,1.2) node $y_m$;

(temdis11)  [above =.33cm of x1] ;
(dis2) [right =-.15cm of temdis11]$\eta$;

[-,thin,black,dotted](x1)edge node[swap,near start]  (y1);

(yn.east |-yn.north)+(.5,+.45) node(cor1);
(x1.west|-x1.south)+(-.35,-.5) node(cor2);
(cor1) rectangle  (cor2);

The ground truth is $X=\{x_1,\ldots,x_m\}$, $m>1$. In Figure <ref>, an algorithm estimates no target, i.e. $Z=\emptyset$, while an algorithm in Figure <ref> estimates $m$ targets, i.e $Y=\{y_1,y_2,\ldots,y_{m}\}$.

Analysis of the three metrics
for Figure <ref>

Metric   Fig. <ref>($Z$)   Fig. <ref> ($Y$)   Which one
is closer to $X$?
OSPA    $c$     $\min(c,\eta)$
$Y$ if $\boldsymbol{\eta<c}$

GOSPA   $c\frac{m}{\alpha}$    $m\min(c,\eta)$

$Y$  if $\boldsymbol{\eta<\frac{c}{\alpha}}$

COSPA   $\dot{c}\left(2-\frac{1}{m}\right)$     $\min(c,\eta)
$   $Y$
Intuition           $Y$

§.§ Effect of Choice of Cut-off and Cardinality  Penalty

In the OSPA metric, the cut-off
distance  $c$ between two vectors is the same as the penalty for each extra vector. For the COSPA metric, this cut-off distance $c$
is smaller than or equal to the penalty for each extra vector, $\dot{c}$. This is not the case for the GOSPA metric (<ref>) where the cut-off distance $c$ between two vectors is larger than or equal to the penalty for each extra vector, $c/\sqrt[p]{\alpha}$.
This choice makes the GOSPA metric intuitively unreliable[If a metric is `unreliable', it occasionally assign a large value to the distance between two sets of vectors that are intuitively close.] for most other scenarios where both sets are not empty.
GOSPA gives the same value for the two scenarios in Figure <ref> if $\alpha = 1$; the smaller or bigger value for scenario <ref> than <ref> if $\alpha < 1$ or $\alpha > 1$ respectively. If $\alpha>1$, which is the preferred distance choice of the author, GOSPA gives the smaller, same or bigger distance for scenario <ref> than <ref> if $c> \Delta> c/\alpha$, $\Delta= c/\alpha$, or $\Delta<c/\alpha$ respectively.

As long as $\Delta<c$, OSPA gives the same value for the two scenarios in Figure <ref> while COSPA gives a smaller value for scenario <ref> than <ref> if $\dot{c}>c$ and the same value if $\dot{c}=c$. This is summarized in Table <ref>. Both OSPA and COSPA give the smaller value for scenario <ref> than for <ref>, agreeing with intuitive thinking. This is shown in Table <ref>. Intuitively, $Y_c$ is the closest to $X$ as given by OSPA and COSPA. However, the GOSPA metric only gives the same conclusion if $\Delta<c/\alpha$. If $\alpha=2$, which the authors [3] claim is the optimal solution for the GOSPA metric, GOSPA gives  $Y_a$ as the closest to $X$ if $\Delta\ge c/\alpha$, which contradicts intuitive thinking. The behaviour of these three metrics for the scenarios in Figure <ref> are summarized in Table <ref>.

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x1);
(0,-.3) node $x_1$;

(0.6,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x2) ;
(0.6,-.3)node $x_2$;

(1.2,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x3) ;
(1.2,-.3) node $x_3$;

(0,.4) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y1);
(0,.6) node $y_1$;

(.6,.4) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y2);
(.6,.6) node $y_2$;

(temdis11)  [above =.05cm of x1] ;
(dis2) [right =-.15cm of temdis11]$\Delta$;

[-,thin,black,dashed](x1)edge node[swap,near start]  (y1);
[-,thin,black,dashed](x2)edge node[swap,near start]  (y2);

(x3.east |-y2.north)+(.5,+.45) node(cor1);
(x1.west|-x1.south)+(-.35,-.5) node(cor2);
(cor1) rectangle  (cor2);

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x1);
(0,-.3) node $x_1$;

(0.6,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x2) ;
(0.6,-.3)node $x_2$;

(1.2,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x3) ;
(1.2,-.3) node $x_3$;

(0,.4) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y1);
(0,.6) node $y_1$;

(.6,.4) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y2);
(.6,.6) node $y_2$;

(2.8,.4) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y2n1);
(3,.6) node $z$;

(temdis11)  [above =.05cm of x1] ;
(dis2) [right =-.15cm of temdis11]$\Delta$;

[-,thin,black,dashed](x1)edge node[swap,near start]  (y1);
[-,thin,black,dashed](x2)edge node[swap,near start]  (y2);
[-,thin,black,dashed](x3)edge node[swap,near start]  (y2n1);

(y2n1.east |-y2n1.north)+(.5,+.45) node(cor1);
(x1.west|-x1.south)+(-.35,-.5) node(cor2);
(cor1) rectangle  (cor2);

$X=\{x_1,x_2,x_3\}$ where $d(x_i,z)\ge c>\Delta$ for $i=1,2,3$.
The OSPA metric $\bar{d}^{(c)}_p(X,Y)=\frac{2\Delta+c}{3}= \bar{d}^{(c)}_p(X,Z)$.
OSPA gives the same value for the distance between $X$ and $Y$ and for the distance between $X$ and $Z$.

Analysis of the three metrics
for	Figure <ref>

Metric   Fig.<ref> ($Y$)   Fig.<ref> ($Z$)   Which one ($Y$ or $Z$) is closer to $X$?   Explanation
OSPA    $\frac{2\Delta^p+c^p}{3}$     $\frac{2\Delta+c}{3}$   They have the same distance to $X$    cut-off $=$ cardinality penalty
GOSPA   $2\Delta+\frac{1}{\alpha}c$    $2\Delta+c$   $Y$ if $\alpha>1$; $Z$ if $\alpha<1$    cut-off $\ne$ cardinality penalty ($c\ne\frac{c}{\alpha}$)
COSPA   $\frac{2\Delta+\dot{c}}{3}$    $\frac{2\Delta+c}{3}$   $Z$ if $c<\dot{c}$    cut-off $\le$ cardinality penalty ($c\le\dot{c}$)

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0.2) node[cross] (x1);

(lx1)at (0,-.2)$x_1$;
(x3)[right =.7cm of x1] ;

(1.1,.2) node[cross] (x2);
[above =2cm of x3] (y3) ;

(y2) at (1.1,1)
(ly3) at (1,2.7);
(ly2) at (1.1,1.3);
(lx2)at (1.1,-.2)$x_2$;

(ly2.east |-ly3.north)+(.14,+.12) node(cor1);
(lx1.west|-lx1.south)+(-.01,-.1) node(cor2);
(cor1) rectangle  (cor2);

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0.2) node[cross] (x1);

[circle,draw=black, fill=black, inner sep=0pt,minimum size=3pt] (y1) at (0,1) ;
(ly1) at (0,1.3)$y_1$;
(lx1)at (0,-.2)$x_1$;
(x3)[right =.7cm of x1] ;

(1.1,.2) node[cross] (x2);
[above =2cm of x3] (y3) ;

(y2) at (1.1,1)
(ly3) at (1,2.7);
(ly2) at (1.1,1.3);
(lx2)at (1.1,-.2)$x_2$;

(temdis11)  [below =.2cm of y1] ;
(dis1) [right =.05cm of temdis11] $\Delta$;
[-,thin,black,dashed](x1)edge node[swap,near start]  (y1);

(ly2.east |-ly3.north)+(.14,+.12) node(cor1);
(lx1.west|-lx1.south)+(-.01,-.1) node(cor2);
(cor1) rectangle  (cor2);

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0.2) node[cross] (x1);

[circle,draw=black, fill=black, inner sep=0pt,minimum size=3pt] (y1) at (0,1) ;
(ly1) at (0,1.3)$y_1$;
(lx1)at (0,-.2)$x_1$;
(x3)[right =.7cm of x1] ;

(1.1,.2) node[cross] (x2);
[above =2cm of x3] (y3) ;

[circle,draw=black, fill=black, inner sep=0pt,minimum size=3pt] (y2) at (1.1,1)
(ly3) at (1,2.7);
(ly2) at (1.1,1.3)$y_2$;
(lx2)at (1.1,-.2)$x_2$;

(temdis11)  [below =.2cm of y1] ;
(dis1) [right =.05cm of temdis11] $\Delta$;
[-,thin,black,dashed](x1)edge node[swap,near start]  (y1);

(temdis)  [below =.2cm of y2] ;
(dis2) [right =.05cm of temdis];
[-,thin,black,dashed](x2)edge node[swap,near start]  (y2);

(ly2.east |-ly3.north)+(.14,+.12) node(cor1);
(lx1.west|-lx1.south)+(-.01,-.1) node(cor2);
(cor1) rectangle  (cor2);

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0.2) node[cross] (x1);

[circle,draw=black, fill=black, inner sep=0pt,minimum size=3pt] (y1) at (0,1) ;
(ly1) at (0,1.3)$y_1$;
(lx1)at (0,-.2)$x_1$;
(x3)[right =.7cm of x1] ;

(1.1,.2) node[cross] (x2);
[circle,draw=black, fill=black, inner sep=0pt,minimum size=3pt,above =2cm of x3] (y3) ;

[circle,draw=black, fill=black, inner sep=0pt,minimum size=3pt] (y2) at (1.1,1)
(ly3) at (1,2.7)$y_3$;
(ly2) at (1.1,1.3)$y_2$;
(lx2)at (1.1,-.2)$x_2$;

(temdis11)  [below =.2cm of y1] ;
(dis1) [right =.05cm of temdis11] $\Delta$;
[-,thin,black,dashed](x1)edge node[swap,near start]  (y1);

(temdis)  [below =.2cm of y2] ;
(dis2) [right =.05cm of temdis];
[-,thin,black,dashed](x2)edge node[swap,near start]  (y2);

(ly2.east |-ly3.north)+(.01,+.01) node(cor1);
(lx1.west|-lx1.south)+(-.01,-.1) node(cor2);
(cor1) rectangle  (cor2);

$\Delta<c$, $X=\{x_1,x_2\}$ and the distances $d(x_1,y_1)=d(x_2,y_2)=\Delta$ and $d(x_i,y_3)>\Delta$ for $i=1,2$.

Analysis of the three metrics and intuitive thinking to evaluate the sets shown in Figure <ref>

Metric   Fig.<ref> ($Y_a$)   Fig.<ref> ($Y_b$)   Fig. <ref> ($Y_c$)   Fig.<ref>  ($Y_d$)   Which one is the closest to $X$?
OSPA    $c$     $\frac{c+\Delta}{2}$    $\Delta$   $\frac{2\Delta+c}{3}$   $Y_c$
GOSPA   $c\frac{2}{\alpha}$    $\Delta+\frac{c}{\alpha}$   $2\Delta$   $2\Delta+\frac{c}{\alpha}$

$Y_c$ if $\boldsymbol{\Delta< c\frac{1}{\alpha}}$ or $Y_a$ if $\boldsymbol{\Delta\ge c\frac{1}{\alpha}}$

COSPA   $\dot{c}\frac{3}{2}$     $\frac{\dot{c}+\Delta}{2}$    $\Delta$   $\frac{2\Delta+\dot{c}}{3}$   $Y_c$
Intuition                   $Y_c$

§.§ Importance of Normalization

The normalization of the distance between two finite sets plays an important role for measuring how close these two finite sets are because it scales the total error, which is the minimum sum of all distances of pairs of vectors and all distances of unpaired vectors, to be within the interval $[0,c]$ for the OSPA metric and to  be within the half-closed interval $[0,\dot{c}+\dot{c}\xi)$ where $0\le\xi\le 1$  for the COSPA metric. Without the normalization, the GOSPA metric is simply the sum. Hence, if one/both of these sets are large and distinguishable, the total error of these two sets is large and hence the GOSPA distance is large. Considering Figures <ref> and <ref>, the GOSPA metric concludes that $Y_d$ is worse than $Y_b$ because their distances to the same set $X$ gives the bigger value to the first compared to the second. This is summarized in Table <ref>. Clearly, the unnormalized OSPA distance is not a proper distance measure because it is proportional to the size of the bigger set. Therefore, the GOSPA metric may not be a suitable tool to measure how close two finite sets of vectors are.

Now consider Figure <ref> as an example. There are two  parallel line segments in Figure <ref> and the distance between these line segments is $\Delta$. If these two line segments are discretized into two sets of two points each, with $X=\{x_1,x_2\}$ and $Y=\{y_1,y_2\}$ as in Figure <ref>, the GOSPA distance will be $\Delta\sqrt[p]{2}$. Similarly, if these two lines are discretized  into two sets of $n$ points each, $X=\{x_1,x_2,\ldots,x_n\}$ and $Y=\{y_1,y_2,\ldots,y_n\}$ as in Figure <ref>, then the GOSPA metric will be $\Delta\sqrt[p]{n}$. This means the distance for the scenario in Figure <ref> is larger than the distance for the scenario in Figure <ref>. This is not intuitively reasonable because all scenarios in Figure <ref> have the same distance, $\Delta$, between the two line segments $l_1$ and $l_2$ in Figure <ref>. OSPA and COSPA give the same distance $\Delta$ for all scenarios in Figure <ref>. The explanation and computation of the three metrics for Figures <ref>-<ref> are summarized in Table <ref>.

Analysis of the three metrics
for Figures <ref> - <ref>

Metric   Fig.<ref>  ($Y_b$)   Fig.<ref> ($Y_d$)   Which one
is closer to $X$?
OSPA     $\frac{c+\Delta}{2}$    $\frac{2\Delta+c}{3}$   $Y_d$
GOSPA    $\Delta+\frac{c}{\alpha}$   $2\Delta+\frac{c}{\alpha}$   $Y_b$
COSPA    $\frac{\dot{c}+\Delta}{2}$    $\frac{2\Delta+\dot{c}}{3}$   $Y_d$
Intuition           $Y_d$

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt]
(0,-.2) node $l_1$;

(0,2) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x2) ;

(1.1,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y1);
(1.1,-.2) node $l_2$;

(1.1,2) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y2);

(temdis11)  [above =.2cm of x1] ;
(dis2) [right =.3cm of temdis11] ;
(dis1) [below =-.1cm of dis2]$\Delta$;
[-,thin,black,dotted](x1)edge node[swap,near start]  (y1);

[-,thick,black](x1)edge node[swap,near start]  (x2);
[-,thick,black](y1)edge node[swap,near start]  (y2);

(y1.east |-x2.north)+(.25,+.35) node(cor1);
(x1.west|-x1.south)+(-.25,-.35) node(cor2);
(cor1) rectangle  (cor2);

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x1);
(0,-.2) node $x_1$;

(0,2) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x2) ;
(0,2.15) node $x_2$;

(1.1,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y1);
(1.1,-.2) node $y_1$;

(1.1,2) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y2);
(1.1,2.15) node $y_2$;

(temdis11)  [above =.2cm of x1] ;
(dis2) [right =.3cm of temdis11] ;
(dis1) [below =-.1cm of dis2]$\Delta$;
[-,thin,black,dotted](x1)edge node[swap,near start]  (y1);

[-,thick,black,dashed](x1)edge node[swap,near start]  (x2);
[-,thick,black,dashed](y1)edge node[swap,near start]  (y2);

(y2.east |-x2.north)+(.2,+.35) node(cor1);
(x1.west|-x1.south)+(-.25,-.35) node(cor2);
(cor1) rectangle  (cor2);

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x1);
(0,-.2) node $x_1$;

(0,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x2) ;
(-0.25,1) node $x_2$;

(0,2) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x3) ;
(0,2.15) node $x_3$;

(1.1,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y1);
(1.1,-.2) node $y_1$;

(1.1,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y2);
(.85,1) node $y_2$;

(1.1,2) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y3);
(1.1,2.15) node $y_3$;

(temdis11)  [above =.2cm of x1] ;
(dis2) [right =.3cm of temdis11] ;
(dis1) [below =-.1cm of dis2]$\Delta$;
[-,thin,black,dotted](x1)edge node[swap,near start]  (y1);

[-,thick,black,dashed](x1)edge node[swap,near start]  (x3);
[-,thick,black,dashed](y1)edge node[swap,near start]  (y3);

(y3.east |-x3.north)+(.2,+.35) node(cor1);
(x2.west|-x1.south)+(-.45,-.35) node(cor2);
(cor1) rectangle  (cor2);

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x1);
(0,-.2) node $x_1$;

(0,.4) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x2) ;
(-0.25,.4) node $x_2$;

(0,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xi) ;
(-0.25,1) node $x_i$;

(0,1.68) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xn-1) ;
(0,1.68) node ;

(0,2) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xn) ;
(0,2.15) node $x_n$;

(1.1,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y1);
(1.1,-.2) node $y_1$;

(1.1,.4) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y2);
(.85,.4) node $y_2$;

(1.1,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](yi);
(.85,1) node $y_i$;

(1.1,1.68) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](yn-1) ;
(1.1,1.68) node ;

(1.1,2) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](yn);
(1.1,2.15) node $y_n$;

(temdis11)  [above =.2cm of x1] ;
(dis2) [right =.3cm of temdis11] ;
(dis1) [below =-.1cm of dis2]$\Delta$;
[-,thin,black,dotted](x1)edge node[swap,near start]  (y1);

[-,thick,black,dashed](x1)edge node[swap,near start]  (xn);
[-,thick,black,dashed](y1)edge node[swap,near start]  (yn);

(yn.east |-xn.north)+(.2,+.35) node(cor1);
(xi.west|-x1.south)+(-.45,-.35) node(cor2);
(cor1) rectangle  (cor2);

The Euclidean distance between two line segments $l_1$ and $l_2$ is $\Delta$ which is shown in Figure <ref>. If these two line segments are discretized into two sets of two points each as in Figure <ref>; two sets of three points each as in Figure <ref>; and two sets of $n$ points each as in Figure <ref>, then the GOSPA metric is $2\Delta$, $3\Delta$ and $n\Delta$ respectively. The OSPA and COSPA metrics are the same and equal to $\Delta$ for all cases.

Analysis of the three metrics
for Figures <ref>-<ref>

Metric   Fig.<ref>   Fig.<ref>   Fig.<ref>   Are the distances
the same?
OSPA    $\Delta$    $\Delta$    $\Delta$   Yes
GOSPA   $2\Delta$    $3\Delta$   $n\Delta$   No
COSPA   $\Delta$    $\Delta$    $\Delta$   Yes
Intuition             Yes

Furthermore, the absence of normalization makes GOSPA inconsistent when comparing the empty set with a non-empty set using another non-empty set as a reference. Take Figure <ref> as an example. The computation and the comparison of the results (of the three metrics) against intuitive thinking are summarized in Table <ref>.

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x1);
(0,-.3) node $x_1$;

(0.4,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x2) ;
(0.4,-.3)node $x_2$;

(1,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xi) ;
(1,-.3) node $x_i$;

(1.68,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xn-1) ;

(2,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xn) ;
(2,-.3) node $x_m$;

(xn.east |-x1.north)+(.35,+.25) node(cor1);
(x1.west|-x1.south)+(-.35,-.5) node(cor2);
(cor1) rectangle  (cor2);

background foreground

cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,

(0,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x1);
(0,-.3) node $x_1$;

(0.4,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](x2) ;
(0.4,-.3)node $x_2$;

(1,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xi) ;
(1,-.3) node $x_i$;

(1.68,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xn-1) ;

(2,0) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](xn) ;
(2,-.3) node $x_m$;

(0,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y1);
(0,1.2) node $y_1$;

(.4,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y2);
(.4,1.2) node $y_2$;

(1,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](yi);
(1,1.2) node $y_i$;

(1.68,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](yn-1) ;

(2,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](yn);
(2,1.2) node $y_m$;

(2.6,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](yn1);

(2.9,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](yn2);

(3.2,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](yn2);

(3.5,1) node[circle,draw=black, fill=black, inner sep=0pt,minimum size=2.5pt](y2n1);
(3.8,1.2) node $y_n$;

(temdis11)  [above =.33cm of x1] ;
(dis2) [right =-.15cm of temdis11]$\Delta$;

[-,thin,black,dotted](x1)edge node[swap,near start]  (y1);

(y2n1.east |-y2n1.north)+(.5,+.45) node(cor1);
(x1.west|-x1.south)+(-.35,-.5) node(cor2);
(cor1) rectangle  (cor2);

The ground truth is $X=\{x_1,\ldots,x_m\}$. In Figure <ref>, an algorithm estimates no target, i.e. $Z=\emptyset$ while in Figure <ref> an algorithm estimates $n$ targets ($n>m$), i.e $Y=\{y_1,y_2,\ldots,y_{n}\}$ where $\Delta<c$.

Analysis of the three metrics
for Figure <ref>

Metric   Fig.<ref> ($Z$)   Fig. <ref> ($Y$)    Which one
is closer to $X$?
OSPA    $c$     $\frac{m\Delta+(n-m)c}{n}$    $Y$
GOSPA   $c\frac{m}{\alpha}$    $m\Delta+\frac{(n-m)c}{\alpha}$
$Y$ if $\boldsymbol{\Delta< c\frac{2m-n}{m\alpha}}$

COSPA   $\dot{c}\left(2-\frac{1}{n}\right)$     $\frac{m\Delta+(n-m)\dot{c}}{n}$     $Y$
Intuition         $Y$

§ EXPERIMENT

We demonstrate the proposed metric by evaluating an multi-target tracking (MTT) algorithm together with OSPA and GOSPA metrics. We use the data and one of the result produced from MTT algorithm in [5] $38$ targets move from top right or middle of the surveillance
area to bottom left, and middle of the surveillance area to top
right. Each target
survives with probability $99\%$ and is detected with $80\%$. The measurements are added noize with  zero mean Gaussian process. The detected measurements are immersed in clutter modeled as a Poisson RFS with the average number of clutter returns per unit volume is $50$.
In this example, we use the OSPA metric [2], the GOSPA metric [3] and the proposed COSPA metric to compute the distances between the truth tracks and the estimated tracks (produced by the MTT algorithm).

Ground truth targets are immersed with their measurements and clutter.

For computation of the three metrics, we chose the cut-off parameter for the three metrics $c=80$, $\alpha=2$, $\xi=1$, $p=1$ and $\dot{c}=c+1$.
In this scenario,  Figure <ref> shows that at any time $t$ ($t=1,\ldots,50$) the GOSPA metric has much bigger error compared to OSPA and COSPA because GOSPA is a sum of all spatial distances and cardinality error at that time while OSPA and COSPA are the average of this sum. While Figure <ref> shows that most of the time OSPA and COSPA have the same errors except at time where the cut-off is applied or one of the two sets but not both is empty. It is because the cut-off distance $c$ and the penalty to each cardinality error are different in COSPA but the same in OSPA. Figure <ref> shows that most of the time all objects in the smaller sets associate correctly with objects in a bigger set except at time $t=21,38,39$ (two objects are correctly associated if their distance is smaller than the cut-off $c$). Indeed, OSPA gives value as a cut-off $c$ to the distance between the two sets at time $t=1$ and $t=50$ at which no track is detected. Furthermore, the distances and cardinality error at time steps $t=21,38$ and $t=39$ are the same in both COSPA and OSPA but OSPA local error is the sum of COSPA local error and COSPA outline error. It is because the distance between two vectors not to be smaller than cut-off parameter $c$ is considered as the wrongly associated tracks in COSPA. It happens because in COSPA the cut-off parameter $c$ is smaller than the penalty for each cardinality error $\dot{c}$.

Error versus time calculated of OSPA, GOSPA and COSPA metric.

Error versus time calculated under OSPA and COSPA metric.

Target numbers are plotted versus time

Figure <ref> shows that the GOSPA cardinality error is proportional to the difference in number between bigger set and the smaller set.
Indeed, the graph of cardinality error is $c/\alpha=80/2$ times the graph of the cardinality number. Similarly, Figure <ref> shows that the GOSPA localization error is proportional to the number of spatial distance across two sets.

Difference in number between two sets and Cardinality Error versus time calculated under GOSPA and COSPA metrics.

Difference in number between two sets and Cardinality Error versus time calculated under GOSPA and COSPA metrics.

§ CONCLUSION

This paper has discussed limitations of the  OSPA metric via some simple scenarios. Furthermore, some certain limitations of the GOSPA metric were discussed. GOSPA tries to overcome the limitations of the OSPA metric by removing the normalization of the OSPA metric and setting the cardinality penalty different from the cut-off parameter by multiplying the cut-off parameter with a positive number $\frac{1}{\alpha}$  that is larger than or equal to $\frac{1}{2}$ ($0<\alpha\le 2$). This alteration of the OSPA metric results in a greater penalty for a distance $\Delta$ between two vectors

than the penalty for each cardinality penalty if $c>\Delta>c/\alpha$. Hence the GOSPA metric will normally favor the empty set over the non-empty set when $\alpha>1$ and these two sets are compared with another non-empty set of the same size as the non-empty set and the distances between pairs of vectors across the two non-empty sets are equal to or larger than the cardinality penalty; or the localization error can be penalized more than the cardinality penalty if a distance $\Delta$ between two vectors is larger than $c/\alpha$ ($\Delta<c$).
Furthermore, the lack of normalization in the GOSPA metric makes it unreliable for measuring the distance between two non-empty sets.

The proposed COSPA metric was developed to overcome the shortcomings mentioned by [2] and also provide a practical assessment of the MTF  or multiple target tracking algorithms at a particular time in terms of missing targets, false targets, incorrectly associated targets and correctly associated target. Furthermore, the identities of missing targets, false targets and pairs of associated targets are provided in the process of calculating this metric (available as Matlab code in [6]). The COSPA metric retains the advantages of the OSPA metric, unlike the GOSPA metric. Thorough analysis of the COSPA metric reveals no major weaknesses,

noting that the penalty for each extra element in two sets is always larger than or equal to the cut-off distance between two vectors. The choice of the cut-off is problem dependent. We analysed the proposed metric together with other two metrics with some simple scenarios which shows the consistency and improvement compared the other two metric OSPA and GOSPA. We also use the set of tracks resulting from a multiple target tracking algorithm to evaluate the proposed metric compared with other two metrics.

§ ACKNOWLEDGEMENT

The author would like to thank Professor Henk Blom for the discussions, suggestions and comments.

[1]
J. Hoffman and R. Mahler, “Multitarget miss distance via optimal
assignment,” Systems, Man and Cybernetics, Part A: Systems and
Humans, IEEE Transactions on, vol. 34, no. 3, pp. 327 – 336, May 2004.

[2]
D. Schuhmacher, B.-T. Vo, and B.-N. Vo, “A Consistent Metric for Performance
Evaluation of Multi-Object Filters,” IEEE Transactions on Signal
Processing, vol. 56, no. 8, pp. 3447–3457, Aug. 2008.

[3]
A. S. Rahmathullah, Á. F. García-Fernández, and L. Svensson,
“Generalized optimal sub-pattern assignment metric,” in 2017 20th
International Conference on Information Fusion (Fusion), July 2017, pp.

[4]
T. Vu, “A complete optimal subpattern assignment (cospa) metric,” in
2020 IEEE 23rd International Conference on Information Fusion
(FUSION), 2020, pp. 1–8.

[5]
T. Vu, B.-N. Vo, and R. Evans, “A particle marginal metropolis-hastings
multi-target tracker,” Signal Processing, IEEE Transactions on,
vol. 62, no. 15, pp. 3953–3964, Aug 2014.

[6]
T. Vu. (2020) Matlab code for cospa metric. [Online]. Available: