# The Polarized Cosmic Hand: IXPE Observations of PSR B1509-58/MSH 15$-$52

Roger W. Romani Department of Physics and Kavli Institute for Particle
Astrophysics and Cosmology, Stanford University, Stanford, California 94305,
USA  Josephine Wong Department of Physics and Kavli Institute for Particle
Astrophysics and Cosmology, Stanford University, Stanford, California 94305,
USA  Niccoló Di Lalla Department of Physics and Kavli Institute for Particle
Astrophysics and Cosmology, Stanford University, Stanford, California 94305,
USA  Nicola Omodei Department of Physics and Kavli Institute for Particle
Astrophysics and Cosmology, Stanford University, Stanford, California 94305,
USA  Fei Xie Guangxi Key Laboratory for Relativistic Astrophysics, School of
Physical Science and Technology, Guangxi University, Nanning 530004, China
INAF Istituto di Astrofisica e Planetologia Spaziali, Via del Fosso del
Cavaliere 100, 00133 Roma, Italy  C.-Y. Ng Department of Physics, The
University of Hong Kong, Pokfulam, Hong Kong  Riccardo Ferrazzoli INAF
Istituto di Astrofisica e Planetologia Spaziali, Via del Fosso del Cavaliere
100, 00133 Roma, Italy  Alessandro Di Marco INAF Istituto di Astrofisica e
Planetologia Spaziali, Via del Fosso del Cavaliere 100, 00133 Roma, Italy
Niccoló Bucciantini INAF Osservatorio Astrofisico di Arcetri, Largo Enrico
Fermi 5, 50125 Firenze, Italy  Dipartimento di Fisica e Astronomia, Università
degli Studi di Firenze, Via Sansone 1, 50019 Sesto Fiorentino (FI), Italy
Istituto Nazionale di Fisica Nucleare, Sezione di Firenze, Via Sansone 1,
50019 Sesto Fiorentino (FI), Italy  Maura Pilia INAF Osservatorio Astronomico
di Cagliari, Via della Scienza 5, 09047 Selargius (CA), Italy  Patrick Slane
Center for Astrophysics — Harvard & Smithsonian, 60 Garden Street, Cambridge,
MA 02138, USA  Martin C. Weisskopf NASA Marshall Space Flight Center,
Huntsville, AL 35812, USA  Simon Johnston Australia Telescope National
Facility, CSIRO, Space and Astronomy, PO Box 76, Epping NSW 1710, Australia
Marta Burgay INAF Osservatorio Astronomico di Cagliari, Via della Scienza 5,
09047 Selargius (CA), Italy  Deng Wei Guangxi Key Laboratory for Relativistic
Astrophysics, School of Physical Science and Technology, Guangxi University,
Nanning 530004, China  Yi-Jung Yang Department of Physics & Laboratory for
Space Research, The University of Hong Kong, Pokfulam, Hong Kong  Shumeng
Zhang Department of Physics, The University of Hong Kong, Pokfulam, Hong Kong
Lucio A. Antonelli INAF Osservatorio Astronomico di Roma, Via Frascati 33,
00078 Monte Porzio Catone (RM), Italy  Space Science Data Center, Agenzia
Spaziale Italiana, Via del Politecnico snc, 00133 Roma, Italy  Matteo Bachetti
INAF Osservatorio Astronomico di Cagliari, Via della Scienza 5, 09047
Selargius (CA), Italy  Luca Baldini Istituto Nazionale di Fisica Nucleare,
Sezione di Pisa, Largo B. Pontecorvo 3, 56127 Pisa, Italy  Dipartimento di
Fisica, Università di Pisa, Largo B. Pontecorvo 3, 56127 Pisa, Italy  Wayne H.
Baumgartner NASA Marshall Space Flight Center, Huntsville, AL 35812, USA
Ronaldo Bellazzini Istituto Nazionale di Fisica Nucleare, Sezione di Pisa,
Largo B. Pontecorvo 3, 56127 Pisa, Italy  Stefano Bianchi Dipartimento di
Matematica e Fisica, Università degli Studi Roma Tre, Via della Vasca Navale
84, 00146 Roma, Italy  Stephen D. Bongiorno NASA Marshall Space Flight Center,
Huntsville, AL 35812, USA  Raffaella Bonino Istituto Nazionale di Fisica
Nucleare, Sezione di Torino, Via Pietro Giuria 1, 10125 Torino, Italy
Dipartimento di Fisica, Università degli Studi di Torino, Via Pietro Giuria 1,
10125 Torino, Italy  Alessandro Brez Istituto Nazionale di Fisica Nucleare,
Sezione di Pisa, Largo B. Pontecorvo 3, 56127 Pisa, Italy  Fiamma Capitanio
INAF Istituto di Astrofisica e Planetologia Spaziali, Via del Fosso del
Cavaliere 100, 00133 Roma, Italy  Simone Castellano Istituto Nazionale di
Fisica Nucleare, Sezione di Pisa, Largo B. Pontecorvo 3, 56127 Pisa, Italy
Elisabetta Cavazzuti ASI - Agenzia Spaziale Italiana, Via del Politecnico snc,
00133 Roma, Italy  Chien-Ting Chen Science and Technology Institute,
Universities Space Research Association, Huntsville, AL 35805, USA  Nicoló
Cibrario Dipartimento di Fisica, Universitá degli Studi di Torino, Via Pietro
Giuria 1, 10125 Torino, Italy  Stefano Ciprini Istituto Nazionale di Fisica
Nucleare, Sezione di Roma ”Tor Vergata”, Via della Ricerca Scientifica 1,
00133 Roma, Italy  Space Science Data Center, Agenzia Spaziale Italiana, Via
del Politecnico snc, 00133 Roma, Italy  Enrico Costa INAF Istituto di
Astrofisica e Planetologia Spaziali, Via del Fosso del Cavaliere 100, 00133
Roma, Italy  Alessandra De Rosa INAF Istituto di Astrofisica e Planetologia
Spaziali, Via del Fosso del Cavaliere 100, 00133 Roma, Italy  Ettore Del Monte
INAF Istituto di Astrofisica e Planetologia Spaziali, Via del Fosso del
Cavaliere 100, 00133 Roma, Italy  Laura Di Gesu ASI - Agenzia Spaziale
Italiana, Via del Politecnico snc, 00133 Roma, Italy  Immacolata Donnarumma
ASI - Agenzia Spaziale Italiana, Via del Politecnico snc, 00133 Roma, Italy
Victor Doroshenko Institut für Astronomie und Astrophysik, Universität
Tübingen, Sand 1, 72076 Tübingen, Germany  Michal Dovčiak Astronomical
Institute of the Czech Academy of Sciences, Boční II 1401/1, 14100 Praha 4,
Czech Republic  Steven R. Ehlert NASA Marshall Space Flight Center,
Huntsville, AL 35812, USA  Teruaki Enoto RIKEN Cluster for Pioneering
Research, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan  Yuri Evangelista INAF
Istituto di Astrofisica e Planetologia Spaziali, Via del Fosso del Cavaliere
100, 00133 Roma, Italy  Sergio Fabiani INAF Istituto di Astrofisica e
Planetologia Spaziali, Via del Fosso del Cavaliere 100, 00133 Roma, Italy
Javier A. Garcia California Institute of Technology, Pasadena, CA 91125, USA
Shuichi Gunji Yamagata University,1-4-12 Kojirakawa-machi, Yamagata-shi
990-8560, Japan  Kiyoshi Hayashida  Osaka University, 1-1 Yamadaoka, Suita,
Osaka 565-0871, Japan  Jeremy Heyl University of British Columbia, Vancouver,
BC V6T 1Z4, Canada  Wataru Iwakiri International Center for Hadron
Astrophysics, Chiba University, Chiba 263-8522, Japan  Ioannis Liodakis
Finnish Centre for Astronomy with ESO, 20014 University of Turku, Finland
Philip Kaaret NASA Marshall Space Flight Center, Huntsville, AL 35812, USA
Vladimir Karas Astronomical Institute of the Czech Academy of Sciences, Boční
II 1401/1, 14100 Praha 4, Czech Republic  Dawoon E. Kim INAF Istituto di
Astrofisica e Planetologia Spaziali, Via del Fosso del Cavaliere 100, I-00133
Roma, Italy  Dipartimento di Fisica, Università degli Studi di Roma “La
Sapienza”, Piazzale Aldo Moro 5, I-00185 Roma, Italy  Dipartimento di Fisica,
Università degli Studi di Roma “Tor Vergata”, Via della Ricerca Scientifica 1,
I-00133 Roma, Italy  Takao Kitaguchi  RIKEN Cluster for Pioneering Research,
2-1 Hirosawa, Wako, Saitama 351-0198, Japan  Jeffery J. Kolodziejczak NASA
Marshall Space Flight Center, Huntsville, AL 35812, USA  Henric Krawczynski
Physics Department and McDonnell Center for the Space Sciences, Washington
University in St. Louis, St. Louis, MO 63130, USA  Fabio La Monaca INAF
Istituto di Astrofisica e Planetologia Spaziali, Via del Fosso del Cavaliere
100, 00133 Roma, Italy  Luca Latronico Istituto Nazionale di Fisica Nucleare,
Sezione di Torino, Via Pietro Giuria 1, 10125 Torino, Italy  Grzegorz Madejski
Department of Physics and Kavli Institute for Particle Astrophysics and
Cosmology, Stanford University, Stanford, California 94305, USA  Simone
Maldera Istituto Nazionale di Fisica Nucleare, Sezione di Torino, Via Pietro
Giuria 1, 10125 Torino, Italy  Alberto Manfreda Istituto Nazionale di Fisica
Nucleare, Sezione di Napoli, Strada Comunale Cinthia, 80126 Napoli, Italy
Frédéric Marin Université de Strasbourg, CNRS, Observatoire Astronomique de
Strasbourg, UMR 7550, 67000 Strasbourg, France  Andrea Marinucci ASI - Agenzia
Spaziale Italiana, Via del Politecnico snc, 00133 Roma, Italy  Alan P.
Marscher Institute for Astrophysical Research, Boston University, 725
Commonwealth Avenue, Boston, MA 02215, USA  Herman L. Marshall MIT Kavli
Institute for Astrophysics and Space Research, Massachusetts Institute of
Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA  Francesco
Massaro Istituto Nazionale di Fisica Nucleare, Sezione di Torino, Via Pietro
Giuria 1, 10125 Torino, Italy  Dipartimento di Fisica, Università degli Studi
di Torino, Via Pietro Giuria 1, 10125 Torino, Italy  Giorgio Matt Dipartimento
di Matematica e Fisica, Università degli Studi Roma Tre, Via della Vasca
Navale 84, 00146 Roma, Italy  Riccardo Middei Space Science Data Center,
Agenzia Spaziale Italiana, Via del Politecnico snc, 00133 Roma, Italy  INAF
Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monte Porzio Catone
(RM), Italy  Ikuyuki Mitsuishi  Graduate School of Science, Division of
Particle and Astrophysical Science, Nagoya University, Furo-cho, Chikusa-ku,
Nagoya, Aichi 464-8602, Japan  Tsunefumi Mizuno Hiroshima Astrophysical
Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima,
Hiroshima 739-8526, Japan  Fabio Muleri INAF Istituto di Astrofisica e
Planetologia Spaziali, Via del Fosso del Cavaliere 100, 00133 Roma, Italy
Michela Negro Department of Physics and Astronomy, Louisiana State University,
Baton Rouge, LA 70803 USA  Stephen L. O’Dell NASA Marshall Space Flight
Center, Huntsville, AL 35812, USA  Chiara Oppedisano Istituto Nazionale di
Fisica Nucleare, Sezione di Torino, Via Pietro Giuria 1, 10125 Torino, Italy
Luigi Pacciani INAF Istituto di Astrofisica e Planetologia Spaziali, Via del
Fosso del Cavaliere 100, 00133 Roma, Italy  Alessandro Papitto INAF
Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monte Porzio Catone
(RM), Italy  George G. Pavlov Department of Astronomy and Astrophysics,
Pennsylvania State University, University Park, PA 16802, USA  Matteo Perri
Space Science Data Center, Agenzia Spaziale Italiana, Via del Politecnico snc,
00133 Roma, Italy  INAF Osservatorio Astronomico di Roma, Via Frascati 33,
00078 Monte Porzio Catone (RM), Italy  Melissa Pesce-Rollins Istituto
Nazionale di Fisica Nucleare, Sezione di Pisa, Largo B. Pontecorvo 3, 56127
Pisa, Italy  Pierre-Olivier Petrucci Université Grenoble Alpes, CNRS, IPAG,
38000 Grenoble, France  Andrea Possenti INAF Osservatorio Astronomico di
Cagliari, Via della Scienza 5, 09047 Selargius (CA), Italy  Juri Poutanen
Department of Physics and Astronomy, University of Turku, FI-20014, Finland
Simonetta Puccetti Space Science Data Center, Agenzia Spaziale Italiana, Via
del Politecnico snc, 00133 Roma, Italy  Brian D. Ramsey NASA Marshall Space
Flight Center, Huntsville, AL 35812, USA  John Rankin INAF Istituto di
Astrofisica e Planetologia Spaziali, Via del Fosso del Cavaliere 100, 00133
Roma, Italy  Ajay Ratheesh INAF Istituto di Astrofisica e Planetologia
Spaziali, Via del Fosso del Cavaliere 100, 00133 Roma, Italy  Oliver J.
Roberts Science and Technology Institute, Universities Space Research
Association, Huntsville, AL 35805, USA  Carmelo Sgró Istituto Nazionale di
Fisica Nucleare, Sezione di Pisa, Largo B. Pontecorvo 3, 56127 Pisa, Italy
Paolo Soffitta INAF Istituto di Astrofisica e Planetologia Spaziali, Via del
Fosso del Cavaliere 100, 00133 Roma, Italy  Gloria Spandre Istituto Nazionale
di Fisica Nucleare, Sezione di Pisa, Largo B. Pontecorvo 3, 56127 Pisa, Italy
Douglas A. Swartz Science and Technology Institute, Universities Space
Research Association, Huntsville, AL 35805, USA  Toru Tamagawa RIKEN Cluster
for Pioneering Research, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan  Fabrizio
Tavecchio INAF Osservatorio Astronomico di Brera, Via E. Bianchi 46, 23807
Merate (LC), Italy  Roberto Taverna Dipartimento di Fisica e Astronomia,
Università degli Studi di Padova, Via Marzolo 8, 35131 Padova, Italy  Yuzuru
Tawara  Graduate School of Science, Division of Particle and Astrophysical
Science, Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8602,
Japan  Allyn F. Tennant NASA Marshall Space Flight Center, Huntsville, AL
35812, USA  Nicholas E. Thomas NASA Marshall Space Flight Center, Huntsville,
AL 35812, USA  Francesco Tombesi Dipartimento di Fisica, Università degli
Studi di Roma ”Tor Vergata”, Via della Ricerca Scientifica 1, 00133 Roma,
Italy  Istituto Nazionale di Fisica Nucleare, Sezione di Roma ”Tor Vergata”,
Via della Ricerca Scientifica 1, 00133 Roma, Italy  Department of Astronomy,
University of Maryland, College Park, Maryland 20742, USA  Alessio Trois INAF
Osservatorio Astronomico di Cagliari, Via della Scienza 5, 09047 Selargius
(CA), Italy  Sergey Tsygankov Department of Physics and Astronomy, University
of Turku, FI-20014, Finland  Roberto Turolla Dipartimento di Fisica e
Astronomia, Università degli Studi di Padova, Via Marzolo 8, 35131 Padova,
Italy  Mullard Space Science Laboratory, University College London, Holmbury
St Mary, Dorking, Surrey RH5 6NT, UK  Jacco Vink Anton Pannekoek Institute for
Astronomy & GRAPPA, University of Amsterdam, Science Park 904, 1098 XH
Amsterdam, The Netherlands  Kinwah Wu Mullard Space Science Laboratory,
University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK
Silvia Zane Mullard Space Science Laboratory, University College London,
Holmbury St Mary, Dorking, Surrey RH5 6NT, UK

###### Abstract

We describe IXPE polarization observations of the Pulsar Wind Nebula (PWN) MSH
15$-$52, the ‘Cosmic Hand’. We find X-ray polarization across the PWN, with B
field vectors generally aligned with filamentary X-ray structures. High
significance polarization is seen in arcs surrounding the pulsar and toward
the end of the ‘jet’, with polarization degree $PD>70\%$, thus approaching the
maximum allowed synchrotron value. In contrast, the base of the jet has lower
polarization, indicating a complex magnetic field at significant angle to the
jet axis. We also detect significant polarization from PSR B1509$-$58 itself.
Although only the central pulse-phase bin of the pulse has high individual
significance, flanking bins provide lower significance detections and, in
conjunction with the X-ray image and radio polarization, can be used to
constrain rotating vector model solutions for the pulsar geometry.

particle acceleration, pulsars, polarization, radiation mechanisms: non-
thermal, pulsars: individual (PSR B1509-58)

††facilities: ATCA, CXO, IXPE

##  1 Introduction

PSR B1509$-$58 (=PSR J1513$-$5809) is a young ($\tau=1600$y), energetic
(${\dot{E}}=1.7\times 10^{37}{\rm erg\,s^{-1}}$), high field ($B_{s}=1.5\times
10^{13}$ G) pulsar embedded in the supernova remnant RCW89/G320.4$-$1.2/MSH
15$-$52 (Caswell et al., 1981). The relativistic particles and fields produced
by this pulsar power a bright X-ray pulsar wind nebula (PWN), whose
spectacular Chandra X-ray Observatory (CXO) image has earned the moniker ‘The
Cosmic Hand’ or ‘The Hand of God’. This structure and the surrounding
supernova remnant are detected from radio (Leung & Ng, 2016) to TeV (Aharonian
et al., 2005) energies with complex morphology, often complementary at
different energy bands. At a distance $d\approx 5$ kpc the $32^{\prime}$
diameter radio shell spans 47 pc. The PWN’s non-thermal X-ray emission extends
$\sim 8^{\prime}$ from the pulsar, making the PWN complex $\sim\,4\times$
larger in angle and $\sim 10\times$ larger in size than the famous Crab PWN.
MSH 15$-$52 shares a ‘torus+ jet’ morphology with the Crab, with a $\sim
10^{\prime\prime}$ sub-luminous X-ray zone around the pulsar representing the
pre-termination shock flow (Yatsu et al., 2009). The two bright X-ray arcs
wrapping the northern side of the pulsar may represent the distorted
equatorial torus of the shocked PWN or may represent field lines wrapped
around the termination shock by ram pressure or backflow in the surrounding
PWN (Figure 1). To the south along the torus axis is a prominent ridge of
X-ray emission extending at least $5^{\prime}$, often referred to as a ‘jet’
(Gaensler et al., 2002). To the northwest, non-thermal X-ray ridges form the
‘thumb’ and ‘fingers’ of the hand. The fingers extend to a region of softer
thermal X-ray emission to the north.

The $P_{s}=150$ ms pulsations are detected in the radio (Manchester et al.,
1982), X-ray (Seward & Harnden, 1982) and $\gamma$-ray (Abdo et al., 2010)
bands. The pulsed spectral energy distribution (SED) is actually quite soft
for a $\gamma$-ray pulsar, peaking at $\sim 10$ MeV, which may be associated
with its relatively large dipole field. As for most young gamma-ray pulsars,
the high energy emission lags the radio peak, here by $\Delta\phi\approx 0.3$.
At X-ray energies, the peak has two overlapping peaks, with separation
$\delta\phi\approx 0.2$ (de Plaa et al., 2003).

Existing polarization information on this system is limited. While the
supernova shell itself is quite bright in the radio, the non-thermal emission
is radio-faint. Radio observations with the Australia Telescope Compact Array
(ATCA) at 3 cm and 6 cm have detected significant linear polarization,
especially in the torus-like arcs (Leung & Ng, 2016, Zhang, et al. in prep.).
Here the inferred magnetic field follows the arcs as they wrap around the
pulsar. To the south, this polarized radio emission brackets the X-ray jet.
The X-ray jet fills a cavity in the radio emission, with little or no radio
flux apparent, as also noted by Gaensler et al. (2002). To the north, radio
emission seems to follow the thumb and finger structures but is rather faint
for reliable polarization maps. Like many young energetic pulsars, B1509-58
shows high linear polarization in the radio Crawford et al. (2001); Rookyard
et al. (2015a). From the ROSAT X-ray PWN structure Brazier & Becker (1997)
qualitatively estimated the viewing angle $i>70^{\circ}$, although a somewhat
smaller value is indicated by CXO data (Figure 3). There is a claim of a
pulsar phase averaged optical polarization of degree $PD\sim$10.4% by Wagner &
Seifert (2000), but the measurement is compromised by a bright field star, and
lacking any error bar or position angle (PA) estimate, this measurement needs
to be confirmed.

Here we report on the first measurements of X-ray polarization from this
complex, with robust detections in both the pulsar and the surrounding PWN,
and describe how these results constrain the system geometry.

Figure 1: Overview of the MSH 15$-$52 complex, the ‘Cosmic Hand’. The
background energy-coded CXO image ($0.5-1.2$ keV=red, $1.2-2$ keV=green, $2-7$
keV=blue) shows low energy (red/yellow) emission from thermal gas in
G320.4$-$1.2, and harder non-thermal emission from the PWN. Superposed are
IXPE 2-8 keV PCUBE-derived bars showing the polarization degree (PD, bar
length is for PD=1) and projected magnetic field direction, on a 30′′ grid.
Yellow bars show $>5\sigma$ polarization detections, thick white bars
$>3\sigma$ and thin white bars $>2\sigma$. In general the magnetic field
appears to follow the thumb, fingers and other linear structures (see Fig. 2
for region labels). Particularly strong polarization is associated with the
arc, especially to the northeast of the pulsar and with the far southeast end
of the jet. We have used a flux cut to trim anomalous vectors near the field-
of-view edge.

##  2 IXPE Observations of PSR B1509-58/MSH 15$-$52

The Imaging X-ray Polarimetry Explorer (IXPE), the first mission devoted to
spatially-resolved polarization measurements in X-rays (Weisskopf et al.,
2022), was successfully launched on December 9 2001. IXPE observed MSH 15$-$52
on 2-16 September 2022, 14-21 February 2023 and 13-19 March 2023 for a total
of $\sim 1.5$ Ms livetime. Data were extracted and analysed according to
standard procedures: HEASOFT 6.30.1 (Nasa High Energy Astrophysics Science
Archive Research Center (2014), Heasarc) was used to perform barycenter
corrections using the DE421 JPL ephemeris. ixpeobssim V30.2.2 (Baldini et al.,
2022) was used to do energy calibration, detector WCS correction, bad aspect-
ratio corrections, and all further analysis, including phase folding at the
pulsar ephemeris.

Background events were cleaned from the data following the procedure of Di
Marco et al. (2023). The residual instrumental background was modeled from
1.5Ms of cleaned IXPE source-free exposure in the fields of several high-
latitude sources (MGC$-$5$-$23$-$16, 1ES 0299+200, PG 1553, PSR B0540$-$69 and
IC 4329A). MSH 15$-$52 lies close to the Galactic ridge so some contribution
from background X-rays is expected as well. However, it covers most of the
IXPE field of view so we cannot extract a local background spectrum directly
from IXPE. Instead we use CXO observations to compute the background flux,
passed through the IXPE instrument response, south of the thumb, finding a
count rate $\sim 1.1\times$ the instrumental background, and so we increase
the background spectrum surface brightness by this factor. This unpolarized
background surface brightness ($8.9\times 10^{-8}$cnts/arcsec2/s/det, 2-5.5
keV; $1.07\times 10^{-7}$cnts/arcsec2/s/det, 2-8 keV) is scaled and subtracted
from the flux of each aperture.

Since DeLaney et al. (2006) have noted temporal variations in the fine
structure of the PWN, especially in knots near the pulsar, but also in the jet
feature, we collected a contemporaneous 28 ks CXO observation of the PWN
(ObsIDs 23540, 27448) to have a current high resolution image for comparison.

Figure 1 gives an overview of the IXPE polarization measurements superimposed
on an energy-coded image from archival CXO exposures (ObsIDs 0754, 3833, 5534,
5535, 6116, 6117 – 204 ks livetime total). Here we show the projected magnetic
field direction (orthogonal to the Electric Vector Position Angle, EVPA)
measured on a $30^{\prime\prime}$ grid, comparable to the resolution of the
IXPE PSF. Complex polarization features extend throughout the nebula.

Figure 2: Morphological regions of interest on the IXPE combined DU1-3 2-8 keV
image. The Hand’s ‘wrist’ includes the bright jet and flanking filament
structures, divided here into several regions.

We can use the deep archival data to define nebula regions of interest (Figure
2 and Table 1). This allows us to discuss the polarization properties of
extended regions too faint for high S/N mapping. The 2022 CXO image does show
small departures from the archival morphology; most are changes in the shock
structure near the pulsar, within the central IXPE resolution element,
although there are also small changes in the jet $1-2^{\prime}$ from the
pulsar. None affect the locations of our extended regions. In the region near
the pulsar, CXO-maps show an inner arc, not resolved by IXPE (Figure 3).
Strong, high significance polarization follows the outer arc, as also seen in
the 6cm radio maps. The arc magnetic field structure extends, with field lines
parallel to the jet, in left and right ‘arc extensions’. It is more faintly
visible in the sheath regions flanking the X-ray bright jet. The field also
clearly follows the curve of the thumb region. All of these features are also
discernible in the radio. In addition we measure fields paralleling the
‘finger’ structures – in the radio, these are lost to bright emission from
shock in the SNR shell. The shell interaction produces the low energy X-ray
emission appearing in red and yellow to the north in Figure 1. This thermal
emission, near the edge of the IXPE field-of-view, is unpolarized.

The general pattern of polarization in Fig. 1 is as expected, with the
magnetic field lines following the filamentary nebula structure. The highest
fractional polarizations, in the outer arc, thumb, and end of the jet, reach
$PD\sim 70\%$ (after background subtraction). By integrating over the regions
of Figure 2, we also see that the magnetic field is aligned with the thin
‘finger’ structures, but with substantial background from the thermal emission
at the fingertips, we suspect that the $PD$ in these regions is
underestimated.

Table 1: CXO spectral measurements and IXPE PCUBE Polarizations for the regions of Figure 2. Region | Flux | $\Gamma_{X}$ | $B_{Eq}$ | Q | U | Q, U err | PD | $\psi$ | Sig
---|---|---|---|---|---|---|---|---|---
|  $0.5-8\,$keV |  |  $\mu$G |  |  |  |  | $\circ$ |
Fingers | 26.98 | $2.35\pm 0.01$ | 21 | $-0.119$ | $0.114$ | 0.022 | $0.165\pm 0.022$ | $68.0\pm 3.8$ | 7.5
Index | 5.64 | $2.22\pm 0.01$ | 24 | $-0.216$ | $0.065$ | 0.044 | $0.225\pm 0.044$ | $81.7\pm 5.6$ | 5.1
Middle | 6.65 | $2.34\pm 0.01$ | 29 | $-0.055$ | $0.183$ | 0.043 | $0.191\pm 0.043$ | $53.3\pm 6.5$ | 4.4
Ring | 4.95 | $2.26\pm 0.02$ | 22 | $-0.177$ | $0.121$ | 0.050 | $0.214\pm 0.050$ | $72.8\pm 6.7$ | 4.3
Thumb† | 6.15 | $1.92\pm 0.01$ | 16 | $0.313$ | $-0.085$ | 0.038 | $0.324\pm 0.038$ | $-7.6\pm 3.4$ | 8.4
Outer Arc† | 6.38 | $1.87\pm 0.01$ | 27 | $0.348$ | $-0.020$ | 0.026 | $0.348\pm 0.026$ | $-1.7\pm 2.1$ | 13.5
ArcExt L | 4.10 | $1.93\pm 0.01$ | 22 | $0.084$ | $0.412$ | 0.043 | $0.421\pm 0.043$ | $39.2\pm 2.9$ | 9.7
ArcExt R | 9.12 | $1.84\pm 0.01$ | 16 | $-0.109$ | $0.075$ | 0.027 | $0.133\pm 0.027$ | $72.7\pm 5.9$ | 4.9
Jet | 15.81 | $1.69\pm 0.01$ | 18 | $0.238$ | $0.089$ | 0.021 | $0.254\pm 0.021$ | $10.2\pm 2.3$ | 12.2
Sheath L | $3.47$ | $1.77\pm 0.02$ | 16 | $-0.001$ | $0.142$ | 0.034 | $0.142\pm 0.034$ | $45.2\pm 6.9$ | 4.2
Sheath R | 9.68 | $1.80\pm 0.01$ | 20 | $0.092$ | $0.256$ | 0.025 | $0.272\pm 0.025$ | $35.1\pm 2.6$ | 11.0
Jet-Sheath∗ |  |  |  | $0.602$ | $-0.146$ | 0.075 | $0.620\pm 0.075$ | $-6.8\pm 3.4$ | 8.3
J1 | 4.58 | $1.64\pm 0.01$ | 25 | $0.012$ | $0.173$ | 0.037 | $0.173\pm 0.037$ | $43.1\pm 6.1$ | 4.7
J2 | 6.54 | $1.74\pm 0.01$ | 19 | $0.276$ | $0.034$ | 0.035 | $0.278\pm 0.035$ | $3.5\pm 3.6$ | 8.0
J3 | 4.59 | $1.67\pm 0.01$ | 16 | $0.466$ | $0.053$ | 0.041 | $0.469\pm 0.041$ | $3.2\pm 2.5$ | 11.5
J1-Sh1∗ |  |  |  | $0.281$ | $0.123$ | 0.109 | $0.307\pm 0.109$ | $11.8\pm 10.1$ | 2.8
J2-Sh2∗ |  |  |  | $0.526$ | $-0.293$ | 0.090 | $0.602\pm 0.090$ | $-14.6\pm 4.3$ | 6.7
J3-Sh3∗ |  |  |  | $0.819$ | $-0.145$ | 0.159 | $0.831\pm 0.159$ | $-5.0\pm 5.3$ | 5.2

The most unusual feature is the X-ray bright hard spectrum ‘jet’, which is
essentially invisible in the radio, implying a low-energy cutoff in the jet
electron spectrum. This may be an intrinsic cut-off in the injected electron
spectrum or the result of limited time available for cooling in the rapid jet
flow. We also note that the overall polarization level is low at the base of
the jet region. Interestingly, the weak polarization that we do see appears to
be at substantial angle to that of the bracketing nebula. These are, of
course, 3-D structures so it seems likely that the jet zone is viewed through
a plasma emitting like the ‘sheath’ zones to either side. If we subtract the
average Stokes I, Q, and U of this sheath (scaled for area) we do indeed see
jet polarization increase to $PD>60\%$, with the measured EVPA implying an
average magnetic field at angles up to 50∘ from the jet axis.

Figure 3: Inner Region of MSH 15$-$52 from $0.5-8$ keV CXO data, centered on
the pulsar and smoothed with a $1^{\prime\prime}$ Gaussian kernel. Note the
sub-luminous zone (marked by an ellipse 2/3 of the zone size), the symmetry
axis (red arrow), the bright polar outflow to the northwest and the fainter
more diffuse outflow to the southeast. The magenta square marks the central
pixel of the PCUBE polarization map, Fig. 1 (approximating IXPE’s resolution);
‘simultaneous fit’ pixels are half this size, Fig. 4. The ‘Inner Arc’ passes
through the northwest polar outflow, while the ‘Outer Arc’ is partly visible
to the north.

##  3 Phase-Resolved Analysis of PSR B1509-58

We obtained contemporaneous Parkes L-band radio observations (2023 2/7, 2/20,
2/26 and 3/1), and folded them with the same ephemeris used for the IXPE X-ray
events to confirm that the $\Delta\phi=0.25$ radio-X-ray phase lag (Kawai et
al., 1991; Cusumano et al., 2001) remains valid.

To extract the nebula map and the pulsar polarization we employ the
‘Simultaneous Fitting’ technique of Wong et al. (2023). This uses the
contemporaneous 2022 CXO image of the nebula, with the point source
subtracted, to define the intensity (and local spectrum) of the extended
emission at the IXPE observation epoch. For the pulsar point source
contribution, we also rely on CXO data, using the ACIS-CC and HRC analysis of
Hu et al. (2017) to define the light curve and phase-varying spectral index of
the pulsar emission. The phase dependent pulsar and spatially dependent nebula
spectra are folded through the IXPE response, using ixpeobsim, to predict the
IXPE counts as a function of position and phase. Note that PSR B1509$-$58 is
relatively bright at minimum, at $\sim 4\%$ of its peak flux – this means that
the phase invariant DC emission contributes $\sim 11\%$ of the pulsed flux. To
model faint PWN regions the uniform background must be included in the
simultaneous fitting model; here we use the instrumental background, as local
photon background is included in the CXO-derived flux.

Simultaneous fitting defines a set of spatial and phase bins, and uses the
predicted IXPE counts from the nebula, background, and PSF-spread pulsar to
define the expected PSR/PWN contributions to each bin. It then executes a
global least-squares fit for the pulsar polarization at each phase and the
phase-independent nebular polarization at each spatial pixel. Here we define a
$13\times 11$, $15^{\prime\prime}$ pixel grid centered on the pulsar. We use
$2-5.5$ keV photons, to best isolate the pulsar polarization signal (which is
slightly softer than the spectrum of the inner PWN) and simple ‘Moments
Ellipticity’ weights to quantify the accuracy of the polarization
reconstruction of each event. Note that with different PSFs for each detector
(as measured from ground calibration images) and different spacecraft
orientations for each of the three IXPE pointings, we have nine measurements
of the combined PSR/PWN polarization signal in each spatial and phase bin, all
of which must be simultaneously fit.

Because errors in the reconstructed photon conversion point are correlated
with the reconstructed polarization vector, bright point sources (and any
sharp flux gradient) will have a ‘halo’ of polarization at scales less than
the PSF FWHM, which can be corrected by an iterative estimate of this so-
called polarization leakage (Bucciantini et al., 2023). Here we apply the
energy-dependent version of this correction, using the detailed ground-
measured PSFs of the three telescope assemblies, as outlined in Wong et al.
(2023). This correction is applied before the simultaneous fit extraction of
the component polarizations. The correction makes modest ($<20\%$) amendments
to the polarization degree in the inner few arcmin, especially associated with
the relatively sharp arcs to the north of the pulsar. Also, when the spatial
bins are smaller than the PSF FWHM and the counts/bin are low, anti-correlated
fluctuations between adjacent pixels increase the scatter and error in fit $q$
and $u$. Here we mitigated this by smoothing the $q$ and $u$ maps by the PSFs
for a decrease in fluctuations, at a cost of some spatial resolution.

Figure 4: 2-5.5 keV polarization bars (orthogonal to the EVPA), as determined
from simultaneous fitting inside the magenta region (cyan bars $>4\sigma$
significance, green bars $>3\sigma$, white bars $>2\sigma$; the labels at
upper right show the bar length for PD=1). Here we employ a
$15^{\prime\prime}$ grid to better resolve the pulsar from inner nebula
structures. Background is the full-band CXO image.

In Figure 5 we show the simultaneous fitting-derived pulsar X-ray EVPA
estimates along with radio polarization measurements and the IXPE X-ray light
curve for reference. Only one X-ray phase bin, near the center of the peak, is
significant with a $PD$ of 17.5% at $3.7\sigma$. The large pulse-minimum bin
formally has a very high $PD\sim 1$ at low significance. However the PCUBE
analysis shows small total polarization in the central pixel – simultaneous
fitting evidently optimizes the central region fit by introducing some $q$ and
$u$ into the faint minimum phase pulsar component and producing canceling PSR
minimum and PWN polarizations in this phase bin. The other bins in the X-ray
peak have low $PD=10-20\%$; at $\sim 2-3\sigma$ significance/bin there is no
definitive polarization detection, although the EVPA values do assume an
intriguing smooth sweep across the X-ray peak.

##  4 Discussion

The background-subtracted, leakage-corrected polarization map (Fig. 1) has
several $>5\sigma$ polarization regions. The most significant (Left Arc
extension) pixel has a background-subtracted $PD=0.72\pm 0.08$. A few low-
count pixels near the nebula edge have higher PD, the most extreme being in
the Left Arc Extension, with $PD=0.87\pm 0.14$. Thus all pixels are consistent
with $PD<0.75$ at the $1\sigma$ level. Other highly polarized pixels are at
the jet end ($PD=0.65\pm 0.12$), the thumb base ($PD=0.66\pm 0.11$) and the
index finger ($PD=0.73\pm 0.20$). The jet as a region is highly polarized
toward its end with $PD=0.83\pm 0.16$ at its far (J3) end, if one subtracts
the adjacent sheath emission as a background. Thus, as also seen in the Vela
PWN (Xie et al., 2022), polarization approaches
$PD=\Gamma_{X}/(\Gamma_{X}+2/3)$, the maximum allowed for synchrotron
polarization at the observed X-ray photon index $\Gamma_{X}$ in a uniform
magnetic field. For example, in the J3-Sh3 region, the maximum allowed value
is $PD=0.72$; the observed polarization is $0.7\sigma$ above this value,
consistent with a statistical fluctuation.

In the inner region, simultaneous fitting lets us map fields closer to the
pulsar (Fig. 4). Here again the strongest polarization follows the outer arc
and the left arc extension, with a peak value of $PD=0.64\pm 0.14$. There is
also a $PD=0.54\pm 0.19$ polarized pixel of modest $2.8\sigma$ significance
located on the ridge of the inner arc. This should be a pure nebula
measurement, as it comes from the nebula portion of the simultaneous fit,
generated with the contemporaneous CXO-defined structure, and has also been
corrected for polarization leakage. However, at only $15^{\prime\prime}$ from
the pulsar PSF peak, some concern about systematic effects persists. At both
scales, polarization at the base of the jet is low.

Table 1 lists the average polarization degree and angle in the larger regions
defined in Figure 2. We can estimate the regions’ magnetic field strength
under the assumption of equipartition. For an optically thin region filled
with relativistic electrons and magnetic field emitting synchrotron radiation
the equipartition field is

$B_{\rm
Eq}=46\left[\frac{J_{\rm-20}(E_{1},E_{2})\sigma}{\phi}\frac{C_{1.5-\Gamma}(E_{m},E_{M})}{C_{2-\Gamma}(E_{1},E_{2})}\right]^{2/7}\mu
G$ (1)

where

$C_{q}(x_{1},x_{2})=\frac{x_{2}^{q}-x_{1}^{q}}{q}.$ (2)

$J_{\rm-20}(E_{1},E_{2})=4\pi f_{X}(E_{1},E_{2})d^{2}/V$ is the observed
emissivity (in $10^{-20}$ erg s-1 cm-3, between $E_{1}$ keV and $E_{2}$ keV),
$\sigma_{B}=w_{B}/w_{e}$ is the magnetization parameter, $\phi$ the filling
factor, and $E_{m}$ and $E_{M}$ the minimum and maximum energies, in keV, of
the synchrotron spectrum with photon index $\Gamma$. We assume that the
structures are cylindrical, with diameter set to the observed region width. We
list the derived equipartition fields in Table 1 for $\sigma=\phi=1$,
$E_{m}=0.01\,\rm{keV}$ and $E_{M}=10\,\rm{keV}$.

MSH 15$-$52 is complex but a few trends can be extracted from Table 1. First,
the fingers region is notably softer than the bulk of the PWN. This may, in
part be due to contamination by the soft thermal emission to the north.
However, the ‘Thumb’ with $\Gamma=1.92$ is free of the thermal emission but
still somewhat softer than the outer Arc. The hardest feature is, of course,
the ‘jet’ as seen in the color image (Fig. 1). This suggests that this feature
contains the freshest electron population and that the outer features have
suffered some synchrotron burn-off. Indeed the jet may represent a site of
$e^{\pm}$ re-acceleration and the low polarization at its base may be, in
part, due to magnetic turbulence and dissipation there. The spectral trends
are broadly consistent with those found by An et al. (2014). These authors,
using NuSTAR, infer a nebula-averaged spectral break at $\sim 6$ keV. Thus,
the average CXO spectral indices shown here should not resolve a full
$\Delta\Gamma=0.5$ cooling break.

Our equipartition field estimates are of course subject to the uncertain 3-D
geometry and $\phi$ fill factor. There does seem to be a trend of higher
fields to the north, which may be associated with compression from interaction
with G320.4$-$1.2. We also note that the equipartition field strength appears
to decrease along the jet, although the field becomes more uniform, as shown
by the $PD$ increase as one moves away from the pulsar.

The field orientations for the morphologial regions support the pattern in
Fig. 1, with the arc and thumb fields well aligned with the curved ridges. The
mean ‘Jet’ field is oriented $\sim 25-35^{\circ}$ from the surrounding
‘Sheath’ regions. If we imagine that these are 3D structures, with the Sheath
surrounding the jet, we can subtract the mean sheath flux, to find that the
offset angle increases to $\sim 40-50^{\circ}$ and the residual polarization
is quite high at $PD=62\pm 8$%. The jet field is not fully transverse to the
jet axis, but the significant orthogonal component might implicate a helical
structure. We subdivided the jet into three regions, finding a strongly
increasing polarization as one moves downstream. The B orientations do not
show a smooth trend, even after subtracting flanking sheath fluxes. The
brightest mid jet region, however has the largest angle to the local jet axis
at $\approx 50^{\circ}$.

Figure 5: IXPE Top: phase bin PA estimates (green; $q$ and $u$ values
available as data behind figure). The blue-circled bin at pulse center is
individually significant ($3.7\sigma$); cyan-circled points are at
$2-3\sigma$. The radio PAs (red) can be fit (excluding the circled point) with
the normal mode (magenta) and orthogonal mode (red) $i\sim 60^{\circ}$ models.
For the $i\sim 120^{\circ}$ model only the orthogonal mode solution is
tenable. Note that the $i\sim 60^{\circ}$ models better accommodate the late
radio point and the X-ray data. Somewhat better match can be found for the
orthogonal mode model when the X-ray emission is at higher altitude, with
modest altitude $h$ (dashed curve). Bottom: IXPE pulsar light curve (blue),
Parked 1.4GHz light curve (red).

Pulsar polarization can also be related to the PWN geometry. Examining the
CXO-measured fine structure in the inner nebula, we see a general symmetry
axis at $\psi=140\pm 5^{\circ}$. PWN features are best described by tangential
views of structure in the MHD flow, as described by Komissarov & Lyubarsky
(2003), but the geometrical inferences from a torus-jet picture with
cylindrical symmetry are robust. As for the Crab, there is a sub-luminous zone
surrounding the pulsar, first described for MSH 15$-$52 by Yatsu et al.
(2009), which marks the equatorial flow prior to the termination shock (marked
in Fig. 3 by an ellipse). For the Crab this zone is bracketed by the inner
ring and wisps. The Crab wisps are brighter to the northwest, and if
interpreted as due to Doppler boosting in mildly relativistic post-shock flow
(Ng & Romani, 2008) this determines the 3-D orientation of the spin axis. For
MSH 15$-$52 the zone has no bright edge and the Doppler boosting is not
obvious. So while the ellipticity of the zone constrains the spin axis
inclination to the Earth line-of-sight, both $i=60\pm 2^{\circ}$ with the
southeast axis out of the plane of the sky (since the ‘jet’ to the southeast
would then approach us, we call this the ‘Jet’ solution) and $i=120\pm
2^{\circ}$ with that axis into the plane of the sky (the ‘C-Jet’ solution) are
viable. One might interpret the blob of polar emission to the northwest in
Fig. 3 as a Doppler boosted ‘jet’. However, it is diffuse and is more likely
outflow analogous to the dome of PWN emission northwest of the Crab, rather
than a collimated relativistic jet flow.

We can compare this geometry with that inferred from radio pulsar polarization
measurements. Figure 5 shows Parkes 1.4GHz EVPA values, referenced to infinite
frequency for a rotation measure of $RM=216.0\,{\rm rad/m^{2}}$ (Johnston &
Weisberg, 2006), where phase bins with linear polarization detected with
$>2.5\sigma$ significance are plotted. Traditionally one fits the EVPA data
$\psi(\phi)$ to the rotating vector model (RVM, Radhakrishnan & Cooke, 1969),
which can be generalized to include the effect of Doppler boosting of the
rotating emission point at height $h=r/R_{\rm LC}$ as (e.g., Poutanen, 2020)

${\rm tan}(\psi-\psi_{0})=\\\ \qquad\frac{{\rm sin}\theta\,{\rm
sin}(\phi-\phi_{0})+h[{\rm sin}i\,{\rm sin}\theta+{\rm cos}i\,{\rm
cos}\theta\,{\rm cos}(\phi-\phi_{0})]}{{\rm cos}i\,{\rm sin}\theta\,{\rm
cos}(\phi-\phi_{0})-{\rm sin}i\,{\rm cos}\theta-h{\rm cos}\theta\,{\rm
sin}(\phi-\phi_{0})}$

where $h\approx 0$ for the low altitude radio emission. Here $i$ is the
inclination of the spin axis to the line of sight, $\theta$ is the angle
between the magnetic and spin axes and the magnetic axis passes closest to the
line of sight at $\phi=\phi_{0}$ with impact parameter $\beta=i-\theta$ and
EVPA $\psi_{0}$. Note that the sign of the denominator addresses the ‘$\psi$
convention problem’ (Everett & Weisberg, 2001). With the limited radio phase
coverage, a simple $h=0$ fit to the radio data is not particularly
constraining (Rookyard et al., 2015b), but if we impose the prior constraints
on $\psi_{0}$ (which is either along or orthogonal to the projected spin axis
at phase $\phi_{0}$, Johnston et al., 2005) and $i$ (two options, above) from
the X-ray image, we obtain fits with well constrained parameters and small
covariance. In Table 2 we show the Markov Chain Monte Carlo fit parameters for
the three viable options (the normal mode orientation with $i\approx
120^{\circ}$ provides no acceptable fit to the radio data). Both orthogonal
mode solutions provide very good fits. The one normal mode solution is worse,
but with a p-value of $0.025$, still acceptable. If the X-ray image
constraints are relaxed, the best-fit solutions remain stable, although errors
of course increase and there is substantial $\phi_{0}-\psi_{0}$ covariance.

Table 2: MCMC Radio RVM fits – see Figure 5. Parameter | Norm Jet | Orth Jet | Orth C-Jet
---|---|---|---
$i(^{\circ})$ | $58.0\pm 1.9$ | $60.0\pm 1.9$ | $119.7\pm 2.0$
$\theta(^{\circ})$ | $87.6\pm 3.8$ | $123.0\pm 2.7$ | $148.5\pm 1.6$
$\phi_{0}(^{\circ})$ | $48.9\pm 3.3$ | $-22.6\pm 3.0$ | $-18.1\pm 2.4$
$\Psi_{0}(^{\circ})$ | $138.1\pm 2.7$ | $49.9\pm 2.9$ | $52.8\pm 2.7$
$\beta(^{\circ})$ | -28.7 | -63.0 | -28.8
$\chi^{2}$/63 DoF | 1.37 | 1.13 | 1.19

The IXPE polarization data (and the late phase radio point) can help us
distinguish between RVM models. The $i\approx 120^{\circ}$ RVM model cannot
explain these points as the model EVPA is far off. The $i\approx 60^{\circ}$
models have EVPA increasing past the radio peak, and so more plausibly account
for these data. In fact, for these models the post radio EVPA increases
slightly for higher altitude emission; the orthogonal $i\approx 60^{\circ}$
model can match the IXPE EVPA and approach the late phase radio point if their
emission is from higher altitude with $h>0$. A fit to the X-ray data formally
gives $h\approx 0.15\pm 0.05$, but there are multiple minima and large
departures at late phases. Note that the $i\approx 60^{\circ}$ orthogonal
model sweep is slowest near the significant IXPE detection. For this case some
loss of polarization signal might be attributed to sweep in the surrounding
bins. Non-zero $h$ does not help the normal mode model, as it already has EVPA
larger than that of the late phase points.

The $i\approx 60^{\circ}$ orthogonal RVM model has the minimum $\chi^{2}$ but
there is a major peculiarity: the radio peak appears when the associated
magnetic pole sweeps $|\beta|=|i-\theta|=63^{\circ}$ from the Earth line-of-
sight, while the opposite pole, sweeping $3^{\circ}$ away at $\phi=0.44$ shows
no radio emission. In contrast the other two models have large, but less
extreme $\beta\approx-29^{\circ}$. Of these the Normal mode model has the
radio pulse leading the radio axis by a substantial $49^{\circ}$, while the
orthogonal solution would have the radio pulse trailing the magnetic axis.
Thus, no solution is ideal and all require a very large, partly filled radio
beam. Additional significant X-ray EVPAs would certainly help the model
discrimination, as would more late-phase radio measurements.

##  5 Conclusions

In sum, the CXO-measured X-ray morphology of the inner PWN does constrain the
3-D spin axis and helps select between otherwise viable RVM fits to the PSR
B1509$-$58 radio polarization. With IXPE we also extract a single phase bin of
pulsar X-ray polarization. This is plausibly interpreted as an extension of
the radio polarization sweep, but with only one significant bin, it its
difficult to make detailed model tests. Further IXPE observations could
promote 2-3 more bins to 3$\sigma$ significance, but would probably require
$\sim 2$ Ms of additional exposure.

We conclude by noting that the rich polarization structure of the MSH 15$-$52
PWN reflects the interplay of axisymetric pulsar outflow and complex, possibly
unstable interaction with the surrounding SNR. Although, unlike the Crab and
Vela PWNe, toroidal symmetry does not dominate the polarization pattern, the
polarization degree of MSH 15$-$52 is similar to that found earlier by IXPE
for the Crab and Vela: polarization is very high in parts of the hard spectrum
emission regions, approaching the maximum $PD$ allowed for synchrotron
emission (Xie et al., 2022). This suggests that these portions of the PWN
contain uniform fields with little turbulence. On the other hand, the base of
the jet which may be re-accelerating particles has a low polarization and
complex field geometry. It seems that if diffusive shock acceleration (DSA)
energizes the PWN particles, then much of the radiation comes from uniform
field zones separate from the acceleration sites. Alternatively a lower-
turbulence mechanism, possibly associated with magnetic reconnection, may be
involved. Full mapping of the field geometry requires higher resolution and
sensitivity than IXPE can provide. But even the present data provide a
visually striking polarization map of the Cosmic Hand’s fields and some
important challenges to MHD PWN models.

## acknowledgments

The Imaging X-ray Polarimetry Explorer (IXPE) is a joint US and Italian
mission. The US contribution is supported by the National Aeronautics and
Space Administration (NASA) and led and managed by its Marshall Space Flight
Center (MSFC), with industry partner Ball Aerospace (contract NNM15AA18C). The
Italian contribution is supported by the Italian Space Agency (Agenzia
Spaziale Italiana, ASI) through contract ASI-OHBI-2017-12-I.0, agreements ASI-
INAF-2017-12-H0 and ASI-INFN-2017.13-H0, and its Space Science Data Center
(SSDC) with agreements ASI-INAF-2022-14-HH.0 and ASI-INFN 2021-43-HH.0, and by
the Istituto Nazionale di Astrofisica (INAF) and the Istituto Nazionale di
Fisica Nucleare (INFN) in Italy. This research used data products provided by
the IXPE Team (MSFC, SSDC, INAF, and INFN) and distributed with additional
software tools by the High-Energy Astrophysics Science Archive Research Center
(HEASARC), at NASA Goddard Space Flight Center (GSFC). Funding for this work
was provided in part by contract NNM17AA26C from the MSFC to Stanford and
80MSFC17C0012 to MIT in support of the IXPE project. Support for this work was
provided in part by the NASA through the Smithsonian Astrophysical Observatory
(SAO)contract SV3-73016 to MIT for support of the Chandra X-Ray Center (CXC),
which is operated by SAO for and on behalf of NASA under contract NAS8-03060.
C.-Y. Ng and Y.-J. Yang are supported by a GRF grant of the Hong Kong
Government under HKU 17305419. N.B. was supported by the INAF MiniGrant
“PWNnumpol - Numerical Studies of Pulsar Wind Nebulae in The Light of IXPE”.

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 714, 927, doi: 10.1088/0004-637X/714/1/927
  * Aharonian et al. (2005) Aharonian, F., Akhperjanian, A. G., Aye, K. M., et al. 2005, A&A, 435, L17, doi: 10.1051/0004-6361:200500105
  * An et al. (2014) An, H., Madsen, K. K., Reynolds, S. P., et al. 2014, ApJ, 793, 90, doi: 10.1088/0004-637X/793/2/90
  * Baldini et al. (2022) Baldini, L., Bucciantini, N., Lalla, N. D., et al. 2022, SoftwareX, 19, 101194, doi: 10.1016/j.softx.2022.101194
  * Brazier & Becker (1997) Brazier, K. T. S., & Becker, W. 1997, MNRAS, 284, 335, doi: 10.1093/mnras/284.2.335
  * Bucciantini et al. (2023) Bucciantini, N., Di Lalla, N., Romani, R. W. R., et al. 2023, A&A, 672, A66, doi: 10.1051/0004-6361/202245744
  * Caswell et al. (1981) Caswell, J. L., Milne, D. K., & Wellington, K. J. 1981, MNRAS, 195, 89, doi: 10.1093/mnras/195.1.89
  * Crawford et al. (2001) Crawford, F., Manchester, R. N., & Kaspi, V. M. 2001, AJ, 122, 2001, doi: 10.1086/323107
  * Cusumano et al. (2001) Cusumano, G., Mineo, T., Massaro, E., et al. 2001, A&A, 375, 397, doi: 10.1051/0004-6361:20010884
  * de Plaa et al. (2003) de Plaa, J., Kuiper, L., & Hermsen, W. 2003, A&A, 400, 1013, doi: 10.1051/0004-6361:20030039
  * DeLaney et al. (2006) DeLaney, T., Gaensler, B. M., Arons, J., & Pivovaroff, M. J. 2006, ApJ, 640, 929, doi: 10.1086/500189
  * Di Marco et al. (2023) Di Marco, A., Soffitta, P., Costa, E., et al. 2023, Handling Background in IXPE Polarimetric Data, arXiv, doi: 10.48550/arXiv.2302.02927
  * Everett & Weisberg (2001) Everett, J. E., & Weisberg, J. M. 2001, ApJ, 553, 341, doi: 10.1086/320652
  * Gaensler et al. (2002) Gaensler, B. M., Arons, J., Kaspi, V. M., et al. 2002, ApJ, 569, 878, doi: 10.1086/339354
  * Hu et al. (2017) Hu, C.-P., Ng, C. Y., Takata, J., Shannon, R. M., & Johnston, S. 2017, ApJ, 838, 156, doi: 10.3847/1538-4357/aa67e3
  * Johnston et al. (2005) Johnston, S., Hobbs, G., Vigeland, S., et al. 2005, MNRAS, 364, 1397, doi: 10.1111/j.1365-2966.2005.09669.x
  * Johnston & Weisberg (2006) Johnston, S., & Weisberg, J. M. 2006, MNRAS, 368, 1856, doi: 10.1111/j.1365-2966.2006.10263.x
  * Kawai et al. (1991) Kawai, N., Okayasu, R., Brinkmann, W., et al. 1991, ApJ, 383, L65, doi: 10.1086/186242
  * Komissarov & Lyubarsky (2003) Komissarov, S. S., & Lyubarsky, Y. E. 2003, MNRAS, 344, L93, doi: 10.1046/j.1365-8711.2003.07097.x
  * Leung & Ng (2016) Leung, W. Y., & Ng, C. Y. 2016, in Supernova Remnants: An Odyssey in Space after Stellar Death, 53
  * Manchester et al. (1982) Manchester, R. N., Tuohy, I. R., & Damico, N. 1982, ApJ, 262, L31, doi: 10.1086/183906
  * Nasa High Energy Astrophysics Science Archive Research Center (2014) (Heasarc) Nasa High Energy Astrophysics Science Archive Research Center (Heasarc). 2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, record ascl:1408.004.  http://ascl.net/1408.004
  * Ng & Romani (2008) Ng, C. Y., & Romani, R. W. 2008, ApJ, 673, 411, doi: 10.1086/523935
  * Poutanen (2020) Poutanen, J. 2020, A&A, 641, A166, doi: 10.1051/0004-6361/202038689
  * Radhakrishnan & Cooke (1969) Radhakrishnan, V., & Cooke, D. J. 1969, Astrophys. Lett., 3, 225
  * Rookyard et al. (2015a) Rookyard, S. C., Weltevrede, P., & Johnston, S. 2015a, MNRAS, 446, 3356, doi: 10.1093/mnras/stu2083
  * Rookyard et al. (2015b) —. 2015b, MNRAS, 446, 3367, doi: 10.1093/mnras/stu2236
  * Seward & Harnden (1982) Seward, F. D., & Harnden, F. R., J. 1982, ApJ, 256, L45, doi: 10.1086/183793
  * Wagner & Seifert (2000) Wagner, S. J., & Seifert, W. 2000, in Astronomical Society of the Pacific Conference Series, Vol. 202, IAU Colloq. 177: Pulsar Astronomy - 2000 and Beyond, ed. M. Kramer, N. Wex, & R. Wielebinski, 315
  * Weisskopf et al. (2022) Weisskopf, M. C., Soffitta, P., Baldini, L., et al. 2022, Journal of Astronomical Telescopes, Instruments, and Systems, 8, 026002, doi: 10.1117/1.JATIS.8.2.026002
  * Wong et al. (2023) Wong, J., Romani, R. W., & Dinsmore, J. T. 2023, ApJ, 953, 28, doi: 10.48550/arXiv.2306.08788
  * Xie et al. (2022) Xie, F., Di Marco, A., La Monaca, F., et al. 2022, Nature, 612, 658, doi: 10.1038/s41586-022-05476-5
  * Yatsu et al. (2009) Yatsu, Y., Kawai, N., Shibata, S., & Brinkmann, W. 2009, PASJ, 61, 129, doi: 10.1093/pasj/61.1.129


# TTT: A Temporal Refinement Heuristic for Tenuously Tractable Discrete Time
Reachability Problems

Chelsea Sidrane  Jana Tumova C. Sidrane and J. Tumova are with the Division of
Robotics, Perception and Learning, Intelligent Systems Department, School of
Electrical Engineering & Computer Science, KTH Royal Institute of Technology,
Brinellvägen 8, 114 28 Stockholm, Sweden<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Reachable set computation is an important tool for analyzing control systems.
Simulating a control system can show that the system is generally functioning
as desired, but a formal tool like reachability analysis can provide a
guarantee of correctness. For linear systems, reachability analysis is
straightforward and fast, but as more complex components are added to the
control system such as nonlinear dynamics or a neural network controller,
reachability analysis may slow down or become overly conservative. To address
these challenges, much literature has focused on spatial refinement, e.g.,
tuning the discretization of the input sets and intermediate reachable sets.
However, this paper addresses a different dimension: temporal refinement. The
basic idea of temporal refinement is to automatically choose _when_ along the
horizon of the reachability problem to execute slow symbolic queries which
incur less approximation error versus fast concrete queries which incur more
approximation error. Temporal refinement can be combined with other refinement
approaches and offers an additional “tuning knob” with which to trade off
tractability and tightness in approximate reachable set computation. Here, we
introduce an automatic framework for performing temporal refinement and we
demonstrate the effectiveness of this technique on computing approximate
reachable sets for nonlinear systems with neural network control policies. We
demonstrate the calculation of reachable sets of varying approximation error
under varying computational budget and show that our algorithm is able to
generate approximate reachable sets with a similar amount of error to the
baseline approach in 20-70% less time.

3-5 keywords

##  1 Introduction

Every controller needs to be analyzed for performance and reliability.
Simulation is useful to get an idea of the general behavior of a control
policy, but only formal analysis can guarantee correctness with respect to a
model of the system and a specification. In this paper, we focus on reach-
avoid specifications which dictate a goal set the system must reach and/or an
avoid set it must not reach. For a given discrete time dynamical system with a
control policy and set of possible starting states the forward reachable set
at time $t$ is the set of states that the system could reach $t$ steps into
the future. Given a final set instead of an initial set, the $t$-step backward
reachable set is the set of states that will reach the final set after $t$
steps. For simple settings like linear dynamical systems and convex polytopic
initial sets reachability analysis is straightforward and fast [1]. However,
many modern control systems operate on nonlinear systems [2], and increasingly
contain neural network components. Reachability methods for such complex
systems may not exist, and if they do they they almost always require under or
over-approximation (also referred to as inner and outer approximation) of the
reachable set as the exact set is not computable. It is desirable to produce
tight reachable sets which have as little approximation error as possible, but
this often trades off with tractability – the computation time needed to
compute the reachable set.

Figure 1: An illustration of overapproximate forward reachability. A symbolic
query implicitly constrains the reachable set at intermediate timesteps
producing tighter final sets. This concept is central to the refinement
algorithm.

There have been many approaches in the reachability literature to address
tractability and tightness. One of the central concepts used is domain
refinement [3, 4]. Domain refinement refers to splitting the domain into
subsets, e.g., through gridding, and computing the reachable set for each
subset. Without domain refinement, the reachable set approximation would be
unusably loose, but if gridded too finely the problem can become intractable
as there is then exponential complexity in the number of state dimensions.

Another concept that has been used in the reachability literature is
decomposition [5, 6, 7]. Decomposition refers to re-writing the system
dynamics into a collection of uncoupled or loosely coupled subsystems. If
e.g., using gridding to produce tight set approximations, fine gridding can
then be used on each low dimensional subsystem without incurring the
exponential cost of gridding the original high dimensional system.

Lastly, many methods make use of pre-specified template sets in order to speed
up computation [8, 9, 10, 11]. In these methods, the tightest circumscribing
(or inscribing) set of a particular type such as a hyperrectangle, ellipsoid,
or zonotope is computed rather than computing the arbitrary non-convex true
reachable set.

In this work we make use of a concept known as symbolic reachability [8] to
define an algorithm for temporal refinement. Many verification methods use
concrete reachability which involves computing the approximate reachable set
at time $t+1$ using only the approximate reachable set at time $t$, or $t-1$
from $t$ in the backward case [12]. While concrete reachability is fast, it
accumulates large approximation error quickly, a phenomenon known as the
“wrapping effect” [13]. Symbolic reachability can be used to mitigate this
effect [8]. Using symbolic reachability, the approximate reachable set at
e.g., time $t+2$ is computed from the set at time $t$ while implicitly
constraining the system trajectories to pass through the intermediate set at
time $t+1$ (fig. 1). This results in tighter sets than if the set at $t+2$ was
computed using only the approximate set at time $t+1$. It is possible to
perform symbolic reachability for limited time horizons (temporal depths),
each iteration computing a set farther into the future (or further in the
past) and implicitly constraining system trajectories to pass through all sets
at timesteps in between, but the problem will quickly approach intractability.
To address this tradeoff between tightness and tractability, the authors of
[8] perform hybrid-symbolic reachability which uses mixtures of concrete and
symbolic computations in a hand-tuned, ad-hoc manner.

In this paper, we contribute a temporal refinement algorithm to automate
hybrid-symbolic reachability analysis for tenuously tractable discrete-time
dynamical systems. The temporal depth and end timestep of each query are
automatically selected to produce tight reachable sets given the remaining
computational budget and projected total computation length. Our algorithm is
an orthogonal contribution to existing strategies for balancing tightness and
tractability because it may be combined with concepts such as domain
refinement, system decomposition, and pre-specified template sets; among
others. 111Note that while temporal refinement may sound similar to timestep
adaptation for continuous systems [14], it is not closely related. In a
continuous time system, the length of time step $\Delta t$ is an algorithm
hyperparameter that maybe adapted according to error requirements but in a
discrete time system, there is either no specific step length, only time
indices, $t=0$, $t=1$ or there is a fixed timestep length $\Delta t$ which
cannot be adapted.

We demonstrate our algorithm for temporal refinement on overapproximate
forward reachability analysis, but it is similarly applicable to backward
reachability analysis and underapproximation. Specifically, we treat a
challenging class of problems: computing forward reachable sets for discrete
time nonlinear dynamical systems with neural network control policies; also
called Neural Feedback Loops (NFLs). For this class of problem, there has been
work on computing forward reachable sets [8], backward reachable sets [15],
performing domain refinement [16, 4, 17], and using exotic template polyhedron
[18] but there remains difficulty balancing tractability and tightness making
it a good class of problem on which to demonstrate our technique. We use
several examples from the literature and show that our approach can produce a
range of overapproximation error values depending on the allotted
computational budget, including computing reachable sets 20-70% faster than
previous hand-tuned approaches for the same amount of error in the reachable
set approximation.

##  2 Problem Definition

Consider a discrete time dynamical system with state $\vec{x}\in R^{d}$
governed by update equation $\vec{x}_{t+1}=f(\vec{x}_{t},\vec{u}_{t})$ with
$\vec{u}\in R^{m}$. In this paper, we focus on nonlinear update functions $f$
but linear systems may also be considered. The control input
$\vec{u}_{t}=c(\vec{x}_{t})$ comes from a feedback control policy.

We seek to compute the forward reachable sets $\mathcal{R}_{t}$ of the
controlled system at future timesteps $t\in 1,2,3,\ldots$. The forward
reachable set at time $t_{n}$, $n$ steps into the future from time $0$ is
defined as follows:

$\displaystyle\mathcal{R}_{n}\triangleq\\{\vec{x}_{n}\mid\vec{x}_{n}=f(\vec{x}_{n-1},c(\vec{x}_{n-1})),$
(1) $\displaystyle\vec{x}_{n-1}=f(\vec{x}_{n-2},c(\vec{x}_{n-2})),$
$\displaystyle\ldots,$
$\displaystyle\vec{x}_{1}=f(\vec{x}_{0},c(\vec{x}_{0})),$
$\displaystyle\vec{x}_{0}\in\mathcal{X}_{0}\\}$

where $\mathcal{X}_{0}$ is the initial set of states. As aforementioned,
methods exist to compute the exact reachable set $\mathcal{R}_{t}$ for classes
of linear or piecewise linear systems, but for arbitrary smooth nonlinear
dynamics, approximation (over or under) must be used. In this paper we denote
the overapproximate reachable set $\tilde{\mathcal{R}_{t}}$ where
$\mathcal{R}_{t}\subseteq\tilde{\mathcal{R}_{t}}$. The overapproximation error
$e$ at a given time $t$ is defined

$e=\frac{m(\tilde{\mathcal{R}_{t}})}{m(\mathcal{R}_{t})}$ (2)

where $m$ is an error metric and the desired value is $e=1$. The error over a
finite time horizon $t\in\\{1\ldots n\\}$ is defined

${e}=\frac{\sum_{t=1}^{n}m(\tilde{\mathcal{R}_{t}})}{\sum_{t=1}^{n}m({\mathcal{R}_{t}})}$
(3)

Thus the problem this paper addresses may be stated as follows. When
performing finite time hybrid-symbolic reachability analysis of $n$ steps for
a discrete time dynamical system, algorithmically select points
$t_{i}\in\\{1,\ldots,n\\}$ at which to perform symbolic reachability queries
and select temporal depths $b_{\text{steps}}\in\\{1\ldots t_{i}\\}$ for those
symbolic queries so as to minimize approximation error (eq. 3) while adhering
to a pre-specified computational budget $b$.

Figure 2: Symbolic reachability scheme for primary refinement. In this example
$b_{\text{steps}}=4$.

##  3 Refinement Algorithm

At a high level, the algorithm begins in the “search” phase and searches for
an estimate of the longest tractable temporal depth (symbolic horizon)
$b_{\text{steps}}$. The algorithm then moves to a “jump” phase and performs
symbolic queries of that temporal depth or shorter until the set at the
desired final time $n$ has been reached. Alg. 1 provides more detail. In alg.
1, $\mathcal{X}_{0}$ is the initial set for reachability analysis, $b$ is the
time budget in seconds, and $q$ is a query object consisting of a dynamics
function $f$, a neural network controller $c(x):=NN(x)$, a final time horizon
$n$, a variable to store current temporal depth $h$, and various parameters
(e.g., hyperparameters controlling tightness of the nonlinear
overapproximation). The method returns $\tilde{\mathbb{R}}$ a vector of
reachable sets $\tilde{\mathcal{R}}_{t},\forall t\in[1\ldots n]$.

The motivation for the structure of the algorithm comes partially from the
initial work on hybrid-symbolic reachability [8] which used long symbolic
jumps with one-step concrete queries to generate the sets in between. Long
symbolic jumps produce tight approximations of the reachable set and beginning
symbolic jumps from tight approximations limits the wrapping effect. However,
in this paper, a reasonable symbolic “jump” size is not known a priori so we
first search for a long jump size. The search procedure also has the benefit
of limiting reachable set growth. Small errors early in the time horizon
compound over many timesteps, and because the search phase produces tight
reachable sets for early steps, it limits error compounding.

Alg. 1 begins in the “search” phase with a temporal depth of
$b_{\text{steps}}=1$ (line 3) beginning from the set $\mathcal{X}_{0}$ at time
$t_{0}$ (lines 4-5). A symbolic reachability query is then performed (lines
7-8), and the new set(s) are then logged (lines 9-13). If there is enough time
budget to continue searching, $b_{\text{steps}}$ will increment by $1$ (line
16). Symbolic reachable set computation queries are performed from the initial
set $\mathcal{X}_{0}$ at time $t_{0}$ to increasingly large temporal depths /
numbers of steps $b_{\text{steps}}$ into the future until it is estimated that
trying a temporal depth one step longer would put the entire reachable set
procedure over budget. Once this occurs, the algorithm will enter the “jump”
phase and perform symbolic queries beginning from the set at
$t_{\text{start}}:=t_{\text{cur}}$ (line 18-19) of temporal depth
$b_{\text{steps}}$ or shorter until the set at the desired final time $n$ has
been reached (line 6). Line 21 keeps $b_{\text{steps}}$ within correct range;
and line 22 updates the early stopping timeout. The early stopping timeout
limits budget overruns and is hand-tuned. The algorithm returns the reachable
sets from time $t=1\ldots n$.

1 function  _RefinedReach(_$q$ , $\mathcal{X}_{0}$, $b$_)_:

2

3 set_initial_timeout!($b$, $q$)

4 phase $\leftarrow$ “search”; $b_{\text{steps}}\leftarrow 1$

5 $t_{\text{start}}\leftarrow$ 0; $t_{\text{cur}}\leftarrow$ 0;

6 $\mathcal{X}_{\text{start}}\leftarrow\mathcal{X}_{0}$;
$\tilde{\mathbb{R}}\leftarrow$ []

7 while  _$t_{\text{cur}} <n$_ do

8 q.h $\leftarrow$ $b_{\text{steps}}$

9 data, $\tilde{\mathbb{R}}_{\text{out}}$ $\leftarrow$
symbolic_reach($q,\mathcal{X}_{\text{start}}$)

10 if  _phase == “search”_ then

11 push!($\tilde{\mathbb{R}}$, $\tilde{\mathbb{R}}_{\text{out}}[end]$)

12

13 else

14 push!($\tilde{\mathbb{R}}$, $\tilde{\mathbb{R}}_{\text{out}}$)

15

16 end if

17 $t_{\text{cur}}\leftarrow t_{\text{start}}+b_{\text{steps}}$

18 $b\leftarrow b~{}-~{}$elapsed()

19 $b_{\text{steps}_{0}}$, phase $\leftarrow$
calc_steps($t_{\text{start}},b_{\text{steps}},b$, data, $t_{\text{est}}$, $n$,
phase, status(q))

20 if  _phase == “jump”_ then

21 $t_{\text{start}}\leftarrow t_{\text{start}}+b_{\text{steps}}$

22 $\mathcal{X}_{\text{start}}\leftarrow\tilde{\mathbb{R}}[t_{\text{start}}$]

23

24 end if

25 $b_{\text{steps}}\leftarrow\min(b_{\text{steps}_{0}},n-t_{\text{start}})$

26 update_timeout!($q$, $b$)

27

28 end while

29 return $\tilde{\mathbb{R}}$

30

Algorithm 1 Reachability for NFLs with automatic temporal refinement

One of the key pieces of Algorithm 1 lies within the subroutine calc_steps
(line 16) which determines whether to continue searching for a longer temporal
depth $b_{\text{steps}}$ or to move to the “jump” phase. The complete
algorithm for calc_steps is included in the Appendix (alg. 2) and described
here. calc_steps first estimates the amount of time needed compute one
symbolic step, $t_{\text{est}}$, assuming that the time to compute a symbolic
query of $m$ steps is roughly $t_{\text{est}}*m$. This is a simplification
that is not strictly accurate – the time a symbolic query takes may be in the
worst case exponential in temporal depth for some problems and timing depends
on the particular input set – but it leads to roughly accurate budgeting given
that $t_{\text{est}}$ is updated at each iteration. The algorithm then
estimates the time that future queries will take and uses this to determine if
there is enough time budget to continue searching. This time estimate includes
a symbolic query from $t_{0}$ to $b_{\text{steps}}+1$ plus the time to finish
finish the temporal horizon in the “jump” phase. If the remaining budget $b$
is large enough, the algorithm increments $b_{\text{steps}}\leftarrow
b_{\text{steps}}+1$ and if not, the algorithm enters the “jump” phase.
Finally, if the solver stops early, the phase is also switched to “jump”. The
early stopping is used to ensure the budget is approximately respected given
that the estimated time per symbolic step $t_{\text{est}}$ is not an upper
bound.

Note that long symbolic queries generally produce tight reachable sets but
shorter queries can produce tighter sets if the starting set of the shorter
query is tighter itself.

##  4 Overapproximate Forward Reachability

Our refinement algorithm is applicable to both forwards and backwards
reachability and both over (outer) and under (inner) approximation but in this
paper, we demonstrate our refinement algorithm on overapproximate forward
reachability of nonlinear neural feedback loops (NFLs). In particular, we use
the reachability approach presented in [8]. To briefly summarize, the approach
computes tight piecewise linear bounds for smooth nonlinear functions in the
dynamics and then encodes these bounds as well as the ReLU neural network
control policy into a mixed integer linear program. A hyperrectangular
template set is then optimized to obtain an overapproximation of the reachable
set. The mixture of symbolic and concrete queries is determined by the
refinement algorithm presented here, but the procedure to approximate
nonlinear functions and compute approximate reachable sets using mixed integer
programs comes from [8] where further details can be found.

##  5 Soundness, Completeness, Complexity

The algorithm is sound, i.e.
$\mathcal{R}_{t}\subseteq\tilde{\mathcal{R}}_{t}$, for all $t$, following
directly from the soundness of [8]. In terms of completeness, the algorithm
presented here does not guarantee that the overapproximate reachable sets
$\tilde{\mathcal{R}}_{t}$ converge to the true reachable sets
$\mathcal{R}_{t}$. The method presented here could in theory be combined with
methods that provably converge to the true reachable set for some classes of
nonlinear systems [19] to overcome this limitation, but completeness in this
sense is beyond the scope of this paper.

The algorithm is guaranteed to produce a sequence of reachable sets as even
with early stopping, the algorithm will extend the compute time until a finite
feasible solution with $<50\%$ relative duality gap has been found. The
algorithm is guaranteed to terminate under the mild assumption that a feasible
solution with $<50\%$ relative duality gap exists.

The time complexity of this algorithm is in the worst case NP-hard as it
involves solving mixed integer linear programs. However, the average case
complexity is in practice more reasonable as is demonstrated in the results
section.

##  6 Numerical Experiments

Numerical experiments were run on examples taken from [8, 20] to assess the
tightness and computational speed of the temporal refinement algorithm. Four
problems were used: the pendulum dynamical system with a controller of 2
layers, 25 neurons per layer (S1); the TORA dynamical system with a controller
of 3 layers, 25 neurons per layer (T1); the car dynamical system with a
controller of 1 layer with 100 neurons (C1) and a controller of 1 layer with
200 neurons (C2). The reader is referred to [8] for the details on the neural
network controlled dynamical systems.

###  6.1 Results

The temporal refinement algorithm was run multiple times to assess how
changing the time budget changes the approximate reachable set error. Figure 3
measures the approximation error of the sets computed by the optimizer using
the ratio of the sum of volumes:

$e_{\text{total}}=\frac{\sum_{t=1}^{n}\text{vol}(\tilde{\mathcal{R}}_{t})}{\sum_{t=1}^{n}\text{vol}(\text{\@text@baccent{$\mathcal{R}_t$}})}$
(4)

In fig. 3, one can observe that the approximate error generally decreases with
increasing budget. Note that the trend is not monotonic because the algorithm
is a heuristic.

We compare our algorithm to the hand tuned hybrid symbolic approach used in
[8] which is represented in fig. 3 as dotted lines. One might expect that the
hand tuned approach would produce approximate reachable sets more quickly for
a given amount of error. A significant amount of work was presumably done
offline to identify good symbolic intervals that must instead be identified at
runtime with our approach. However, it turns out that our heuristic can
actually produce approximate reachable sets with similar amounts of error as
the hand-tuned approach in 20-70% less time, as shown in fig. 3. Furthermore,
our approach is able to produce reachable sets for some problems such as the
pendulum problem S1 that have 40% less error than the hand-tuned approach when
allowed to compute for longer than the time taken by the hand-tuned approach.

Figure 3: Approximate reachable sets of various fidelity are produced
depending on allotted time and are competitive with hand-tuned baseline
indicated with dashed lines.

We also compare sets produced by our algorithm to concrete one step reachable
sets. Using one step sets is fast but leads to large approximation error.
Examples of the reachable sets computed by our algorithm and baselines can be
seen in fig. 4 and fig. 5. Set are compared both using eq. 4 as well as the
ratio of the sum of radii:

$e_{\text{total}}=\frac{\sum_{t=1}^{n}\big{(}\sum_{i=1}^{d}\text{radius}_{i}(\tilde{\mathcal{R}}_{t})\big{)}}{\sum_{t=1}^{n}\big{(}\sum_{i=1}^{d}\text{radius}_{i}(\text{\@text@baccent{$\mathcal{R}_t$}})\big{)}}$
(5)

One step concrete sets are labeled ’naive’ and one can observe that every set
produced by our approach is equal to or tighter than those produced with the
naive approach. To compare to the hand-tuned approach, we select a run from
fig. 3 that had roughly equivalent or less error in a smaller amount of time
and calculate the increase in speed by our approach. For the runs selected in
fig. 4 and fig. 5, not every set produced by our approach is tighter than sets
produced by the hand tuned approach, but when error is calculated over all
sets in the trajectory, the total error is roughly equal or smaller. For
pendulum problem S1 featured fig. 4, our algorithm is 20.9% faster than the
hand-tuned approach for less error by both radius and volume metrics, and for
the TORA problem T1 featured in fig. 5, our algorithm is 24% faster for an
$\approx$ equal amount of error (1.5% less error according to the ratio metric
and 1.2% more error according to the volume metric).

In fig. 6, we display numerical results for the remaining problems, car
problems C1 and C2. For C1 and C2, we select runs that have less error by both
radius and volume metrics and are respectively 70.8% faster and 65.2% faster,
which translates to 4.7 minutes faster and 4.1 minutes faster.

| Naive | Hand-Tuned | Ours
---|---|---|---
Error (volume) | 11.9 | 2.27 | 1.94
Error (radius) | 3.79 | 1.61 | 1.59
Time | 7.84s | 19.3s | 15.3s

Figure 4: Reachable sets for pendulum problem S1 with two layer, 25 neuron-
per-layer neural network controller show less error than the naive baseline
and greater speed than the hand-tuned baseline for less error.

| Naive | Hand-Tuned | Ours
---|---|---|---
Error (volume) | 31.7 | 7.67 | 7.78
Error (radius) | 1.81 | 1.41 | 1.38
Time | 6.62s | 63.0s | 47.8s

Figure 5: Dimensions 1 and 2 of the reachable sets for TORA problem T1 with three layer, 25 neuron-per-layer neural network controller which has less error than the naive baseline and greater speed than the hand-tuned baseline for similar error. Car (C1) | One-Step | Hand-Tuned | Ours
---|---|---|---
Error (volume) | 37.7 | 4.90 | 3.89
Error (radius) | 1.88 | 1.31 | 1.28
Time | 10.3s | 400s | 117s
Car (C2) |  |  |
Error (volume) | 37.9 | 5.29 | 4.43
Error (radius) | 1.85 | 1.30 | 1.29
Time | 11.1s | 378s | 131s
Figure 6: Reachable sets for the Car problems C1 and C2 are computed
significantly faster for comparable error using our algorithm as compared to
the hand-tuned baseline.

Lastly, our code will be released as an open source tool and be available
here.

###  6.2 Limitations

The approach described here has limitations related to budgeting. A poor
choice of initial solution from optimizer can lead to budget overruns as the
algorithm will continue to compute if a finite solution has not yet been
found. Further, budgeting with the linear time approximation as described is
not accurate for all reachability problems, leading to budget overruns and/or
conservative solutions for some problems; which may be exacerbated by long
horizon queries. As previously mentioned, our algorithm is a heuristic. There
does not exist a guarantee that given two budgets $b_{A}>b_{B}$ the sets
produced with budget $b_{A}$ will be tighter than those produced using
$b_{B}$.

##  7 Conclusion

In this paper we have introduced an algorithm for temporal refinement of
tenuously tractable discrete time reachability problems. The algorithm
computes reachable sets using constraints from several timesteps at once in
order to produce tighter sets. The temporal depth of timestep constraints
considered for a given reachable set is calculated based upon remaining
computational budget and estimated future computation. The algorithm
introduced here makes it possible to do hybrid-symbolic reachability without
hand tuning concrete and symbolic intervals. We demonstrated this algorithm on
a difficult class of reachability problems: nonlinear dynamical systems with
neural network control policies, and demonstrated that it is able to produce
reachable sets of varying error given varying computational budget.
Additionally, it is able to produce reachable sets 20-70% faster then hand
tuned approaches for the same amount of error. Ultimately, this algorithm and
associated open-source code represent tools to enable more efficient
reachability analysis for complex control systems.

## References

  * [1] G. Frehse, C. Le Guernic, A. Donzé, S. Cotton, R. Ray, O. Lebeltel, R. Ripado, A. Girard, T. Dang, and O. Maler, “Spaceex: Scalable verification of hybrid systems,” in _Computer Aided Verification: 23rd International Conference, CAV 2011, Snowbird, UT, USA, July 14-20, 2011. Proceedings 23_. Springer, 2011, pp. 379–395.
  * [2] J.-J. E. Slotine, “Applied nonlinear control,” _PRENTICE-HALL google schola_ , vol. 2, pp. 1123–1131, 1991.
  * [3] S. Bansal, M. Chen, S. Herbert, and C. J. Tomlin, “Hamilton-jacobi reachability: A brief overview and recent advances,” in _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_. IEEE, 2017, pp. 2242–2253.
  * [4] M. Everett, R. Bunel, and S. Omidshafiei, “Drip: domain refinement iteration with polytopes for backward reachability analysis of neural feedback loops,” _IEEE Control Systems Letters_ , 2023.
  * [5] M. Chen, S. L. Herbert, M. S. Vashishtha, S. Bansal, and C. J. Tomlin, “Decomposition of reachable sets and tubes for a class of nonlinear systems,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 11, pp. 3675–3688, 2018.
  * [6] I. M. Mitchell and C. J. Tomlin, “Overapproximating reachable sets by hamilton-jacobi projections,” _journal of Scientific Computing_ , vol. 19, pp. 323–346, 2003.
  * [7] X. Chen and S. Sankaranarayanan, “Decomposed reachability analysis for nonlinear systems,” in _2016 IEEE Real-Time Systems Symposium (RTSS)_. IEEE, 2016, pp. 13–24.
  * [8] C. Sidrane, A. Maleki, A. Irfan, and M. J. Kochenderfer, “OVERT: An algorithm for safety verification of neural network control policies for nonlinear systems,” _The Journal of Machine Learning Research_ , vol. 23, no. 1, pp. 5090–5134, 2022.
  * [9] A. B. Kurzhanski and P. Varaiya, “Ellipsoidal techniques for reachability analysis: internal approximation,” _Systems & control letters_, vol. 41, no. 3, pp. 201–211, 2000.
  * [10] W. Kühn, “Rigorously computed orbits of dynamical systems without the wrapping effect,” _Computing_ , vol. 61, pp. 47–67, 1998.
  * [11] M. Althoff and B. H. Krogh, “Zonotope bundles for the efficient computation of reachable sets,” in _2011 50th IEEE conference on decision and control and European control conference_. IEEE, 2011, pp. 6814–6821.
  * [12] J. A. Siefert, T. J. Bird, J. P. Koeln, N. Jain, and H. C. Pangborn, “Successor sets of discrete-time nonlinear systems using hybrid zonotopes,” in _2023 American Control Conference (ACC)_. IEEE, 2023, pp. 1383–1389.
  * [13] A. Neumaier, _The wrapping effect, ellipsoid arithmetic, stability and confidence regions_. Springer, 1993.
  * [14] M. Wetzlinger, A. Kulmburg, A. Le Penven, and M. Althoff, “Adaptive reachability algorithms for nonlinear systems using abstraction error analysis,” _Nonlinear Analysis: Hybrid Systems_ , vol. 46, p. 101252, 2022.
  * [15] N. Rober, S. M. Katz, C. Sidrane, E. Yel, M. Everett, M. J. Kochenderfer, and J. P. How, “Backward reachability analysis of neural feedback loops: Techniques for linear and nonlinear systems,” _IEEE Open Journal of Control Systems_ , 2023.
  * [16] N. Rober, M. Everett, S. Zhang, and J. P. How, “A hybrid partitioning strategy for backward reachability of neural feedback loops,” in _2023 American Control Conference (ACC)_. IEEE, 2023, pp. 3523–3528.
  * [17] W. Xiang, H.-D. Tran, X. Yang, and T. T. Johnson, “Reachable set estimation for neural network control systems: A simulation-guided approach,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 32, no. 5, pp. 1821–1830, 2020.
  * [18] D. Manzanas Lopez, T. Johnson, H.-D. Tran, S. Bak, X. Chen, and K. L. Hobbs, “Verification of neural network compression of acas xu lookup tables with star set reachability,” in _AIAA Scitech 2021 Forum_ , 2021, p. 0995.
  * [19] E. Asarin, T. Dang, and A. Girard, “Reachability analysis of nonlinear systems using conservative approximation,” in _International Workshop on Hybrid Systems: Computation and Control_. Springer, 2003, pp. 20–35.
  * [20] T. Johnson, D. Lopez, P. Musau, H. Tran, T. Botoeva, F. Leofante, A. Maleki, C. Sidrane, J. Fan, and C. Huang, “Arch-comp20 category report: Artificial intelligence and neural network control systems (AINNCS) for continuous and hybrid systems plants,” _EPiC Series in Computing_ , vol. 74, pp. 107–139, 2020.

###  .1 Calculating temporal depth for Algorithm LABEL:alg:initial_reach

Algorithm 2, calc_steps, calculates the next symbolic query for alg. 1. In
line 2 of alg. 2, calc_steps calculates the time index of the current step,
$t_{\text{cur}}$, and then in line 3, calc_steps makes a conservative estimate
of the amount of time needed compute one symbolic step, $t_{\text{est}}$. If
previously, phase = “search”, (line 4), and we have not stopped early (line
5), we will estimate if there is enough time budget to continue searching
(line 6).

To estimate if there is enough time budget to continue searching, we must
estimate the time that future symbolic queries will take. Assuming that the
time it takes to perform 1 symbolic step is roughly constant
($t_{\text{est}}$) and that we are in the search phase, the last symbolic
query was from step 0 of length $b_{\text{steps}}$, and the time cost of a
potential next query from time 0 to time $b_{\text{steps}}+1$ is $\approx
t_{\text{est}}*(b_{\text{steps}}+1)$. At any given time, the total remaining
computational cost of primary_refinement is the cost of the next search step
plus the cost of symbolic jumps to finish the time horizon in the “jump” phase
which can be estimated as taking time $(n-(t_{\text{cur}}+1))*t_{\text{est}}$.
As $t_{\text{cur}}=b_{\text{steps}}$ during the search phase, the total time
estimate for one further symbolic query from $0$ to $b_{\text{steps}}+1$ plus
jumps from $b_{\text{steps}}+1$ to $n$ is then
$t_{\text{est}}*(b_{\text{steps}}+1)+(n-(b_{\text{steps}}+1))*t_{\text{est}}=n*t_{\text{est}}$
(see line 6).

This estimate assumes that the time to compute $m$ symbolic steps is roughly
$t_{\text{est}}*m$ regardless of whether those steps are computed in 2
symbolic queries of size $\frac{t_{\text{est}}*m}{2}$ or 10 symbolic queries
of size $\frac{t_{\text{est}}*m}{10}$. This is a simplification we make that
is not strictly accurate – the time a symbolic query takes may be exponential
in temporal depth for some problems and timing depends on the particular input
set – but it leads to roughly accurate budgeting given that $t_{\text{est}}$
is updated at each iteration.

If the remaining budget $b$ does allow for continued search,
$b_{\text{steps}}$ is incremented (alg. 2 line 7). If not, we instead
calculate a reasonable temporal query depth for finishing the remaining steps
of the reachability problem in the jump phase (lines 8-13). One could use
jumps of size $b_{\text{steps}}$ but if $n$ is not divisible or close to
divisible by $b_{\text{steps}}$ this could lead to a short jump at the end of
the horizon and a loose final set. Instead, we calculate the num_jumps needed
to finish the time horizon in jumps of size $b_{\text{steps}}$ (line 9-10) and
instead use jumps of size
$\text{ceil}(\frac{s_{\text{left}}}{\text{num\\_jumps}})$ (line 11) which
makes all the jumps closer in size.

If the solver has stopped early (line 14), the phase is switched to “jump”
(line 18) and a similar calculation is performed to estimate the jump size
(lines 15-17) with the difference being num_jumps is calculated using
$b_{\text{steps}}-1$ instead of $b_{\text{steps}}$ (line 16). If we are
already in the “jump” phase (line 20) we do not change the jump size of
$b_{\text{steps}}$ unless we stop early (line 21) in which case we decrement
$b_{\text{steps}}$ by 1 (line 22). The early stopping is used to ensure the
budget is approximately respected given that the estimated time per symbolic
step $t_{\text{est}}$ is not an upper bound.

1 function  _calc_steps(_$t_{\text{start}},b_{\text{steps}},b$ , data, $n$,
phase, status_)_:

2 $t_{\text{cur}}\leftarrow t_{\text{start}}+b_{\text{steps}}$

3
$t_{\text{est}}\leftarrow\max(t_{\text{est}},\text{data}[t_{\text{cur}}].\text{time}/\text{data}[t_{\text{cur}}].\text{steps})$

4 if  _phase == “search”_ then

5 if  _status == “nominal”_ then

6 if  _$n*t_{\text{est}} <b$_ then

7 $b_{\text{steps}}\leftarrow b_{\text{steps}}+1$

8

9 else

10 $s_{\text{left}}\leftarrow n-b_{\text{steps}}$

11 num_jumps = $\text{ceil}(\frac{s_{\text{left}}}{b_{\text{steps}}})$

12
$b_{\text{steps}}\leftarrow\text{ceil}(\frac{s_{\text{left}}}{\text{num\\_jumps}})$

13 phase $\leftarrow$ “jump”

14 end if

15

16 else if  _status == “stopped early”_ then

17 $s_{\text{left}}\leftarrow n-b_{\text{steps}}$

18 num_jumps = $\text{ceil}(\frac{s_{\text{left}}}{b_{\text{steps}}-1})$

19
$b_{\text{steps}}\leftarrow\text{ceil}(\frac{s_{\text{left}}}{\text{num\\_jumps}})$

20 phase $\leftarrow$ “jump”

21 end if

22

23 else if  _phase == “jump”_ then

24 if  _status == “stopped early”_ then

25 $b_{\text{steps}}=\max(b_{\text{steps}}-1,1)$

26 end if

27

28 end if

29 return $b_{\text{steps}}$, phase

30

Algorithm 2 Procedure to calculate symbolic steps during primary refinement.