# Three-dimensional simulations of X-ray cavities inflated by radio galaxies.

Michael D. Smith 1 & Justin Donohoe 1
1Centre for Astrophysics & Planetary Science, The University of Kent,
Canterbury, Kent CT2 7NH, U.K E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted ….. Received ….. ; in original form …..)

###### Abstract

Vast cavities in the intergalactic medium are excavated by radio galaxies. The
cavities appear as such in X-ray images because the external medium has been
swept up, leaving a hot but low density bubble surrounding the radio lobes. We
explore here the predicted thermal X-ray emission from a large set of high-
resolution three dimensional simulations of radio galaxies driven by
supersonic jets. We assume adiabatic non-relativistic hydrodynamics with
injected straight and precessing jets of supersonic gas emitted from nozzles.
Images of X-ray Bremsstrahlung emission tend to generate oval cavities in the
soft keV bands and leading arcuate structures in hard X-rays. However, the
cavity shape is sensitive to the jet-ambient density contrast, varying from
concave-shaped at $\eta=0.1$ to convex for $\eta=0.0001$ where $\eta$ is the
jet/ambient density ratio. We find lateral ribs in the soft X-rays in certain
cases and propose this as an explanation for those detected in the vicinity of
Cygnus A. In bi-lobed or X-shaped sources and in curved or deflected jets, the
strongest X-ray emission is not associated with the hotspot but with the relic
lobe or deflection location. This is because the hot high-pressure and dense
high-compression regions do not coincide. Directed toward the observer, the
cavity becomes a deep round hole surrounded by circular ripples. With short
radio-mode outbursts with a duty cycle of 10% , the intracluster medium
simmers with low Mach number shocks widely dissipating the jet energy in
between active jet episodes.

###### keywords:

hydrodynamics – galaxies: active – galaxies: jets – radio continuum: galaxies

††pagerange: Three-dimensional simulations of X-ray cavities inflated by radio
galaxies.–References††pubyear: 2019

##  1 Introduction

It is well established that powerful Active Galactic Nuclei are capable of
inflating bubbles in the intergalactic and intracluster media. This includes
vast cavities created when gas is driven out by powerful radio galaxies. The
cavities are observed as depressions in the surface brightness of X-ray
emission which encompass the radio lobes. Such cavities are found in diverse
cluster environments and surround a variety of radio galaxy types (Carilli et
al., 1994; Giacintucci et al., 2011; Pandge et al., 2012; Vagshette et al.,
2017; Pandge et al., 2019); with evidence for low Mach number bow shocks being
driven into, and capable of supporting, the environment (Wilson et al., 2006;
Stawarz et al., 2014). This phenomenon takes on extra significance as we probe
the regulatory role it has on galaxy formation through feedback from Active
Galactic Nuclei into the cluster gas (Fabian, 2012; Hardcastle et al., 2019).

With the advent of new X-ray technology, we can now expect detailed
measurements of the structure of these cavities at multiple frequencies. Our
main purpose here is to predict how this X-ray structure would appear if the
radio galaxy is driven by supersonic hydrodynamic jets. We aim to account for
the Mach number, precession and orientation in order to provide a rich bank of
simulated images. The physical properties and radio morphologies have already
been published for this set of three-dimensional simulations in Paper 1
(Donohoe & Smith, 2016) and Paper 2 (Smith & Donohoe, 2019). We thereby wish
to probe the implications for the intracluster medium. Specifically, these
cavities may provide the link between the regulation of star formation in
central massive galaxies and cluster cooling flows. In this respect, we
already confirmed in Paper 1 that a high fraction of the jet energy is
transferred and dissipated as thermal energy into the intracluster gas.

An X-ray cavity was first associated with the radio source 3C 84 in the
Perseus elliptical galaxy NGC 1275 by Boehringer et al. (1993). The Perseus
intracluster medium was found to be displaced by the radio lobes although the
minimum lobe pressure falls below that of the thermal gas. The brightest
regions are located just outside the radio lobes. In addition, the jets appear
to change angle from the south-east towards the south with the lobes extending
towards the west. While this may indicate jet deflection, it could also be
directly related to jet precession. Note that the X-ray ROSAT image is from
the soft range between 0.1 and 2.4 keV.

Another X-ray cavity was then associated with the radio galaxy Cygnus A by
Carilli et al. (1994). Such X-ray cavities are likely to endure, remaining
intact for over 100 Myr, well after the high energy relativistic electrons
have cooled within the lobes. Superimposed X-ray structure can also take the
form of ripples, edges, bow shocks and cold fronts (e.g. Bîrzan et al., 2004;
Croston et al., 2011; Shin et al., 2016; Sonkamble et al., 2015).

Numerical simulations of jet-driven radio galaxies have shown that detectable
X-ray cavities should surround the radio lobes. A three-dimensional simulation
by Clarke et al. (1997) of a straight light jet was used to model Cygnus A.
Background subtracted integrations for various energies were presented
assuming jet motion in the sky plane. A number of two-dimensional computations
by Hardcastle & Krause (2013) of very high Mach number jets confirmed the
cavity-shell structures for jets in the sky plane. The X-ray cavities
associated with a pair of full magnetohydrodynamic simulations were presented
by (Mendygral et al., 2011). Although full two-sided radio galaxies with
evolving jet powers were considered in this work, no asymmetries resulted due
to the imposed axisymmetric initial conditions. It was found that cavity age
and, therefore, cavity power are sensitive to the accuracy of the estimated
inclination angle of the jet axis. Cavity age and power estimates within a
factor of two of the actual values are possible given an accurate inclination
angle. However, even more crucial to the cavity shape is the range of ambient
and jet properties. For example, the jet-ambient density contrast is able to
generate top-wide, centre-wide and bottom-wide X-ray structures as found in
the two-dimensional simulations of Guo (2015) although only the total
Bremsstrahlung cooling rate was modelled.

X-ray surface brightness maps were also generated by Turner & Shabala (2020).
However, these were based on semi-analytical calculations for mock FR II radio
galaxies, enabled after several assumptions including self-similarity and
axial symmetry. Although a crude model, this work reminds us that the X-ray
emission from distant radio galaxies should become dominated by inverse
Compton emission off background photons at high redshifts.

We are now able to perform high resolution three dimensional simulations with
sufficient efficiency to systematically investigate the influence of jet
dynamics and geometry. This approach will assume a uniform intracluster medium
within which the driving galactic nucleus is stationary. We will here study
the large-scale X-ray cavities and shocks/edges. We mainly analyse the
distribution of the hot gas through thermal Bremsstrahlung radiation but
neglect the internal lobe X-rays generated from synchrotron, synchrotron self-
Compton.and inverse Compton processes. Free-free or Bremsstrahlung radiation
generally has a low emissivity and so is associated with the large quantities
of hot gas in cooling flows (Jaffe & Perola, 1974).

Synchrotron emission at X-ray energies is detected from within radio jets and
hot spots where strong particle acceleration occurs across shock fronts (e.g.
Atoyan & Dermer, 2004).

Inverse Compton emission is also deduced from within extended radio jets. In
this case, the supply of photons to be upscattered can be produced internally
via synchrotron emission, hence the term synchrotron self-Compton.
Alternatively, the photons can be from the cosmic microwave background, galaxy
starlight or radiation from the active nucleus (e.g. Wilson & Yang, 2002;
Kataoka & Stawarz, 2005).

We will here also investigate the ability of the radio galaxy to disrupt a
cluster cooling flow. To achieve this, we determine the speed of advance of
the shock wave that propagates into the cluster environment. Besides the Mach
number of the advancing bow shock, it is more relevant to provide the Mach
number associated with the temperature jump since that is the observable
quantity. For example, Croston et al. (2011) report evidence for a large
spheroidal shock in the cluster surrounding the powerful radio galaxy 3C 444.
A temperature jump of approximately1.7 corresponds to a Mach number of 1.7 for
a steady hydrodynamic shock. However, Croston et al. (2011) remark that this
is likely to be an underestimate since an average over a large region is made.

Figure 1:  X-ray images generated from the data cube for Model ZDA without the
uniform background subtraction. The 1 keV emission is displayed for the jet
axis in the sky plane (top panel) and 25∘ out of the plane (lower panel) for
jet-ambient density ratio $\eta=0.1$, Mach 6 straight jet with 1∘ precession
angle. This is the free-free X-ray surface brightness; all other images in
this paper are sky/background subtracted.  Figure 2:  Background subtracted
images for $\eta=0.1$, Mach 6 jet with 1∘ precession angle for Model ZDA. The
top three panels display the free-free X-ray surface brightness at 1 keV, 10
keV and 40 keV while the lower panel displays a radio image based on a simple
synchrotron model. The jet axis is in the plane of the sky. Figure 3:
Background subtracted images for $\eta=0.01$, Mach 6 jet with 1∘ precession
angle (for Model ZDB). The top three panels display the free-free X-ray
surface brightness at 1 keV, 10 keV and 40 keV while the lower panel displays
a radio image based on a simple synchrotron model. The jet axis is in the
plane of the sky . Figure 4:  Background subtracted images for $\eta=0.001$,
Mach 6 jet with 1∘ precession angle (for Model ZDC). The top three panels
display the free-free X-ray surface brightness at 1 keV, 10 keV and 40 keV
while the lower panel displays a radio image based on a simple synchrotron
model. The jet axis is in the plane of the sky. Figure 5:  Background
subtracted images for a very light $\eta=0.0001$, Mach 6 jet with 1∘
precession angle(for Model ZDM). The top three panels display the free-free
X-ray surface brightness at 1 keV, 10 keV and 40 keV while the lower panel
displays a radio image based on a simple synchrotron model. The jet axis is in
the plane of the sky . Figure 6:  The Mach 6 straight jet with the jet axis at
an angle of 25∘ to the plane of the sky. These are background subtracted
images produced by a jet with a jet-ambient density ratio of 0.01 and a slow
precession with half-angle of 1∘ (for Model ZDB). The top three panels display
the free-free X-ray surface brightness at 1 keV, 10 keV and 40 keV while the
lower panel displays a radio image based on a simple synchrotron model.
Figure 7:  The Mach 6 straight jet with the jet axis at an angle of 75∘ to the
plane of the sky. These are background subtracted images produced by a jet
with a jet-ambient density ratio of 0.01 and a slow precession with half-angle
of 1∘ (for Model ZDB). The top three panels display the free-free X-ray
surface brightness at 1 keV, 10 keV and 40 keV while the lower panel displays
a radio image based on a simple synchrotron model. Figure 8:  The
$\eta=0.0001$ Mach 6 straight jet with the jet axis at an angle of 75∘ to the
plane of the sky. These are background subtracted images produced by a jet
with a jet-ambient density ratio of 0.0001 and a slow precession with half-
angle of 1∘ (for Model ZDM). The top three panels display the free-free X-ray
surface brightness at 1 keV, 10 keV and 40 keV while the lower panel displays
a radio image based on a simple synchrotron model. Figure 9:  A Mach 2
straight jet (fModel M2), viewing angle 25∘ and $\eta=0.1$ with 1∘ precession.
Background subtracted images. The top three panels display the free-free X-ray
surface brightness at 1 keV, 10 keV and 40 keV while the lower panel displays
a radio image based on a simple synchrotron model. The jet axis is at an angle
of 25∘ to the plane of the sky. The length scales are expressed in
kiloparsecs. Figure 10:  A Mach 24 straight jet, viewing angle 0∘ and
$\eta=0.1$ (for Model IC). Background subtracted images. The top three panels
display the free-free X-ray surface brightness at 1 keV, 10 keV and 40 keV
while the lower panel displays a radio image based on a simple synchrotron
model. The jet axis is in the plane of the sky.  Figure 11:  A slow precessing
jet with 20∘ precession angle (for Model ZDG). image is from the projected X-Y
plane. Images produced by a Mach 6 jet with a jet-ambient density ratio of
0.1. The top three panels display the free-free X-ray surface brightness at 1
keV, 10 keV and 40 keV while the lower panel displays a radio image based on a
simple synchrotron model. The jet axis is at an angle of 25∘ to the plane of
the sky. Figure 12:  A slow precessing jet with 20∘ precession angle. image is
from the projected X-Z plane. Images produced by a Mach 6 jet with a jet-
ambient density ratio of 0.1 (for Model ZDG). The top three panels display the
free-free X-ray surface brightness at 1 keV, 10 keV and 40 keV while the lower
panel displays a radio image based on a simple synchrotron model. The jet axis
is at an angle of 25∘ to the plane of the sky.

Table 1: Summary of jet properties as expressed for a domain scaled to 150 kpc. The run identification code is consistent with previous papers in this series with the label Z denoting a reflection boundary condition on the injection plane outside the nozzle and the D corresponding to larger 3003 domains. Figure | Code used | Jet Mach | Density | Prec. | Jet Sound Speed | Jet Power
---|---|---|---|---|---|---
number | in Paper 1 | number | ratio $\eta$ | angle | cm s-1 | erg s-1
1/2 | ZDA | 6 | 0.1 | 1∘ | 3.24$\times$108 | 6.42$\times$1044
3 | ZDB | 6 | 0.01 | 1∘ | 1.02$\times$109 | 2.03$\times$1045
4 | ZDC | 6 | 0.001 | 1∘ | 3.24$\times$109 | 6.42$\times$1045
5 | ZDM | 6 | 0.0001 | 1∘ | 1.02$\times$1010 | 2.03$\times$1046
6 | ZDB | 6 | 0.01 | 1∘ | 3.24$\times$108 | 6.42$\times$1044
7 | ZDB | 6 | 0.01 | 1∘ | 3.24$\times$108 | 6.42$\times$1044
8 | ZDM | 6 | 0.0001 | 1∘ | 1.02$\times$1010 | 2.03$\times$1046
9 | M2 | 2 | 0.1 | 1∘ | 3.24$\times$108 | 2.38$\times$1043
10 | IC | 24 | 0.1 | 1∘ | 1.02$\times$109 | 1.30$\times$1047
11/12 | ZDG | 6 | 0.1 | 20∘ | 3.24$\times$108 | 6.42$\times$1044
14-xz | ZDG4 | 6 | 0.1 | 20∘ fast | 3.24$\times$108 | 6.42$\times$1044

##  2 Method

###  2.1 Simulation and display parameters

The simulations analysed here have been introduced in Paper 1 (Donohoe &
Smith, 2016) for their physical attributes and in Paper 2 (Smith & Donohoe,
2019). for their radio morphology and classification. The set of defining
parameters and corresponding figures are provided in Table 1.

The analysis was performed with algorithms incorporated into IDL software.
Before inspecting the figures, some common properties should be noted in
advance (in order to maximise the panel sizes and reduce the caption lengths).
The length scales are expressed in kiloparsecs although the pixel unit is 0.5
kiloparsecs which corresponds to the zone sizes of the fixed-grid Eulerian
Pluto code. Luminosities are in solar units with the total Bremsstrahlung
cooling (total), full emission from the specific band per keV (full) , and
subtracted surface brightness (sub) where the background has been subtracted
(not subtracted for the first figure). See Table 2 for a summary.

The simulations were performed with PLUTO, a grid-based code incorporating
modern Godunov-type shock-capturing schemes (Mignone et al., 2007). After
comparing the results of numerous options, we chose a fast linear
interpolation time-stepping (denoted HLLC) scheme as developed by Harten, Lax
and Van Leer and detailed by Toro et al. (1994). The simulations were mainly
performed on grids of 3003 unless otherwise stated. Hence the voxel and pixel
sizes shown in the figures corresponds to 0.5 kpc.

We have performed a resolution study with grids of 753, 1503, 2253 and 3003.
We presented results in Paper 1 that showed there was a gradual convergence
with 1503 sufficient for most purposes and 2253 oroviding reliable hot spot
properties.

###  2.2 Parameters for the ambient medium

The temperature of the environment is taken as $T=4.56\times 10^{7}$ K,
corresponding to 4 keV, with a hydrogen nuclei density of $n=1.21\times
10^{-3}$ cm-3 to simulate one half of a Radio Galaxy within a cube which
stretches 150 kpc. This is adjustable to other densities and length scales
quite easily since the flows are adiabatic with a specific heat ratio of 5/3.
It would also be possible to scale the derived emission provided the ratio of
the X-ray energy to the temperature were held fixed.

A uniform ambient medium is taken as opposed to the density profiles commonly
assumed to model galactic atmospheres such as by Clarke et al. (1997). We
model the intracluster medium in this instance but realise that cluster
cooling flows and galactic motions would need us to introduce detailed
structure.

###  2.3 Parameters for the jet

We configured the code for non-relativistic hydrodynamics. We do this in order
to explore ranges in the dynamics in full three-dimensions with a reasonable
resolution. Full details were published in Paper 1 for the physical set up and
in Paper 2 for the radio emission.

The jet is initially in pressure balance with the ambient medium but the sound
speed in the jet is considerably higher given the low density jet which is
expressed by the density ratio $\eta$. It is important to allow the jet
direction to vary, not only to break symmetry but also because there is
considerable evidence for this (Krause et al., 2019).Therefore, the jet
direction is precessed through sinusoidal velocity changes in the y and z
directions as detailed in Donohoe & Smith (2016). The standard slow precession
rate has an angular frequency of $\pi$/2 in units of $c_{s}/r_{j}$ where
$c_{s}$ is the ambient sound speed and $r_{j}$ is the initial jet radius. This
corresponds to a precession period of 19.1 Myr with the present units.

###  2.4 Parameters for the emission

To calculate the total Bremsstrahlung emission from the simulations, we assume
full ionisation, with 10% helium by number,. From the theoretical foundations
(Tucker, 1975; Rybicki & Lightman, 1979), there are a number of approximations
fo the formula summarised by Yarza et al. (2020) from which we employ

$L_{\text{brem}}=1.42\times 10^{-27}\cdot n_{e}^{2}\cdot
T^{1/2}\text{erg}~{}\text{s}^{-1}~{}\text{cm}^{-3},$ (1)

where $T$ is the temperature in Kelvin and $n_{e}$ is the electron density
measured in cm-3.

We present emission maps at soft kTo=1 keV ($2.42\times 10^{17}$ Hz),
intermediate kTo = 10 keV and hard kTo = 40 keV X-rays. The formula employed
is

$L_{\text{o}}=6.8\times 10^{-38}\cdot n_{e}^{2}\cdot T^{-1/2}\cdot
exp(-T_{o}/T)\text{erg}~{}\text{s}^{-1}~{}\text{cm}^{-3}~{}\text{Hz}^{-1}.$
(2)

Note that we will consider the emission at three specified energies but
multiply up to band passes of width 1 keV and provide values per square
kiloparsec which generates meaningful luminosities. Scales are provided on
bars but are adjusted by appropriate factors as indicated in their titles for
display purposes.

##  3 Results

###  3.1 Density dependence

The density of the jet relative to the ambient medium is shown to have a major
influence on the X-ray distribution. We define $\eta$ as the jet-ambient
density ratio and consider light jets with $\eta$ between 0.1 and 0.0001.

With $\eta=0.1$, a substantial depression in the soft X-rays is revealed as
displayed in the top panel of Fig. 1. This particular figure is not background
subtracted and the data set is confined to within a 150 kpc $\times$ 150 kpc
region filled with 3002 pixels. The display has been chosen so that when the
3003 data cube is rotated (see below) , the image will remain in the window
with the origin of the jet remaining at a fixed point.

Probably more instructive are the equivalent background subtracted images
displayed in Fig. 2. The central cavity is enveloped by a thick elliptical
shell. As expected, the highest surface brightness coincides with the leading
edge where the heating and compression are strongest. However, the total
emission is dominated by the extended flanks which consist of locally shocked
gas and adiabatically-expanding hot gas flowing within.

This shell is much dimmer at intermediate X-ray energies shown in the second
panel while this is then reduced to an arc at hard X-rays in the third panel.
It is clear that the hard X-ray energy falls well into the exponential part of
the emission function for all except the very hot leading part of the bow. In
addition, the elliptical shell is distorted with stronger emission to the
’south’ in all bands.

A different structure is apparent for $\eta=0.01$ shown in Fig. 3. Firstly,
the shell is better described as conical when the jet axis lies in the sky
plane. Ribs of emission are now very prominent, running across the shell
perpendicular to the jet axis. Four ribs can be seen at the stage shown. The
separation between is approximately 60 kpc. The nature of the ribs will be
elucidated below through viewing the structure at another angle.

At this lower density, the leading shock is stronger and wider. The jet has
penetrated across the grid in a shorter time of 20 Myr, displaying a
streamlined structure rather than the blunt appearance of the $\eta=0.1$ case
above.

With the jet density further reduced, as shown in Fig. 4 and 5, the ribs
disappear. Instead, there is a prominent inner cavity in soft X-rays, which is
roughly oblate. This is enshrouded by bright shoulders. The visible knotty
radio jet now penetrates through the inner cavity although the hot cocoon
immediately surrounding it does not seem very distinct. At the lowest density
ratio, the entire disturbed region contains two cavities: the inner oblate at
all energies and an outer circular bubble at intermediate and hard energies.

###  3.2 Orientation of the jet axis

The integrated line of sight emission with the jet axis at various angles is
presented with the panel widths having been chosen so that the entire
projected data cube is visible and the jet origin remains at a fixed location,
the origin (0, 0).

The X-ray structure is altered systematically as the orientation out of the
plane of the sky is increased. This effect is illustrated in Figs. 6 at 25∘and
Fig. 7 at 75∘. to the plane of the sky.
At soft X-ray energies, the cone turns into a spindle and the ribs become
round ripples. At progressively high energies, the deep cavity becomes smaller
and shifts towards the front. Note also that secondary radio hotspots are
present at these angles due to projection effects of filamentary structures.

A remarkable short knotty radio jet is revealed in the lower panel of Fig. 8
for the lowest density jet at a low inclination to the line of sight. This jet
is surrounded by an extended radio halo. Moreover, a similar sized thick X-ray
front exists at the edge of the halo. Therefore, this simulation may have
relevance to a number of quasar jets with halo properties.

###  3.3 Jet Mach number

At low Mach numbers, the radio galaxy displays considerable turbulent
structure as shown in the lower panel of Fig. 9. Although on these scales the
jet has not completely disintegrated, there is no advancing hotspot. The soft
X-ray cavity precisely follows the radio lobe. Ahead of the lobe is a very
wide and thick shell of enhanced thermal emission. Hence, low Mach number jets
can energise a large volume of the intracluster medium.

In hard X-rays, the structure is reversed with only the very hottest regions
contributing. This includes a sheath of hot gas surrounding the jet out to 60
kiloparsec, detected here as two thin parallel filaments. This is followed by
a turbulent zone coinciding with the jet disruption. Hence, Bremsstrahlung
emission can arise internal to the radio lobes at low Mach numbers, partly due
to the stronger dissipation within the radio galaxy and weaker transmission to
the intracluster medium (see Column 11 of Table 2.

At high Mach numbers, a contrasting pattern emerges. As shown in Fig. 10 for a
Mach number of 24 and $\eta=0.001$, a thimble-like cavity is commensurate with
the streamlining due to the high Mach number. This is consistent with the
short grid crossing time of 7 Myr. Note that the injected jet is pressure-
matched to the environment. That implies that a low-density jet will possess a
high sound speed. This compensates in terms of the momentum flow for a given
Mach number. Here, however, the advanced speed should be proportional to the
Mach number, as indeed found.

###  3.4 Precession angle

The Mach numbers for shocks propagating ahead of radio galaxies have been
estimated through their X-ray properties. The shocks tend to be very weak with
Mach numbers below 3. Moreover, the radio lobes are generally not greatly
streamlined. Together, these facts suggest that the speed of propagation into
the intracluster medium is only mildly supersonic despite the very high Mach
numbers of the driving jets. One way to achieve this in the present context is
through a slow precession of the jet axis so that the jet momentum is
distributed over a wide surface area.

Precession can manifest itself as multiple radio sub-lobes as shown in Fig.
11, reminiscent of X-shaped or double-double radio galaxies. In X-rays, this
simulation has a unique signature: the high surface brightness is associated
with the relic radio component rather than with the present hotspot. This
holds for the soft and intermediate X-ray energies but the newly advancing
lobe is more apparent at hard X-rays. Hence, the Bremsstrahlung emission can
remain dominated by the large amount of swept-up and compressed material from
the previously compressed environment while the jet rotates and progresses
through an earlier ejection of lobe material.

Rotated about the jet axis by 90∘, yields a contrasting radio morphology as
displayed in the lower panel of Fig. 12. The jet begins quite straight snd
becomes slightly curved before deflection at a secondary hotspot. In this
case, the peak soft X-ray emission does not coincide with either hotspot but
appears where the overspill from the primary hotspot has created a new lobe to
the top of the image. A further thick rim is present on the opposite side.

On comparison of the two images, it is apparent that there are two circular
cavities. The smaller inner cavity corresponds to the region associated with
the relic lobe while the large cavity is associated with the overall evacuated
lobe.

###  3.5 Luminosity & energetics

The total Bremsstrahlung luminosity is indicated directly on the third panel
of the figures. In the absence of the radio galaxy, the luminosity from the
150 kpc cube is 3.61 $\times$ 108 L⊙. It should be noted that scaling to a
sphere of radius 150 kpc with a hydrogen nuclei density of 3 $\times$ 10-3
cm-3 would yield an X-ray luminosity of 1.35 $\times$ 1010 L⊙. Cooling times
remain longer than the age of the Universe.

The introduction of the jet increases the total luminosity in all cases
despite the cavity formation. The cavity forms through the compression of
material. and the shell or rim is generally responsible for enhancing the
Bremsstrahlung emissivity. The one exception is the Mach 2 jet in which
significant amounts of ambient material have been expelled from the grid.

The total luminosity increases as the jet density is reduced but not by a
large factor, approximately doubling over the entire density range of a 1,000.
This reflects the fact that a lower density is matched by a higher jet sound
speed, which maintains the same momentum flow rate and approximately the same
advance speed of the radio hotspot and interface.

Also shown in each panel are the full Bremsstrahlung fluxes in solar units per
keV. Underneath this, the background subtracted fluxes are given. This resets
the zero point on the colour bar, allowing the cavities to be emphasized as
negative values. In this respect, the soft X-ray cavity does deepen as the jet
density decreases while, on the other hand, the shell surface brightness
increases. Overall, upon background subtraction, the net integrated flux
combining enhanced and reduced regions is relatively small. It is also clear
that cavities are prevalent in soft X-rays but not in hard X-rays where the
background is almost negligible and the emission is mainly generated from the
leading edge of the ambient bow shock.

The size and depth of the cavities are quantified and presented in Column 9 of
Table 2. The radio lobe area is also listed in Column 8. The detectable lobe
is here given a threshold of 1% of the peak value. Hence the radio lobe and
the cavity area both generally increase as the jet-ambient density ratio
decreases. It can also be noted that the radio area and cavity area can
greatly depend on the viewing angle. The depth of the X-ray cavity is measured
as the total surface brightness of the cavity within the area, defined as 0.3
below the background in the figures. It is evident that cavities are deeper
when the jet axis is closer to the line of sight. However, several factors
contribute. For example, the high Mach number straight jet has blown a deep
cavity although the streamlined shape means that it has not vacated a large
volume. On the other hand, the low Mach number jet has also cleared a deep
cavity coincident with the radio lobe.

Table 2: Properties of X-ray images of simulated radio galaxies. Fig. | Mach | Density | Prec. | Angle to | Interface | Bow | Radio | Cavity | Cavity depth | Energy | X-ray | X-ray
---|---|---|---|---|---|---|---|---|---|---|---|---
| number | ratio $\eta$ | angle | sky plane | Mach No. | Mach No. | area kpc2 | area kpc2 | L${}_{X}/10^{6}$L⊙ | transfer | Temp.107K | Mach No.
1/2 | 6 | 0.1 | 1∘ | 0∘ | 4.35 | 4.48 | 142 | 4,592 | -1.99 | 0.71 | 5.23 | 1.12
3 | 6 | 0.01 | 1∘ | 0∘ | 5.22 | 5.42 | 214 | 4,634 | -2.11 | 0.86 | 6.00 | 1.30
4 | 6 | 0.001 | 1∘ | 0∘ | 6.10 | 6.29 | 464 | 4,187 | -2.91 | 0.92 | 7.48 | 1.60
5 | 6 | 0.0001 | 1∘ | 0∘ | 5.89 | 6.18 | 10,292 | 6,940 | -15.78 | 0.94 | 10.25 | 2.09
6 | 6 | 0.01 | 1∘ | 25∘ | 5.22 | 5.42 | 380 | 4,533 | -2.42 | 0.87 | 6.00 | 1.30
7 | 6 | 0.01 | 1∘ | 75∘ | 5.22 | 5.42 | 543 | 4,575 | -5.00 | 0.87 | 6.00 | 1.30
8 | 6 | 0.0001 | 1∘ | 75∘ | 5.89 | 6.18 | 15,785 | 11,745 | -8.03 | 0.94 | 10.25 | 2.09
9 | 2 | 0.1 | 1∘ | 25∘ | 1.07 | 1.19 | 5,740 | 19,348 | -4.34 | 0.59 | 4.71 | 1.01
10 | 24 | 0.1 | 1∘ | 25∘ | 35.50 | 37.25 | 2,765 | 8,566 | -6.83 | 0.86 | 6.40 | 1.38
11-xy | 6 | 0.1 | 20∘ | 25∘ | 0.80 | 3.30 | 3,024 | 3,769 | -2.02 | 0.79 | 5.23 | 1.13
12-xz | 6 | 0.1 | 20∘ | 25∘ | 0.80 | 3.30 | 1,960 | 3,989 | -1.99 | 0.79 | 5.23 | 1.13
14-xz | 6 | 0.1 | 20∘ fast | 50∘ | 1.44 | 2.37 | 7,217 | 6,532 | -5.79 | 0.82 | 5.50 | 1.19

Figure 13:  The Mach number of the bow shock in the cluster environment as
derived from the maximum jump in pressure taking steady shock jump conditions.
The value is derived along each y-zone in the z mid-plane.

###  3.6 Intracluster shock Mach number

There are several useful Mach numbers that can be derived from the
simulations. The first is from the physical properties. Here, we present the
values for the time-averaged Mach numbers along the jet axis which are easily
derived for both the interface and the shock front running ahead of it. These
are listed in Columns 6 & 7 of Table 2. Note that these are generally high
values but appropriate for estimating source ages. The bow cap is typically
just 4–5% ahead of the interface apart from the Mach 2 jet which has a 1%
stand-off and the anomalous precessing jet which appears as two distinct off-
axis lobes encompassed by a surrounding large-scale bow.

The Mach number attributed to the shock driven into the cluster medium can be
derived from spectral measurements in the X-ray bands. One measure of the Mach
number is calculated from the ratio of X-ray emission in the 1 and 10 keV
bands. This generates an X-ray temperature which can be compared to that
expected from the jump conditions for a hydrodynamic shock. These values,
listed in the final column of Table ,2, are very low since the extended low
Mach number flanks of the bow dominate.

It can be concluded that the derived Mach number depends on the resolution.
This is demonstrated in Fig. 13 where shock hydrodynamic jump conditions have
been employed to convert the maximum change in temperature or pressure along
the x-axis into an equivalent Mach number for a steady shock. While a
temperature jump may be measured observationally, in numerical simulation with
an extreme density contrast, the pressure provides a better solution. This
figure demonstrates that any observation of the curved bow will be strongly
dependent on the spatial resolution. Therefore, extreme caution must be
applied when using derived Mach numbers to deduce radio galaxy energetics or
dynamics.

###  3.7 Evolution of X-ray flux density

Can we identify candidate X-ray cavities from the total emission in the
ambient medium? To determine this, we calculate the total Bremsstrahlung
emission from the data cubes in the soft, medium and hard X-ray bands.
Moreover,by presenting the time evolution, we gain a further potential
diagnostic of radio galaxy age.

Figure 15 shows that the total X-ray emission in each band increases despite
the reduction within the cavity. This applies to all bands and is a result of
the density squared dependence of the emissivity. Hence, when the radio lobes
sweep up the material, they create significant density variations out of the
ambient medium while the total mass remains unaltered. Therefore, the local
depression is more than compensated by the enveloping region of enhanced
emission.

This increase in total emission is mostly detectable in the harder bands. At
40 keV, there is minimal emission before the jet is initiated and the emission
strongly increases due to the direct impact. In this case, the emission is
sensitive to the rise in temperature across the impact region.

As the source ages, ambient material is ejected from the grid. This was found
in the simulations with very low jet densities, as shown by the turn over in
the emission curves. This is due to the large lateral expansion of the cavity.
On the other hand, although a wide radio source is created in the case of a
high jet precession angle, the X-ray emission is not significantly increased
in comparison to the straight jet with $\eta=0.1$. We thus conclude that for
the enhancement in integrated X-ray emission from the Bremsstrahlung process,
the jet density is the crucial factor.

##  4 Discussion

A major result is that the cavity morphology is a function of the jet-ambient
density ratio. Note that this treats the jet Mach number and jet-ambient
pressure ratio as the two other independent variables. We study here only the
pressure-matched case and, for a fixed Mach number, this implies a constant
momentum flow rate independent of how light the jet is. The pressure balance
condition might be critical to the jet stability and will be relaxed in a
follow-up study. Here, we find that for $\eta=0.1$, an ellipsoidal cavity
demarcated by a shocked shell appears as an elliptical X-ray cavity and bright
shell. As shown in Fig. 1. the shell is typically 20%–25% thick and is likely
to be asymmetrical, ultimately caused by a preferred direction of deflection
of the jet in the leading impact region.

These cavity results differ from those of Guo (2015). A comparison is
difficult since their simulations are two dimensional, involve viscosity and
consider subsonic jets that generate top-wide cavities. The main difference,
however, is that their simulations only considered cavities well after the jet
had been switched off, allowing a buoyant bubble to rise in a galactic
potential.

As $\eta$ is reduced in the simulation sequence, the shell shape becomes
progressively more conical with the high surface brightness regions located
further into the flanks. Most remarkable is the appearance of ribs which cut
across the cavity when the jet is in the plane of the sky. The ribs show up as
distinct brightenings along the shell and finally as ripples at viewing angles
close to the line of sight (fig. 7). This indicates that the structure is
caused by oscillatory motion in the impact zone which has been known to
experience periodic vortex shedding associated with periodic transitions
between double hotspots and single hotspots on the axis (Smith et al., 1985).
The shed vortex feeds back on to the approaching jet, squeezing it before the
gas is shocked at the hotspot. The shedding time scale is estimated to be the
hotspot radius divided by the ambient sound speed. This yields four complete
cycles for the $\eta=0.01$ simulation, which is consistent with the simulated
structure. This is confirmed through closer inspection of the radio image in
the lower panel of Fig. LABEL:ZDB-yes-023-00 which displays previously shed
vortices on both sides of the radio cocoon, corresponding to the location of
the X-ray ribs.

Enhanced X-ray structure within the Cygnus A thermal cocoon has been
identified by Smith et al. (2002) as belts. With Chandra data, the description
was refined to cross-cutting thick ribs which appear to be wrapped on the
inside of an ellipsoidal shell (Duffy et al., 2018). It is certainly
reasonable to associate the observed ribs with the rib cage which occurs
naturally in these simulations with the minimum of necessary ingredients.
Further examples would be expected as we gain high resolution data for other
classical doubles.

For the lightest jets, the shell develops wide shoulders on the flanks and a
deep ovoid cavity. The radio jet penetrates through the cavity and propagates
across the domain. These X-ray attributes are very similar to those associated
with 3C 444 (Vagshette et al., 2017) and 3C 320 (Vagshette et al., 2019) which
possess double inner cavities with an adjacent rim plus outer shocks.

Radio galaxies generated by jets within a wide precession cone would extend
the lateral reach of feedback into the intracluster medium. This is evident
from the low Mach numer of the shocks into the cluster environment for the 20∘
precession angle cases in Table 2. However, the fraction of jet energy
transferred into the surroundings remains at a substantial level. Therefore,
precession does not actually help to increase external support but may help
spread this support of the gas across an entire cooling flow (Riffel et al.,
2019) and so inhibit galaxy growth, two major issues recently summarised by
Hardcastle et al. (2019). Precession may prove to be a prominent dynamical
mechanism in distant radio galaxies as they form in proto-clusters. Therefore,
we concentrate below on the effects of precession and the influence this may
have on their apparent morphological classification.

The dynamical times of the simulated radio galaxies are significantly shorter
than the lifetime of a hot dense cooling cluster flow. This would indicate
that the present active phase of an AGN would not be directly responsible for
the onset of measurable feedback and increase of entropy in the cluster gas
(Belsole et al., 2007). Hence, it may be crucial to study how an environment
responds after a radio galaxy power supply has been switched off. Indeed,
cavities without enveloping shock waves are observed which suggest multiple
episodes (Choi et al., 2004) with shock decay in between.

An excellent template for an observed X-ray cavity and shell belongs to the
environment of the fat double FRİ 3C 310 (Kraft et al., 2012). it is a prime
candidate for a radio galaxy formed by a fast precessing jet which has been
able to draw out the radio hotspot into a ring in the southern lobe and a
broken round arc in the northern lobe. In this case, the energetics have been
clearly outlined by Kraft et al. (2012). The X-ray cavity and radio galaxy are
powered by an ongoing outburst releasing thermal energy into the gas within
the shocked inner 180 kpc of 1061 ergs and which emits $\sim$ 2 $\times$ 1043
erg s-1 while the outburst energy is $\sim$ 1061 ergs. The Mach number of 1.9
for the expansion is derived from a model for the shock density jump. While
these numbers are highly uncertain, this leads to an estimated time-averaged
outburst power of 2 $\times$ 1045 erg s-1 and suggests that the AGN could
easily support the cooling flow with a duty cycle of just 10%.

Figure 14:  A fast precessing Mach 6 jet with 20∘ precession angle (for Model
ZDG4). This image is from the projected X-Z plane. The jet-ambient density
ratio is 0.1. The top three panels display the free-free X-ray surface
brightness at 1 keV, 10 keV and 40 keV while the lower panel displays a radio
image based on a simple synchrotron model. The jet axis is at an angle of 50∘
to the plane of the sky.

Figure 15:  Evolution of the flux density. The total X-ray emission in each
band of spectral width 1 keV over the entire 150 kpc cube as a function of
time.

Here we employ our model for the radio emission of a fast precessing jet. The
jet has a 20∘ precession angle and a period of one quarter of the above
standard models (i.e. a period of $r_{j}/c_{s}$ where $r_{j}$ is the jet
radius and $c_{s}$ is the ambient sound speed. Figure 14 displays a fat radio
lobe with a ring-like structure in the lobe apparent for this angle to the
line of sight. The total X-ray power of the simulation of 2 $\times$ 1042 erg
s-1 and jet power of 6.4 $\times$ 1044 erg s-1 for this one half of the radio
galaxy within the 150 kpc cube compares favourably without any need for
parameter manipulation. Note the outburst age of 48 Myr and a derived time-
averaged Mach number of the leading edge of 1.44 for the interface and 2.37
for the disturbed ambient medium. These values are commensurate with those
derived from the jump conditions which yield a maximum Mach number of 2.15.
Note that this simulation was performed with extensions in the form of a
staggered grid beyond the uniform displayed portion.

X-ray cavities can take many forms and provide a means of investigating galaxy
and cluster formation at high redshift since cluster cooling flows are strong
X-ray sources. MACS J1447.4+0827 (z = 0.3755) provides a recent example of jet
evacuated cavities (Prasow-Émond et al., 2020). Unlike Cygnus A, in this case
two separate cavities on a scale of 100 kpc are detected. Moreover, the
cavities appear to open out away from the central galaxy. While this may be
simulated here, the strong emission from the soft X-ray band occurs on the
inner section of each side of the source. Thus we do not reproduce the
distinct double cavity. This may well be remedied in simulations with, for
example, ellipsoidal density profiles. On the other hand, the double X-ray
cavity associated with the Phoenix galaxy cluster studied by (Akahori et al.,
2020) consists of two well defined depressions that suggest an interpretation
in terms of a low Mach number jet such as shown in Fig 9.

A number of cavities at low redshift do possess a double bubble structure with
the LoFAR telescope revealing the driving radio galaxy (Bîrzan et al., 2020).
Examples are A1361 and MS0375). The single ovoid or thimble structures
predicted here from the high Mach number jets are also found in this survey
(e.g. A2052, Zw8276) although less definitive. There is thus the indication
that double bubbles are associated with samples derived from cluster surveys
while ovals will be found ina association with powerful radio galaxies.

A cluster evolution scenario in which AGN feedback regulates the intracluster
medium and so controls the development of central galaxies has received much
favourable attention (Fabian, 2012; Hardcastle et al., 2019). Here, we find
quantitative support for the radio-mode in which jet-induced outbursts provide
repeated heating events in the environment (Gupta et al., 2020). Although we
find that the jet power is effectively transferred into the cluster gas, one
10–50 Myr outburst is only sufficient to provide a simmering of the
immediately surrounding gas. Nevertheless, this is adequate to raise the
energy in the gas which can be transported through weak shocks and
magnetohydrodynamic waves into the entire cluster medium on a ten times longer
time scale in between the outbursts. Repeated bursts then act to raise the
energy until the simmering gas is able to effectively cool and a cooling
inflow develops. Without the repeated gentle simmering, the radio burst would
blow away the cluster medium. This implies that the cluster cooling flows we
now see as strong X-ray sources have reached the point of balance between the
simmering process and the radiative cooling.

##  5 Conclusions

We have presented simulated X-ray images of radio galaxies formed through
supersonic hydrodynamic jets. Only synchrotron and thermal Bremsstrahlung
processes have been considered without consideration for radiative losses or
particle acceleration. Three dimensional simulations are necessary not only
for the precession but also for the correct form of dissipation of turbulent
energy through cascade downwards in wavelength.

Standard Mach 6 jets drive shells and cavities of roughly ellipsoidal shape
and show a decrement in surface brightness of at most a few tens of per cent.
We have found that both cavities/shells and bows/arcs are generated in the
soft and hard X-ray bands, respectively.

Straight FR II jets are associated with elliptical shells when the jet-ambient
density ratio is high and with (projected as triangular) conical X-ray lobes
when the jets are very light. Remarkable thick rib structures are evident when
the jet-ambient density ratio is near 0.01.

Precession generates bi-lobed or X-shaped radio sources and in curved or
deflected radio jets. In these cases, the strongest X-ray emission is found to
be not associated with the hotspot but with the relic lobe or deflection
location. This is because the regions of hot high-pressure gas and regions of
dense high-compression do not coincide. As the viewing angle is changed
towards the line of sight, the cavity becomes a deep round hole surrounded by
circular ripples.

Low Mach number shocks extend into the cluster environment but a single
episode of radio galaxy formation is not sufficient for long-term support of a
cooling flow. We have found that the total radio luminosity of a simulated
radio galaxy does not increase with source age or size. This is because the
integrated emission is dominated by the higher pressure regions which are of
relatively constant volume, rather than the low pressure lobes. The jet
kinetic power is mainly transferred into the ambient medium via shocks and
compression as shown in Paper 1. In contrast, the X-ray luminosity steadily
increases, especially in the higher energy bands, until ambient gas is
expelled from the simulated region.

A new result is the appearance of a rib cage in soft X-rays which may be
related to the structure found in Cygnus A as presented in detail from Chandra
observations (Duffy et al., 2018). We have only considered thermal X-ray
emission in this work and so have investigated the shape of cavities. It
should also be noted that inverse Compton losses to photons from the Cosmic
Microwave Background may well dominate the X-ray image as well as alter the
radio structure and the morphological type as the relative contributions of
the hotspot and diffuse lobe may systematically change (Ker et al., 2012).
Furthermore, a non-thermal X-ray jet in Cygnus A is mis-aligned from the radio
jet suggesting a precession scenario (Steenbrugge et al., 2008). Relativistic
jet flows are generally expected to yield similar flow patterns but light
travel time effects may introduce new details. However, X-ray properties
should not be altered. As noted above, the relativistic electrons which
generate the radio emission are a source of energy which is transferred to
cosmic microwave background photons via inverse Compton scattering throughout
the lobes. This can fill in the cavities and reduce the number of radio
galaxies which display cavities.

Finally, using 3C 310 as an example, we demonstrate that a radio galaxy is
efficient at simmering the intracluster medium. A periodic gentle heating over
several short episodes could support the external medium through the
propagation of weak shocks.

These conclusions are founded on several model assumptions. We assume a
uniform ambient medium which will evidently never be the case and rarely valid
as an approximation. The density within isolated galaxies and clusters are
observed to decrease with distance from the bulge or core for nearby examples
(Croston et al., 2008). Moreover, for the distant radio galaxies being
considered here, the propensity of progressing galaxy and cluster mergers are
found to significantly disturb the density profiles and X-ray structures
(Bartalucci et al., 2017).

Adopting a non-uniform radial profile would have a profound influence on
precessing jets as has been found for the near straight jet variety (Clarke et
al., 1997). The jet would progress from light to heavy and the corresponding
dominant radio structure change from the artist’s brush to a ballistic pattern
in which the jet path is delineated. For the external X-ray emission, the
cocoon would become broader ad the radio galaxy penetrates into less dense
fields. This reduces the emission fron the bow caps and creates considerably
more open X-ray cavity structures, as often observed and simulated (Krause,
2005).

In addition, it should be reiterated that we also assume no active influence
from a magnetic field in either medium, non-relativistic speeds, no shear or
spray on exit from the jet nozzle and an initial thermal pressure balance on
entry. All these assumptions will be tested in future works.

##  6 Acknowledgements

We wish to thank SEPnet for supplying infrastructure funding.

##  7 Data availability

Data available on request. The data underlying this article will be shared on
reasonable request to the corresponding author.

## References

  * Akahori et al. (2020) Akahori T., et al., 2020, PASJ, 72, 62
  * Atoyan & Dermer (2004) Atoyan A., Dermer C. D., 2004, ApJ, 613, 151
  * Bartalucci et al. (2017) Bartalucci I., Arnaud M., Pratt G. W., Démoclès J., van der Burg R. F. J., Mazzotta P., 2017, A&A, 598, A61
  * Belsole et al. (2007) Belsole E., Worrall D. M., Hardcastle M. J., Croston J. H., 2007, MNRAS, 381, 1109
  * Bîrzan et al. (2004) Bîrzan L., Rafferty D. A., McNamara B. R., Wise M. W., Nulsen P. E. J., 2004, ApJ, 607, 800
  * Bîrzan et al. (2020) Bîrzan L., et al., 2020, MNRAS, 496, 2613
  * Boehringer et al. (1993) Boehringer H., Voges W., Fabian A. C., Edge A. C., Neumann D. M., 1993, MNRAS, 264, L25
  * Carilli et al. (1994) Carilli C. L., Perley R. A., Harris D. E., 1994, MNRAS, 270, 173
  * Choi et al. (2004) Choi Y.-Y., Reynolds C. S., Heinz S., Rosenberg J. L., Perlman E. S., Yang J., 2004, ApJ, 606, 185
  * Clarke et al. (1997) Clarke D. A., Harris D. E., Carilli C. L., 1997, MNRAS, 284, 981
  * Croston et al. (2008) Croston J. H., et al., 2008, A&A, 487, 431
  * Croston et al. (2011) Croston J. H., Hardcastle M. J., Mingo B., Evans D. A., Dicken D., Morganti R., Tadhunter C. N., 2011, ApJ, 734, L28
  * Donohoe & Smith (2016) Donohoe J., Smith M. D., 2016, MNRAS, 458, 558
  * Duffy et al. (2018) Duffy R. T., et al., 2018, MNRAS, 476, 4848
  * Fabian (2012) Fabian A. C., 2012, ARA&A, 50, 455
  * Giacintucci et al. (2011) Giacintucci S., et al., 2011, ApJ, 732, 95
  * Guo (2015) Guo F., 2015, ApJ, 803, 48
  * Gupta et al. (2020) Gupta N., et al., 2020, MNRAS,
  * Hardcastle & Krause (2013) Hardcastle M. J., Krause M. G. H., 2013, MNRAS, 430, 174
  * Hardcastle et al. (2019) Hardcastle M. J., et al., 2019, A&A, 622, A12
  * Jaffe & Perola (1974) Jaffe W., Perola G. C., 1974, A&A, 31, 223
  * Kataoka & Stawarz (2005) Kataoka J., Stawarz Ł., 2005, ApJ, 622, 797
  * Ker et al. (2012) Ker L. M., Best P. N., Rigby E. E., Röttgering H. J. A., Gendre M. A., 2012, MNRAS, 420, 2644
  * Kraft et al. (2012) Kraft R. P., et al., 2012, ApJ, 749, 19
  * Krause (2005) Krause M., 2005, A&A, 431, 45
  * Krause et al. (2019) Krause M. G. H., et al., 2019, MNRAS, 482, 240
  * Mendygral et al. (2011) Mendygral P. J., O’Neill S. M., Jones T. W., 2011, ApJ, 730, 100
  * Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228
  * Pandge et al. (2012) Pandge M. B., Vagshette N. D., David L. P., Patil M. K., 2012, MNRAS, 421, 808
  * Pandge et al. (2019) Pandge M. B., Sonkamble S. S., Parekh V., Dabhade P., Parmar A., Patil M. K., Raychaudhury S., 2019, ApJ, 870, 62
  * Prasow-Émond et al. (2020) Prasow-Émond M., et al., 2020, AJ, 160, 103
  * Riffel et al. (2019) Riffel R. A., et al., 2019, MNRAS, 485, 5590
  * Rybicki & Lightman (1979) Rybicki G. B., Lightman A. P., 1979, Radiative processes in astrophysics
  * Shin et al. (2016) Shin J., Woo J.-H., Mulchaey J. S., 2016, ApJS, 227, 31
  * Smith & Donohoe (2019) Smith M. D., Donohoe J., 2019, MNRAS, 490, 1363
  * Smith et al. (1985) Smith M. D., Norman M. L., Winkler K.-H. A., Smarr L., 1985, MNRAS, 214, 67
  * Smith et al. (2002) Smith D. A., Wilson A. S., Arnaud K. A., Terashima Y., Young A. J., 2002, ApJ, 565, 195
  * Sonkamble et al. (2015) Sonkamble S. S., Vagshette N. D., Pawar P. K., Patil M. K., 2015, Ap&SS, 359, 21
  * Stawarz et al. (2014) Stawarz Ł., et al., 2014, ApJ, 794, 164
  * Steenbrugge et al. (2008) Steenbrugge K. C., Blundell K. M., Duffy P., 2008, MNRAS, 388, 1465
  * Toro et al. (1994) Toro E. F., Spruce M., Speares W., 1994, Shock Waves, 4, 25
  * Tucker (1975) Tucker W., 1975, Radiation processes in astrophysics
  * Turner & Shabala (2020) Turner R. J., Shabala S. S., 2020, MNRAS,
  * Vagshette et al. (2017) Vagshette N. D., Naik S., Patil M. K., Sonkamble S. S., 2017, MNRAS, 466, 2054
  * Vagshette et al. (2019) Vagshette N. D., Naik S., Patil M. K., 2019, MNRAS, 485, 1981
  * Wilson & Yang (2002) Wilson A. S., Yang Y., 2002, ApJ, 568, 133
  * Wilson et al. (2006) Wilson A. S., Smith D. A., Young A. J., 2006, ApJ, 644, L9
  * Yarza et al. (2020) Yarza R., Wong G. N., Ryan B. R., Gammie C. F., 2020, ApJ, 898, 50


# Target Clustering Based Multi-Bernoulli Filter for Superpositional Sensors

Wang Sen This paragraph of the first footnote will contain the date on which
you submitted your paper for review. _(Corresponding author: Wang Sen)_ The
authors are with National University of Defense Technology, Xi’an 710106,
China (e-mail<EMAIL_ADDRESS>
###### Abstract

The sensor whose output is a function of the sum of contributions from targets
present in the surveillance area is called superpositional sensor. In this
letter, target clustering based multi-Bernoulli filter for superpositional
sensors is proposed. Targets are clustered according to the set of resolution
cells illuminated by them. Single target posterior density is strictly
derived, and densities of all the targets are combined to a approximate multi-
target posterior, which makes the multi-Bernoulli density is conjugate with
respect to the likelihood of superpositional sensors. The Gaussian
implementation of the proposed algorithm is also presented, where the
multidimensionality and the nonlinearity of update equation are handled by
sigma point transformation. The simulation results illustrate that the
proposed algorithm is effective confronted with the interaction of multiple
targets and long term overlapping of two targets.

###### Index Terms:

Multi-Bernoulli filter, superpositional sensors, target clustering, sigma
point transformation, track before detect (TBD).

##  I Introduction

Multi-target tracking (MTT) is concerned with estimating the number of targets
and their individual states from a sequence of measurements. The detect before
track (DBT) method approaches point measurements obtained by a detection
process, and may not be adequate for applications with low signal to noise
ratio (SNR). The track before detect (TBD) paradigm, which directly processes
raw measurements, can make use of all information contained in image
observations to improve tracking performance. When tracking multiple closely
spaced targets, regions of target influence in image observations are
overlapping, and the sensor output is a function of the sum of contributions
from individual targets, which is called the superpositional sensor. Examples
belonging to superpositional sensors are direction of arrival estimation for
antenna array [1], acoustic amplitude sensor [2], Doppler radar [3], and radio
frequency tomography [4].

For superpositional sensors, the update equation of cardinalized probability
hypothesis density (CPHD) filter is derived [5], and it is theoretically
rigorous but computationally intractable. Based on probability hypothesis
density (PHD) approximation, the additive likelihood moment (ALM) filter and
its particle implementation are designed [6]. And then, the ALM filter is
generalized to CPHD filtering recursion [7, 8]. Unlike the PHD/CPHD
recursions, the multi-Bernoulli recursion approximately propagates the full
multi-target distribution, and avoids postprocessing procedures such as
clustering to obtain individual target state. Under the assumption that the
regions of the observation influenced by individual targets do not overlap,
the multi-Bernoulli-TBD (MB-TBD) filters for various prior distributions are
derived [9]. In the context of superpositional measurements, the multi-
Bernoulli prior is not conjugate, and approximation is inevitable. In [10,
11], the multi-Bernoulli posterior distribution is propagated by updating
predefined conditional PHD. In [12], to obtain the single target likelihood,
other interacting targets are treated as interferences. In [13], the predicted
measurement information of each target is shared with the rest of the targets.

In this letter, target clustering based multi-Bernoulli (TC-MB) filter for
superpositional sensors is proposed. The set of resolution cells illuminated
by the target is defined. Accordingly, the targets are clustered, and each
cluster will be updated jointly. The single target posterior density is found
by integrating out the other targets, and it is a Bernoulli density. Combining
all the Bernoulli densities produces a multi-Bernoulli density, which can
approximate the real multi-target posterior. In the Gaussian implementation,
the multidimensionality and the nonlinearity of update equation are handled by
sigma point transformation.

##  II Background

The multi-target motion model based on random finite sets (RFS)
comprehensively considers target birth, death, spawning, and motion. Suppose
the multi‐target state at time $k-1$ is represented as a RFS
${{\bf{X}}_{k-1}}=\left\\{{{\bf{x}}_{k-1}^{1},{\bf{x}}_{k-1}^{2},\cdots,{\bf{x}}_{k-1}^{{M_{k-1}}}}\right\\}$,
where ${M_{k-1}}$ is the number of targets, regardless of spawned targets, the
multi-target state at time $k$ is ${{\bf{X}}_{k}}=\cup
S\left({{\bf{x}}_{k-1}^{1}}\right)\cdots\cup
S\left({{\bf{x}}_{k-1}^{{M_{k-1}}}}\right)\cup{\Gamma_{k}}$, where
$S\left({{\bf{x}}_{k-1}^{m}}\right)$ is the state of target
${\bf{x}}_{k-1}^{m}$ at time $k$, $m=1,2,\cdots,{M_{k-1}}$, determined by the
survival probability ${p_{S}}\left({{\bf{x}}_{k-1}^{m}}\right)$ and single
target Markov transition density
$f\left({{\bf{x}}_{k}^{m}|{\bf{x}}_{k-1}^{m}}\right)$, and ${\Gamma_{k}}$ is
the RFS of completely new targets.

At time $k$, the measurement provided by a superpositional sensor is a random
array ${{\bf{Z}}_{k}}=\left[{z_{k}^{n}}\right]$, where $n=1,2,\cdots,N$, with
$N$ being the number of cells. One resolution cell could be illuminated by
multiple targets. The set of targets contributing to cell $n$ is denoted as
$\Theta\left(n\right)$. Consequently, the measurement $z_{k}^{n}$ has the
following expression

$z_{k}^{n}=\left\\{{\begin{array}[]{*{20}{l}}{\sum\limits_{{\bf{x}}\in\Theta\left(n\right)}{{h^{n}}\left({\bf{x}}\right)}+{w^{n}}}&{\Theta\left(n\right)\neq\emptyset}\\\
{{w^{n}}}&{\Theta\left(n\right)=\emptyset}\end{array}}\right.$ (1)

where ${h^{n}}\left({\bf{x}}\right)$ represents the contribution of target
${\bf{x}}$ to the measurement of cell $n$, and ${w^{n}}$ is zero-mean Gaussian
white noise with covariance $R$.

Define
${{\bf{Z}}_{1:k}}=\left[{{{\bf{Z}}_{1}},{{\bf{Z}}_{2}},\cdots,{{\bf{Z}}_{k}}}\right]$,
the posterior probability density
$p\left({{{\bf{X}}_{k}}|{{\bf{Z}}_{1:k}}}\right)$ can be computed using Bayes
rule

$p\left({{{\bf{X}}_{k}}|{{\bf{Z}}_{1:k-1}}}\right)=\int{p\left({{{\bf{X}}_{k}}|{{\bf{X}}_{k-1}}}\right)p\left({{{\bf{X}}_{k-1}}|{{\bf{Z}}_{1:k-1}}}\right)\delta{{\bf{X}}_{k-1}}}$
(2)
$p\left({{{\bf{X}}_{k}}|{{\bf{Z}}_{1:k}}}\right)=\frac{{p\left({{{\bf{Z}}_{k}}|{{\bf{X}}_{k}}}\right)p\left({{{\bf{X}}_{k}}|{{\bf{Z}}_{1:k-1}}}\right)}}{{\int{p\left({{{\bf{Z}}_{k}}|{{\bf{X}}_{k}}}\right)p\left({{{\bf{X}}_{k}}|{{\bf{Z}}_{1:k-1}}}\right)\delta{{\bf{X}}_{k}}}}}$
(3)

A multi-Bernoulli RFS ${\bf{X}}$ is the disjoint union of a fixed number of
independent Bernoulli RFSs, and has density

${p_{{\rm{MB}}}}\left({\bf{X}}\right)=\sum\limits_{{{\bf{X}}^{1}}\uplus\cdots\uplus{{\bf{X}}^{M}}={\bf{X}}}{\prod\limits_{m=1}^{M}{{f^{m}}\left({{{\bf{X}}^{m}}}\right)}}$
(4)

where $\uplus$ denotes disjoint union, and the Bernoulli RFS ${{\bf{X}}^{m}}$
has density ${f^{m}}$ with existence probability ${r^{m}}$ and existence-
conditioned probability density function (PDF) ${p^{m}}\left({\bf{x}}\right)$,
$m=1,\cdots,M$.

##  III Target Clustering Based Multi-Bernoulli Filter

The proposed target clustering based multi-Bernoulli filter for
superpositional sensors propagates over time the parameters of the multi-
Bernoulli RFS representing the multi-target state.

###  III-A Multi-Bernoulli Filter Prediction

The superpositional sensors do not change the multi-Bernoulli prediction
equations [14, 9]. Suppose the posterior multi-Bernoulli density at time $k-1$
is characterized by parameters
$\left({r_{k-1|k-1}^{m},p_{k-1|k-1}^{m}\left({\bf{x}}\right)}\right),m=1,\cdots,{M_{k-1|k-1}}$
and completely new targets at time $k$ is characterized by parameters
$\left({r_{k}^{b,m},p_{k}^{b,m}\left({\bf{x}}\right)}\right),m=1,\cdots,M_{k}^{b}$,
the prior multi-Bernoulli density at time $k$ can be written as

$\begin{array}[]{l}p\left({{{\bf{X}}_{k}}|{{\bf{Z}}_{1:k-1}}}\right)=\\\
\sum\limits_{{\bf{X}}_{k|k-1}^{1}\uplus\cdots\uplus{\bf{X}}_{k|k-1}^{{M_{k|k-1}}}={{\bf{X}}_{k}}}{\prod\limits_{m=1}^{{M_{k|k-1}}}{f_{k|k-1}^{m}\left({{\bf{X}}_{k|k-1}^{m}}\right)}}\end{array}$
(5)

where ${M_{k|k-1}}={M_{k-1|k-1}}+M_{k}^{b}$ is the number of Bernoulli
components, and the Bernoulli RFS ${\bf{X}}_{k|k-1}^{m}$ has density

$f_{k|k-1}^{m}\left({{\bf{X}}_{k|k-1}^{m}}\right)=\left\\{{\begin{array}[]{*{20}{l}}{1-r_{k|k-1}^{m}}&{{\bf{X}}_{k|k-1}^{m}=\emptyset}\\\
{r_{k|k-1}^{m}p_{k|k-1}^{m}\left({\bf{x}}\right)}&{{\bf{X}}_{k|k-1}^{m}=\left\\{{\bf{x}}\right\\}}\\\
0&{\left|{{\bf{X}}_{k|k-1}^{m}}\right|>1}\end{array}}\right.$ (6)

For the surviving Bernoulli components, $m=1,\cdots,{M_{k-1|k-1}}$, the
parameters are

$r_{k|k-1}^{m}=r_{k-1|k-1}^{m}\left\langle{p_{k-1|k-1}^{m},{p_{S}}}\right\rangle$
(7)
$p_{k|k-1}^{m}\left({\bf{x}}\right)=\frac{{\int{p_{k-1|k-1}^{m}\left({\bf{x}}\right){p_{S}}\left({\bf{x}}\right)f\left({\cdot|{\bf{x}}}\right)d{\bf{x}}}}}{{\left\langle{p_{k-1|k-1}^{m},{p_{S}}}\right\rangle}}$
(8)

where
$\left\langle{a,b}\right\rangle=\int{a\left({\bf{x}}\right)b\left({\bf{x}}\right)d{\bf{x}}}$
represents the inner product of $a\left({\bf{x}}\right)$ and
$b\left({\bf{x}}\right)$. For Bernoulli components of completely new targets,
$m={M_{k-1|k-1}}+1,\cdots,{M_{k|k-1}}$, the parameters are

$r_{k|k-1}^{m}=r_{k}^{b,m-{M_{k-1|k-1}}}$ (9)
$p_{k|k-1}^{m}\left({\bf{x}}\right)=p_{k}^{b,m-{M_{k-1|k-1}}}\left({\bf{x}}\right)$
(10)

###  III-B Target Clustering

For the superpositional measurement model, the targets are interacting with
each other in the measurement space inevitably, and the corresponding single
target likelihood is not available. In the proposed algorithm, the targets are
clustered, and each cluster will be updated jointly.

The set of resolution cells illuminated by target ${\bf{X}}_{k|k-1}^{m}$ is
denoted as $\Phi\left({{\bf{X}}_{k|k-1}^{m}}\right)$, for example
$\Phi\left({{\bf{X}}_{k|k-1}^{m}}\right)$ could be the set of cells whose
centers fall within a certain distance from the position of the target. If
there exists resolution cell ${n^{*}}$ being illuminated by target
${\bf{X}}_{k|k-1}^{{m_{1}}}$ and ${\bf{X}}_{k|k-1}^{{m_{2}}}$, namely
${n^{*}}\in\Phi\left({{\bf{X}}_{k|k-1}^{{m_{1}}}}\right)\cap\Phi\left({{\bf{X}}_{k|k-1}^{{m_{2}}}}\right)$,
target ${\bf{X}}_{k|k-1}^{{m_{1}}}$ and ${\bf{X}}_{k|k-1}^{{m_{2}}}$ are
clustered into one cluster, and they are called the interacting targets;
otherwise, target ${\bf{X}}_{k|k-1}^{{m_{1}}}$ and
${\bf{X}}_{k|k-1}^{{m_{2}}}$ are clustered into two different clusters, and
they are called the well-separated targets. Accordingly, multiple targets
${{\bf{X}}_{k}}$ will be divided into ${C_{k}}$ clusters

${{\bf{X}}_{k}}={\bf{X}}_{k}^{\left(1\right)}\uplus\cdots{\bf{X}}_{k}^{\left(c\right)}\uplus\cdots{\bf{X}}_{k}^{\left({{C_{k}}}\right)}$
(11)

One cluster may contain more than two targets. For example, if
$\Phi\left({{\bf{X}}_{k|k-1}^{{m_{1}}}}\right)\cap\Phi\left({{\bf{X}}_{k|k-1}^{{m_{2}}}}\right)\neq\emptyset$
and
$\Phi\left({{\bf{X}}_{k|k-1}^{{m_{1}}}}\right)\cap\Phi\left({{\bf{X}}_{k|k-1}^{{m_{3}}}}\right)\neq\emptyset$,
target ${\bf{X}}_{k|k-1}^{{m_{1}}}$, ${\bf{X}}_{k|k-1}^{{m_{2}}}$, and
${\bf{X}}_{k|k-1}^{{m_{3}}}$ belong to one cluster. One cluster also may
contain one target, when it is well separated from others. Naturally, the set
of resolution cells illuminated by targets in ${\bf{X}}_{k}^{\left(c\right)}$
is defined as

$\Phi\left({{\bf{X}}_{k}^{\left(c\right)}}\right)\buildrel\Delta\over{=}\bigcup\limits_{{\bf{X}}_{k|k-1}^{m}\in{\bf{X}}_{k}^{\left(c\right)}}{\Phi\left({{\bf{X}}_{k|k-1}^{m}}\right)},c=1,\cdots,{C_{k}}$
(12)

and the set of resolution cells illuminated by targets at time $k$ is defined
as

$\Phi\left({{{\bf{X}}_{k}}}\right)\buildrel\Delta\over{=}\bigcup\limits_{{\bf{X}}_{k|k-1}^{m}\in{{\bf{X}}_{k}}}{\Phi\left({{\bf{X}}_{k|k-1}^{m}}\right)}$
(13)

The multi-target likelihood function can be written as

$\begin{array}[]{l}p\left({{{\bf{Z}}_{k}}|{{\bf{X}}_{k}}}\right)=\prod\limits_{n=1}^{N}{{\rm{{\cal
N}}}\left({z_{k}^{n};0,R}\right)}\times\\\
\prod\limits_{c=1}^{{C_{k}}}{\prod\limits_{n\in\Phi\left({{\bf{X}}_{k}^{\left(c\right)}}\right)}{\frac{{{\rm{{\cal
N}}}\left({z_{k}^{n};\sum\limits_{{\bf{x}}\in\Theta\left(n\right)}{{h^{n}}\left({\bf{x}}\right)},R}\right)}}{{{\rm{{\cal
N}}}\left({z_{k}^{n};0,R}\right)}}}}\end{array}$ (14)

where ${\rm{{\cal N}}}\left({\cdot;\mu,\Sigma}\right)$ denotes the probability
density function of the Gaussian distribution with mean $\mu$ and covariance
matrix $\Sigma$.

###  III-C Multi-Bernoulli Filter Update

Substituting (5) and (14) into (3) produces the multi-target posterior density

$p\left({{{\bf{X}}_{k}}|{{\bf{Z}}_{1:k}}}\right)\propto\sum\limits_{{\bf{X}}_{k|k-1}^{1}\uplus\cdots\uplus{\bf{X}}_{k|k-1}^{{M_{k|k-1}}}={{\bf{X}}_{k}}}{\prod\limits_{c=1}^{{C_{k}}}{p\left({{\bf{X}}_{k}^{\left(c\right)}|{{\bf{Z}}_{1:k}}}\right)}}$
(15)

where the posterior density of targets in ${\bf{X}}_{k}^{\left(c\right)}$ is

$\begin{array}[]{l}p\left({{\bf{X}}_{k}^{\left(c\right)}|{{\bf{Z}}_{1:k}}}\right)\propto\\\
\prod\limits_{{\bf{X}}_{k|k-1}^{m}\in{\bf{X}}_{k}^{\left(c\right)}}{f_{k|k-1}^{m}\left({{\bf{X}}_{k|k-1}^{m}}\right)}\prod\limits_{n\in\Phi\left({{\bf{X}}_{k}^{\left(c\right)}}\right)}{\frac{{{\rm{{\cal
N}}}\left({z_{k}^{n};\sum\limits_{{\bf{x}}\in\Theta\left(n\right)}{{h^{n}}\left({\bf{x}}\right)},R}\right)}}{{{\rm{{\cal
N}}}\left({z_{k}^{n};0,R}\right)}}}\end{array}$ (16)

Equation (15) indicates that the multi-target posterior density is the product
of the posterior density of targets in each cluster.

The multi-target likelihood ratio in (16) can be expressed as

${L^{\left(c\right)}}\left({{b_{1}},\cdots,{b_{\left|{{\bf{X}}_{k}^{\left(c\right)}}\right|}}}\right)=\prod\limits_{n\in\Phi\left({{\bf{X}}_{k}^{\left(c\right)}}\right)}{\frac{{{\rm{{\cal
N}}}\left({z_{k}^{n};\sum\limits_{{\bf{x}}\in\Theta\left(n\right)}{{h^{n}}\left({\bf{x}}\right)},R}\right)}}{{{\rm{{\cal
N}}}\left({z_{k}^{n};0,R}\right)}}}$ (17)

All the arguments in ${L^{\left(c\right)}}$ are booleans, representing the
existence of targets. For example, ${b_{1}}=1$ indicates that the first target
in ${\bf{X}}_{k}^{\left(c\right)}$ exists, while ${b_{1}}=0$ means that the
first target in ${\bf{X}}_{k}^{\left(c\right)}$ does not exist.

The posterior density of target ${\bf{X}}_{k}^{m}$ (${\bf{X}}_{k|k-1}^{m}$) is
found by integrating out the other targets

$p\left({{\bf{X}}_{k}^{m}|{{\bf{Z}}_{1:k}}}\right)=\int{p\left({{\bf{X}}_{k}^{\left(c\right)}|{{\bf{Z}}_{1:k}}}\right)\delta{\bf{X}}_{k}^{\left(c\right)\backslash
m}}$ (18)

where ${\bf{X}}_{k}^{\left(c\right)\backslash
m}={\bf{X}}_{k}^{\left(c\right)}\backslash{\bf{X}}_{k}^{m}$.

The single target posterior density (18) is a Bernoulli density, and the
multi-target posterior density (15) can be approximated by the following
multi-Bernoulli density

$p\left({{{\bf{X}}_{k}}|{{\bf{Z}}_{1:k}}}\right)\approx\sum\limits_{{\bf{X}}_{k}^{1}\uplus\cdots\uplus{\bf{X}}_{k}^{{M_{k}}}={{\bf{X}}_{k}}}{\prod\limits_{m=1}^{{M_{k}}}{p\left({{\bf{X}}_{k}^{m}|{{\bf{Z}}_{1:k}}}\right)}}$
(19)

where ${M_{k}}={M_{k|k-1}}$. Thus, it is achieved that the multi-Bernoulli
density is conjugate with respect to the likelihood of superpositional
sensors.

##  IV Gaussian Implementation

This section describes the Gaussian implementation of the proposed algorithm.

###  IV-A Prediction

The Gaussian implementation of prediction is obtained when there are a
constant survival probability
${p_{S}}\left({{\bf{x}}_{k-1}^{m}}\right)\equiv{p_{S}}$ and Gaussian/linear
single target Markov transition density

$f\left({{\bf{x}}_{k}^{m}|{\bf{x}}_{k-1}^{m}}\right)={\rm{{\cal
N}}}\left({{\bf{x}}_{k}^{m};{\bf{Fx}}_{k-1}^{m},{\bf{Q}}}\right)$ (20)

where ${\bf{F}}$ and ${\bf{Q}}$ are the state transition matrix and the
process noise covariance, respectively. Suppose the posterior state density of
the Bernoulli RFS ${\bf{X}}_{k-1|k-1}^{m},m=1,\cdots,{M_{k-1|k-1}}$ is
approximated as

$p_{k-1|k-1}^{m}\left({\bf{x}}\right)={\rm{{\cal
N}}}\left({{\bf{x}};{\bf{\bar{x}}}_{k-1|k-1}^{m},{\bf{P}}_{k-1|k-1}^{m}}\right)$
(21)

the prior mean and covariance of surviving targets at time $k$ are provided by
the Kalman filter prediction

${\bf{\bar{x}}}_{k|k-1}^{m}={\bf{F\bar{x}}}_{k-1|k-1}^{m}$ (22)
${\bf{P}}_{k|k-1}^{m}={\bf{FP}}_{k-1|k-1}^{m}{{\bf{F}}^{\rm{T}}}+{\bf{Q}}$
(23)

which constitute the prior state density

$p_{k|k-1}^{m}\left({\bf{x}}\right)={\rm{{\cal
N}}}\left({{\bf{x}};{\bf{\bar{x}}}_{k|k-1}^{m},{\bf{P}}_{k|k-1}^{m}}\right)$
(24)

and the prior existence probability is

$r_{k|k-1}^{m}=r_{k-1|k-1}^{m}{p_{S}}$ (25)

###  IV-B Target Clustering

In the implementation of target clustering, the set of resolution cells
illuminated by target ${\bf{X}}_{k|k-1}^{m}$ can be specified by the prior
mean ${\bf{\bar{x}}}_{k|k-1}^{m}$ when it exists, and be defined as

$\Phi\left({{\bf{X}}_{k|k-1}^{m}}\right)=\Phi\left({{\bf{\bar{x}}}_{k|k-1}^{m}}\right)$
(26)

According to III-B, all the prior targets can be divided into ${C_{k}}$
clusters, seen (11).

For clarity and conciseness, time index and cluster index of prior parameters
are omitted. Specifically, the number of targets in
${\bf{X}}_{k}^{\left(c\right)}$ is expressed as
$M=\left|{{\bf{X}}_{k}^{\left(c\right)}}\right|$, the Bernoulli RFS
constituting ${\bf{X}}_{k}^{\left(c\right)}$ is expressed as
${{\bf{X}}^{m}},m=1,\cdots,M$, and the corresponding existence probability,
state density, mean, and covariance are expressed as ${r^{m}}$,
${p^{m}}\left({{{\bf{x}}^{m}}}\right)$, ${{{\bf{\bar{x}}}}^{m}}$, and
${{\bf{P}}^{m}}$, respectively. Therefore, (16) can be expressed as

$\begin{array}[]{l}p\left({{\bf{X}}_{k}^{\left(c\right)}|{{\bf{Z}}_{1:k}}}\right)\propto\\\
\left\\{{\begin{array}[]{*{20}{l}}{\left({1-{r^{1}}}\right)\left({1-{r^{2}}}\right)\cdots\left({1-{r^{M}}}\right){L^{\left(c\right)}}\left({0,0,\cdots,0}\right)}\\\
{\left({1-{r^{1}}}\right)\left({1-{r^{2}}}\right)\cdots{r^{M}}{p^{M}}\left({{{\bf{x}}^{M}}}\right){L^{\left(c\right)}}\left({0,0,\cdots,1}\right)}\\\
\vdots\\\
{{r^{1}}{p^{1}}\left({{{\bf{x}}^{1}}}\right){r^{2}}{p^{2}}\left({{{\bf{x}}^{2}}}\right)\cdots{r^{M}}{p^{M}}\left({{{\bf{x}}^{M}}}\right){L^{\left(c\right)}}\left({1,1,\cdots,1}\right)}\end{array}}\right.\end{array}$
(27)

The condition of piecewise function (27) is omitted, because it can be easily
inferred by the arguments in ${L^{\left(c\right)}}$.

###  IV-C Update

The key of update is to obtain posterior target parameters using set integral
(18), which is computationally intractable. Using (27), (18) can be converted
into vector integral

$\begin{array}[]{l}p\left({{\bf{X}}_{k}^{m}|{{\bf{Z}}_{1:k}}}\right)\propto\\\
\left\\{{\begin{array}[]{*{20}{l}}{\left({1-{r^{1}}}\right)\left({1-{r^{2}}}\right)\cdots\left({1-{r^{M}}}\right){L^{\left(c\right)}}\left({0,0,\cdots,0}\right)}\\\
{\int{\left({1-{r^{1}}}\right)\left({1-{r^{2}}}\right)\cdots{r^{M}}{p^{M}}\left({{{\bf{x}}^{M}}}\right){L^{\left(c\right)}}\left({0,0,\cdots,1}\right)d{{\bf{x}}^{M}}}}\\\
\vdots\\\
{\int{{r^{1}}{p^{1}}\left({{{\bf{x}}^{1}}}\right)\cdots{r^{M}}{p^{M}}\left({{{\bf{x}}^{M}}}\right){L^{\left(c\right)}}\left({1,1,\cdots,1}\right)d{{\bf{x}}^{1}}\cdots
d{{\bf{x}}^{M}}}}\end{array}}\right.\end{array}$ (28)

It is worth noting that there is no $d{{\bf{x}}^{m}}$ in (28). Vector integral
(28) still encounters difficulties of the multidimensionality and the
nonlinearity, and this subsection uses sigma point transformation to implement
it.

The implementation is demonstrated with the last piece of (28), and the rest
pieces have simpler computations. Several sigma points and weights are
generated from the prior state density, and can approximate it as following

$\begin{array}[]{l}{p^{u}}\left({{{\bf{x}}^{u}}}\right)={\rm{{\cal
N}}}\left({{{\bf{x}}^{u}};{{{\bf{\bar{x}}}}^{u}},{{\bf{P}}^{u}}}\right)\approx\sum\limits_{{i_{u}}=1}^{{N_{S}}}{{w_{{i_{u}}}}\delta\left({{{\bf{x}}^{u}}-{\bf{\chi}}_{{i_{u}}}^{u}}\right)}\\\
u=1,\cdots,m-1,m+1,\cdots,M\end{array}$ (29)

where ${N_{S}}$ is the number of sigma points. The multi-target likelihood
ratio in the last piece of (28) can be represented as a function of prior
state variables

$l\left({{{\bf{x}}^{1}},\cdots,{{\bf{x}}^{m}},\cdots,{{\bf{x}}^{M}}}\right)\buildrel\Delta\over{=}{L^{\left(c\right)}}\left({1,1,\cdots,1}\right)$
(30)

Substituting (29) and (30) into the last piece of (28) produces

$\sum\limits_{{i_{1}}=1}^{{N_{S}}}{\cdots\sum\limits_{{i_{M}}=1}^{{N_{S}}}{{r^{1}}{w_{{i_{1}}}}\cdots{r^{M}}{w_{{i_{M}}}}l\left({{\bf{\chi}}_{{i_{1}}}^{1},\cdots,{\bf{\chi}}_{{i_{M}}}^{M}}\right)\delta\left({{{\bf{x}}^{m}}-{\bf{\chi}}_{{i_{m}}}^{m}}\right)}}$
(31)

which is a function of ${{\bf{x}}^{m}}$. Similarly, half pieces of (28) yield
a function of ${{\bf{x}}^{m}}$, and the other half pieces yield a scalar.
Suppose the sum of the former is denoted as
$\sum\nolimits_{i=1}^{{N_{S}}}{{{w^{\prime}}_{i}}\delta\left({{{\bf{x}}^{m}}-{\bf{\chi}}_{i}^{m}}\right)}$,
and the sum of the latter is denoted as ${r^{\prime}}$, then the posterior
density (28) has the form of Bernoulli density

$p\left({{\bf{X}}_{k}^{m}|{{\bf{Z}}_{1:k}}}\right)=\left\\{{\begin{array}[]{*{20}{l}}{1-\frac{{\sum\nolimits_{j=1}^{{N_{S}}}{{{w^{\prime}}_{j}}}}}{{r^{\prime}+\sum\nolimits_{j=1}^{{N_{S}}}{{{w^{\prime}}_{j}}}}}}&{{\bf{X}}_{k}^{m}=\emptyset}\\\
{\sum\limits_{i=1}^{{N_{S}}}{\frac{{{{w^{\prime}}_{i}}}}{{r^{\prime}+\sum\nolimits_{j=1}^{{N_{S}}}{{{w^{\prime}}_{j}}}}}\delta\left({{{\bf{x}}^{m}}-{\bf{\chi}}_{i}^{m}}\right)}}&{{\bf{X}}_{k}^{m}=\left\\{{{{\bf{x}}^{m}}}\right\\}}\end{array}}\right.$
(32)

Therefore, posterior parameters including existence probability, mean, and
covariance of the Bernoulli RFS ${\bf{X}}_{k}^{m}$ are

$r_{k|k}^{m}=\frac{{\sum\nolimits_{j=1}^{{N_{S}}}{{{w^{\prime}}_{j}}}}}{{r^{\prime}+\sum\nolimits_{j=1}^{{N_{S}}}{{{w^{\prime}}_{j}}}}}$
(33)
${\bf{\bar{x}}}_{k|k}^{m}=\frac{{\sum\nolimits_{i=1}^{{N_{S}}}{{{w^{\prime}}_{i}}{\bf{\chi}}_{i}^{m}}}}{{\sum\nolimits_{i=1}^{{N_{S}}}{{{w^{\prime}}_{i}}}}}$
(34)
${\bf{P}}_{k|k}^{m}=\frac{{\sum\nolimits_{i=1}^{{N_{S}}}{{{w^{\prime}}_{i}}\left({{\bf{\chi}}_{i}^{m}-{\bf{\bar{x}}}_{k|k}^{m}}\right){{\left({{\bf{\chi}}_{i}^{m}-{\bf{\bar{x}}}_{k|k}^{m}}\right)}^{\mathop{\rm
T}\nolimits}}}}}{{\sum\nolimits_{i=1}^{{N_{S}}}{{{w^{\prime}}_{i}}}}}$ (35)

##  V Simulation

In this section, the performance of the proposed target clustering based
multi-Bernoulli (TC-MB) filter for superpositional sensors is investigated and
compared with the multi-Bernoulli-TBD (MB-TBD) filter [9] in multi-target
tracking application.

In the simulations, the sensor field of view is a two-dimensional region
$\left[{0,128}\right]\times\left[{0,128}\right]$, and the total time of
simulation is $K=70$. There are five targets present in the sensor field of
view, and single target state is defined as
${\bf{x}}=\left[{x;\dot{x};y;\dot{y}}\right]$, representing target position
and velocity along the x-axis and y-axis. The target motion is modeled by a
nearly constant velocity model

${{\bf{x}}_{k}}={\bf{F}}{{\bf{x}}_{k-1}}+{{\bf{n}}_{k-1}}$ (36)
${\bf{F}}={{\bf{I}}_{2}}\otimes\left[{\begin{array}[]{*{20}{c}}1&T\\\
0&1\end{array}}\right]$ (37)

where ${{\bf{I}}_{2}}$ is the $2\times 2$ identity matrix, $T=1$ is the
sampling period, $\otimes$ represents the Kronecker product, and the
covariance of process noise ${{\bf{n}}_{k-1}}$ is

${\bf{Q}}={10^{-4}}\times{{\bf{I}}_{2}}\otimes\left[{\begin{array}[]{*{20}{c}}{0.25}&{0.5}\\\
{0.5}&1\end{array}}\right]$ (38)

Motion parameters of targets are shown in Table I, initial state of each
target obeys Gaussian distribution, and initial state covariance is set as
${\bf{Q}}$. At each time step, the measurement is a $128\times 128$ image
consisting of an array of cells with a scalar intensity, giving each cell side
lengths of $\Delta x=\Delta y=1$. The array index is treated as an ordered
pair of integers $n=\left({i,j}\right),1\leq i,j\leq 128$. The point spread
function in equation (1) is

${h^{n}}\left({\bf{x}}\right)=I\left({\bf{x}}\right)\times\exp\left({-\frac{{{{\left({i\Delta
x-x}\right)}^{2}}+{{\left({j\Delta
y-y}\right)}^{2}}}}{{\sigma_{h}^{2}}}}\right)$ (39)

where $I\left({\bf{x}}\right)$ is the constant intensity of target ${\bf{x}}$,
being $10$, $7$, $8$, $9$, and $10$ for the five targets respectively, and
$\sigma_{h}^{2}=2$ is blurring factor. When ${h^{n}}\left({\bf{x}}\right)$ is
greater than preset threshold ${h^{{\rm{th}}}}=1$, it is considered that cell
$n$ is illuminated by target ${\bf{x}}$. The covariance of measurement noise
${w^{n}}$ in equation (1) is $R=1$.

To ensure a fair comparison, the multi-Bernoulli-TBD filter [9] is implemented
by sigma points transformation, and it shares the same parameters with the
proposed target clustering based multi-Bernoulli filter for superpositional
sensors. The survival probability is ${p_{S}}=0.99$. At each time step, target
with posterior existence probability greater than ${r^{{\rm{th}}}}=0.99$ is
extracted, and that with posterior existence probability less than
${r^{\rm{p}}}={10^{-4}}$ is pruned. The number of sigma points is ${N_{S}}=9$,
and sigma points and weights are determined using [15]. In the Gaussian
implementation of the proposed TC-MB filter, the set of resolution cells
illuminated by prior target is determined by its point spread function and
preset threshold ${h^{{\rm{th}}}}$, which is used for target clustering.

Cardinality and Optimal Sub-Pattern Assignment (OSPA) [16] distance between
ground truth of multi-target state set and estimation of multi-target state
set are utilized as performance evaluation criterions of filters, where the
cut-off factor and the order used in OSPA are $c=10$ and $p=1$, respectively.

Simulation scenario of one trial is depicted in Fig.1. The resolution cell is
called superpositional position (the red circle) if it is illuminated by
multiple targets simultaneously. Near $k=21$, target 1, 2, 3, and 4 almost
reach the same cell (almost cell $\left({50,50}\right)$) at the same time. In
addition, the resolution cells illuminated by target 1 and 5 are overlapping
for a long time.

Fig.2 depicts cardinality and OSPA performances of MB-TBD filter and TC-MB
filter versus time over $1000$ Monte Carlo trials. Due to the interaction of
four targets at $k=21$, mean cardinality of MB-TBD filter drops dramatically
and corresponding mean OSPA deteriorates to $8$ until target 2,3, and 4 die.
In the second event, long term overlapping of two targets significantly
deteriorates mean OSPA of MB-TBD filter and increases standard deviation (std)
of its cardinality. The proposed TC-MB filter performs well throughout the
simulation.

TABLE I: Motion Parameters of Targets Target |  Mean of Initial State |  Birth Time |  Death Time
---|---|---|---
1 |  $\left[{20;1.5;20;1.5}\right]$ |  1 |  65
2 |  $\left[{80;-1.5;20;1.5}\right]$ |  1 |  40
3 |  $\left[{20;2;50;0}\right]$ |  6 |  38
4 |  $\left[{50;0;20;2}\right]$ |  6 |  36
5 |  $\left[{83;1.5;84;1.5}\right]$ |  43 |  60

Figure 1: Simulation scenario of one trial

Figure 2: Cardinality and OSPA performances of MB-TBD filter and TC-MB filter

##  VI Conclusion

In this letter, target clustering based multi-Bernoulli filter for
superpositional sensors is proposed. Single target posterior density is
strictly derived, and multi-target posterior is approximated by the product of
multiple Bernoulli densities. The Gaussian implementation of the proposed
algorithm is also presented. The simulation results illustrate that the
proposed algorithm can effectively track multiple targets when raw
measurements are superpositional.

## References

  * [1] X. Dong, J. Zhao, M. Sun, X. Zhang, and Y. Wang, “A modified $\delta$-generalized labeled multi-bernoulli filtering for multi-source doa tracking with coprime array,” _IEEE Transactions on Wireless Communications_ , pp. 1–1, 2023.
  * [2] A.-A. Saucan, Y. Li, and M. Coates, “Particle flow superpositional GLMB filter,” in _Signal Processing, Sensor/Information Fusion, and Target Recognition XXVI_ , I. Kadar, Ed., vol. 10200, International Society for Optics and Photonics. SPIE, 2017, p. 102000F. [Online]. Available: https://doi.org/10.1117/12.2263236
  * [3] W. Dong, K. Da, Z. Song, and Q. Fu, “Multi-target joint detection and estimation of doppler radar based on superpositional sensors,” _JOURNAL OF SIGNAL PROCESSING_ , vol. 38, no. 5, pp. 964–972, May 2022.
  * [4] I.-G. Ri, U.-G. Pak, S.-J. Han, and C.-J. Yun, “Crosshole radio electromagnetic tomography using apparent source intensity and logarithmic barrier,” _Journal of Applied Geophysics_ , vol. 212, p. 104987, 2023. [Online]. Available: https://www-sciencedirect-com-s.libyc.nudt.edu.cn:443/science/article/pii/S0926985123000654
  * [5] M. Ronald, “Cphd filters for superpositional sensors,” _Proc. SPIE 7445,Signal and Data Processing of Small Targets 2009, 74450E_ , vol. 7445, Sep. 2009.
  * [6] F. Thouin, S. Nannuru, and M. Coates, “Multi-target tracking for measurement models with additive contributions,” in _14th International Conference on Information Fusion_ , July 2011, pp. 1–8.
  * [7] M. Ronald and E.-F. Adel, “An approximate cphd filter for superpositional sensors,” in _Proc. SPIE 8392, Signal Processing, SensorFusion, and Target Recognition XXI, 83920K_ , May 2012.
  * [8] S. Nannuru, M. Coates, and R. Mahler, “Computationally-tractable approximate phd and cphd filters for superpositional sensors,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 7, no. 3, pp. 410–420, June 2013\.
  * [9] B.-N. Vo, B.-T. Vo, N.-T. Pham, and D. Suter, “Joint detection and estimation of multiple objects from image observations,” _IEEE Transactions on Signal Processing_ , vol. 58, no. 10, pp. 5129–5141, Oct 2010.
  * [10] S. Nannuru and M. Coates, “Multi-bernoulli filter for superpositional sensors,” in _Proceedings of the 16th International Conference on Information Fusion_ , July 2013, pp. 1632–1637.
  * [11] ——, “Particle filter implementation of the multi-bernoulli filter for superpositional sensors,” in _2013 5th IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP)_ , Dec 2013, pp. 368–371.
  * [12] D. Y. Kim, B. Ristic, R. Guan, and L. Rosenberg, “A bernoulli track-before-detect filter for interacting targets in maritime radar,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 57, no. 3, pp. 1981–1991, June 2021.
  * [13] E. S. Davies and . F. García-Fernández, “A multi-bernoulli gaussian filter for track-before-detect with superpositional sensors,” in _2022 25th International Conference on Information Fusion (FUSION)_ , July 2022, pp. 01–08.
  * [14] R. Mahler, _Statistical Multisource-Multitarget Information Fusion_. Artech, 2007. [Online]. Available: https://ieeexplore.ieee.org/document/9100524
  * [15] S. Julier and J. Uhlmann, “Unscented filtering and nonlinear estimation,” _Proceedings of the IEEE_ , vol. 92, no. 3, pp. 401–422, March 2004.
  * [16] 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.