# Black hole triggered star formation in the dwarf galaxy Henize 2-10

Zachary Schutte
eXtreme Gravity Institute, Department of Physics
Montana State University
Bozeman, MT 59715, USA
<EMAIL_ADDRESS>
& Amy Reines
eXtreme Gravity Institute, Department of Physics
Montana State University
Bozeman, MT 59715, USA
<EMAIL_ADDRESS>

Black hole driven outflows have been observed in some dwarf galaxies with
active galactic nuclei manzano2019agn , and likely play a role in heating and
expelling gas (thereby suppressing star formation), as they do in larger
galaxies fabian2012observational . The extent to which black hole outflows can
trigger star formation in dwarf galaxies is unclear, because work in this area
has hitherto focused on massive galaxies and the observational evidence is
scarce gaibler2012jet ; maiolino2017star ; gallagher2019widespread . Henize
2-10 is a dwarf starburst galaxy previously reported to have a central massive
black holereines2011actively_Henize ; reines2012parsec ; reines2016deep ;
riffel2020evidence  though that interpretation has been disputed since some
aspects of the observational evidence are also consistent with a supernova
remnant hebbar2019x ; cresci2017muse . At a distance of $\sim$9 Mpc, it
presents an opportunity to resolve the central region and determine if there
is evidence for a black hole outflow impacting star formation. Here we report
optical observations of Henize 2-10 with a linear resolution of a few parsecs.
We find a $\sim$150 pc long ionized filament connecting the region of the
black hole with a site of recent star formation. Spectroscopy reveals a
sinusoid-like position-velocity structure that is well described by a simple
precessing bipolar outflow. We conclude that this black hole outflow triggered
the star formation.

Radio observations of Henize 2-10 using very long baseline interferometry
reveal a nuclear, compact, non-thermal source with a luminosity of
L${}_{R}\sim 4\times 10^{35}$ erg s-1 and a physical size $<$ 3 pc × 1 pc
reines2012parsec . High-resolution X-ray observations unveil a point source
with L${}_{X}\sim 10^{38}$ erg s-1 that is spatially coincident with the
compact nuclear radio source reines2016deep . There are two possible
explanations for these radio and X-ray observations alone - a highly sub-
Eddington massive black hole (i.e., a low-luminosity AGN) or a very young
supernova remnantreines2012parsec ; hebbar2019x . However, there are other
observational results to consider regarding the origin of the nuclear
radio/X-ray source in Henize 2-10. We summarize these results in Figure 4
(Extended Data Table 1) and demonstrate that a highly sub-Eddington massive
black hole is consistent with all the available observations including new
results presented here, while a supernova remnant is not.

We observed Henize 2-10 at optical wavelengths using the Space Telescope
Imaging Spectrograph (STIS) on the Hubble Space Telescope (HST). We obtained
observations of the central regions of Henize 2-10 with the 0.2"-slit in two
orientations. The first, referred to as the EW orientation, is centered on the
nuclear radio/X-ray source and aligned with the filamentary ionized structure
between the two bright, extended regions of ionized gas previously identified
in narrowband H$\alpha$ imaging with HST reines2011actively_Henize  (Figure
1). The second, referred to as the NS orientation, is centered on the nuclear
source and rotated 90∘ with respect to the EW observation. We obtained high-
dispersion observations (velocity resolution of $\sim 40$ km s-1) at both slit
positions using the G750M and G430M gratings, which cover the strong emission
lines of interest (e.g., H$\alpha$, [NII], [SII], [OI] and H$\beta$, [OIII],
respectively).

Figure 1: HST optical image of the dwarf starburst galaxy Henize 2-10 (Image
credit: NASA/STScI). The inset shows a narrowband H$\alpha+$continuum image of
the central 6” x 4” region. At the distance of Henize 2-10 ($\sim$9 Mpc), 1”
corresponds to $\sim$ 44 pc. White contours indicate compact radio emission
detected with very long baseline interferometry and mark the location of the
central massive black hole reines2012parsec , which is also detected in X-rays
reines2016deep . Our HST spectroscopy demonstrates that the black hole is
driving an outflow that is triggering the formation of young massive star
clusters. The main image is 25” ($\sim$1.1 kpc) across.

The kinematics of the ionized gas provide evidence for an outflow originating
from a nuclear massive black hole. First, we detect substantially broadened
emission lines at the location of the central source in both slit
orientations. In particular, the [OI]6300 emission line has a broad component
with a full width at half maximum (FWHM) of 497 km s-1 as measured in the EW
slit (and 445 km s-1 in the NS slit) using a 3-pixel extraction region in the
spatial dimension. While the [OI] line is too weak to be detected all along
the EW filament (or NS slit), it is detected at the location of the bright
star-forming region $\sim$70 pc to the east with a FWHM = 103 km s-1, much
less than that of the central source (e.g., see top left panel of Figure 2).
The [OIII]5007 line is much stronger than [OI]6300 and detectable all along
the EW filament. The FWHM of the broad component of [OIII] at the location of
the nuclear source is 271 km s-1, which is somewhat broader than the median
value per spatial pixel along the EW slit (175 km s-1 with a standard
deviation of 53 km s-1). We do not observe such broad emission beyond the
location of the central source in the NS slit (see top right panel of Figure
2). We emphasize that the line widths at the location of the nuclear source
are consistent with a low-velocity outflow from a massive black hole, but
would be anomalously low for a very young supernova remnant with typical FWHMs
of thousands of km s-1 mathewson1980new ; borkowski2017asymmetric .

In addition to broadened emission lines at the location of the central source,
we also find Doppler shifted velocities along the EW slit that exhibit a
coherent sinusoid-like pattern that is relatively smooth, especially in
comparison to the position-velocity diagram along the NS slit orientation that
shows no evidence of a coherent pattern in the Doppler velocities of strong
emission lines (see Figure 2 and Methods). Moreover, a simple model of a
precessing bipolar outflow broadly reproduces the observed sinusoid-like
velocity pattern along the EW slit that is aligned with the ionized filament
seen in narrowband H$\alpha$ imaging with HST (see Methods for a description
of the model). Precessing jets have been observed in many AGNs, although they
are typically found in more luminous quasars and radio galaxies
gower1982precessing ; dunn2006precession . Theories for the origin of
precessing jets/outflows include accretion disk warping, jet instabilities,
and the presence of massive black hole binaries pringle1996self ;
nixon2013jets . On the other hand, the coherent velocity pattern we observe
over $\sim$150 pc along the ionized filament centered on the central source is
incompatible with a supernova remnant origin since supernova remnants do not
drive quasi-linear outflows on such large scales.

Figure 2: Optical spectra and ionized gas kinematics for the central region of
Henize 2-10. Top panels show 2D continuum-subtracted HST STIS spectra of key
emission line regions along the EW and NS slit orientations. The middle panels
show position-velocity diagrams along each slit position where Doppler shifts
are relative to the systemic velocity of the galaxy (873 km s-1)
kobulnicky1995aperture . The velocity pattern in the EW slit is much more
coherent than that in the NS slit. The bottom panels show extracted 1D spectra
at the location of the nuclear massive black hole. The region surrounding the
nucleus (0 pc in both slit positions) shows strong broadened emission lines,
including [OI]6300 and [OIII]5007.

There is also evidence that the black hole outflow is triggering the formation
of star clusters in the central region of Henize 2-10. HST imaging shows that
the ionized filament extends eastward from the massive black hole to a bright
knot of ionized gas and site of recent star formation located 1.5" ($\sim$70
pc) away from the black hole (Figure 3). Given that our HST spectroscopy along
this filament exhibits a continuous velocity pattern, which can be tracked
from the black hole to the eastern star-forming region and is well described
by a precessing bipolar outflow model, this strongly suggests that the outflow
driven by the black hole is causally connected to the region of recent star
formation. There is also a secondary, blue-shifted peak (offset by 154 km s-1)
detected in the emission lines at the location of the bright star-forming
knot, suggesting the outflow is pushing the line-emitting gas clouds and
influencing their kinematics (Figure 3). The double peaked lines would
naturally arise as the outflow intercepted dense gas and primarily pushed it
in the lateral direction rather than ahead of the flow. The handful of young
star clusters associated with the star-forming knot have ages of $\sim$4 Myr
(see Methods) and are predominately aligned in the north-south direction,
consistent with a scenario in which they formed from gas moving in opposite
directions due to the impact of the black hole outflow kharb2017double . There
is also a local peak in the gas density at the location of the star-forming
knot with n${}_{e}\geq 10^{4}$ cm-3 (see Figure 7), consistent with the
outflow compressing gas clouds and enhancing star formation in this region. It
is notable that we also detect double peaked emission line profiles to the
west of the massive black hole at the boundary of a dark cloud of gas and dust
(region 5, see Figure 8), suggesting the other side of the bipolar outflow is
also intercepting dense clouds and pushing them perpendicular to the outflow
direction. Just to the west of this region (region 6, Figure 8), the
equivalent width of the H$\alpha$ emission line indicates stellar ages $\leq$
3 Myr, indicating a very early stage of star formation in which the infant
star clusters have not had enough time to destroy or disperse the dense
molecular gas from which they formed, leading to high levels of extinction
beck2018dense .

The ionization conditions of the gas provide additional support for a massive
black hole with an outflow. First, our HST spectroscopy at $\sim$0.1”
resolution clearly reveals non-stellar ionization at the location of the
central source based on the flux ratio of [OI]/H$\alpha$, which is consistent
with gas photoionized by an accreting massive black hole (Figure 9).
Additionally, supernova remnants are known to be strong [SII] emitters with
log([SII]/H$\alpha$) $>$ -0.521, which is not seen in the STIS spectrum of the
nuclear source. The ionization conditions probed by optical diagnostics using
[OIII]/H$\beta$ versus [NII]/H$\alpha$ and [SII]/H$\alpha$ do not clearly
indicate a (luminous) AGN, which at first glance my seemingly rule out the
presence of a massive black hole. However, there are a number of potential
reasons for this apparent discrepancy. First, [NII] and [SII] are not as
sensitive to the hardness of the ionizing radiation field as [OI] and are
therefore less reliable at identifying AGNs. Second, the massive black hole in
Henize 2-10 is accreting at a tiny fraction of its Eddington luminosity ($\sim
10^{-6}$ LEdd) reines2011actively_Henize , which can the impact emission lines
such that they do not look like those of luminous AGNs with higher Eddington
ratios trump2011accretion . The apparent contradiction of emission line
diagnostics can also arise when there is underlying star formation
contributing to the spectrum. Finally, we note that there are other examples
of radio-selected massive black holes in dwarf galaxies with enhanced [OI]
that do not look like optical AGNs in the other diagnostic diagrams
reines2020new ; molina2021outflows .

Figure 3: Visualization of the bipolar outflow model and star forming regions.
Top left: Position-velocity diagram of the [OIII]5007 emission line along the
EW slit where the Doppler shifts are relative to the systemic velocity of the
galaxy (873 km s-1) kobulnicky1995aperture . The data are well described by a
simple precessing outflow model with a precession frequency of f $\sim$ 5
revolutions Myr-1 and a precession angle of $\theta$ $\sim$ 4∘ (see Methods).
Bottom left: HST narrowband H$\alpha$+continuum image showing that the EW slit
position is aligned with the ionized filamentary structure connecting the
massive black hole and site of recent star formation 70 pc to the east. White
contours indicate compact radio emission from very long baseline
interferometry reines2012parsec , and show the location of the massive black
hole. Bottom right: HST 0.8 micron ($\sim$I-band) broadband image with the
same field of view. Young star clusters in the eastern star-forming region are
highlighted with an ellipse. Top right: Spectra of the eastern star-forming
region (-70 pc) show strong broad emission lines with a secondary blue shifted
peak most clearly visible in the [SII] emission.

We also find that the ionization conditions probed by strong emission line
ratios at the location of the central source and along the filament are well
described by theoretical models of shocks propagating in a high density medium
allen2008mappings  (see Methods). Specifically, emission line ratios in the
central $\sim$50 pc along the EW slit orientation are consistent with models
of a shock with a velocity of $\sim$200 km s-1 traveling through a gas with an
electron density of $\sim$1000 cm-3 and include a ‘precursor’ of ionizing
photons that travel upstream from the shock front pre-ionizing the gas. The
model shock velocities that match the strong emission line ratios agree with
those obtained from direct measurements of the line widths from our kinematic
study described above (e.g., [OIII]), and the model gas densities are
consistent with our measurements using the density sensitive line ratio
[SII]6716/6731 (see Methods). Moreover, these conditions are well explained by
an AGN-driven outflow mechanically exciting the interstellar medium in these
regions in addition to the precursor component contributing to the
photoionization of the interstellar medium. There is also evidence for shocked
emission at the location of the eastern star-forming region, particularly from
the secondary blue-shifted peak, in addition to photoionization from young hot
stars (see Methods). Combining these results with our kinematic study
indicates that the bipolar outflow generated by the central black hole is
shocking the interstellar medium in the central regions of Henize 2-10,
creating conditions that are favorable for positive AGN feedback
silk2009global ; silk2013unleashing .

AGN driven outflows have been discovered in a small sample of dwarf galaxies,
though the discovery of a black hole outflow in Henize 2-10 provides the first
example that is robustly spatially resolved. Moreover, Henize 2-10 differs
from these other systems in several other ways. The majority of black hole
driven outflows in dwarf galaxies have been found in galaxies with well-
defined nuclei and optically-selected AGNs with relatively high accretion
rates. These observations suggest the AGNs play a role in heating and
expelling gas in the galaxies and quenching star formation, a phenomenon known
as negative feedback penny2018sdss . This is in stark contrast with Henize
2-10, which has an irregular central morphology, is intensely forming stars,
and is experiencing positive feedback from a weakly accreting black hole that
is luminous at radio, rather than optical, wavelengths.

Indeed, “jet-mode” feedback is often associated with radio-loud AGNs accreting
at low Eddington ratios, and attributed to the unbound nature of radiatively
inefficient accretion flows trump2011spectropolarimetric ; santoro2020agn .
Warm ionized gas outflows are observed to accompany the radio jets/outflows in
some cases, particularly in young radio galaxies where nascent jets are
expanding through the interstellar medium in the central regions of their
hosts santoro2020agn . In these systems, the spatial extents of the warm
outflows (traced by emission lines) are similar to the radio morphologies. A
similar phenomenon is observed in Henize 2-10, which is most readily seen in a
comparison between the radio emission detected with the Very Large Array and
Pa$\alpha$ emission detected by HST in the central few hundred pc of the
galaxy reines2011actively_Henize . However, an important difference between
Henize 2-10 and powerful young radio galaxies is that the emission lines are
dominated by AGN photoionization in the more massive and luminous systems,
whereas the extended emission line regions in Henize 2-10 are dominated by
star formation (enhanced/triggered by the black hole outflow). Therefore,
Henize 2-10 may be a low-mass, low-power analog of young radio galaxies.

##  1 Methods

###  1.1 The Controversy – Massive Black Hole or Supernova Remnant:

In recent years, evidence has been mounting for a massive black hole powering
a low-luminosity active galactic nucleus (AGN) at the center of Henize 2-10
reines2011actively_Henize ; reines2012parsec ; reines2016deep ;
riffel2020evidence , although a supernova remnant has been proposed as an
alternative by some authors hebbar2019x ; cresci2017muse . As discussed in the
main text, the radio and X-ray point source luminosities are consistent with
both. A recent study hebbar2019x  argues for a supernova remnant based on
their findings that the X-ray spectrum is better fit by a hot plasma model
(typically used for supernova remnants) than a power-law model (typically used
for luminous AGNs). However, the soft X-ray spectrum of the nuclear source in
Henize 2-10 does resemble massive black holes accreting at very low Eddington
fractions constantin2009probing  including Sagittarius A∗ in the Milky Way
baganoff2003chandra . Another study using ground-based spectroscopy favored a
supernova remnant origin for the central source based on a lack of any AGN
ionization signatures cresci2017muse . However, the ground-based observations
used in that work had an angular resolution of $\sim$0.7", which is not
sufficient to cleanly isolate the weakly accreting black hole from nearby
young ($<$5 Myr) massive (M${}_{*}\sim 10^{5}M_{\odot}$) star clusters that
dominate the line ratios at this relatively course angular resolution. There
are other observational results to consider regarding the origin of the
nuclear radio/X-ray source in Henize 2-10. For example, a recent study using
adaptive optics integral field spectroscopy provides evidence for gas excited
by an AGN and an enhanced stellar velocity dispersion at the location of the
nuclear source consistent with a $\sim 10^{6}M_{\odot}$ black hole, favoring
the low-luminosity AGN interpretation riffel2020evidence . There is also
evidence for moderately significant variability on hour-long timescales in the
X-ray light curve, which is incompatible with a supernova remnant
reines2016deep ; hebbar2019x . Moreover, it is reasonable to expect that
Henize 2-10 hosts a massive black hole since its overall structure resembles
an early-type galaxy (albeit with a central starburst) and its stellar mass
may be as high as M${}_{*}\sim 10^{10}M_{\odot}$ nguyen2014extended , a regime
where the black hole occupation fraction is near unity greene2020intermediate
. The central starburst complicates the identification of the weakly accreting
black hole, yet a variety of multiwavelength observational results taken
collectively strongly support its presence. These results are summarized in
Figure 4 (Extended Data Table 1). A highly sub-Eddington massive black hole is
consistent with all of the observations, including the new work presented
here, while a supernova remnant is not. The present study not only adds to the
evidence for a massive black hole in Henize 2-10, it also demonstrates that a
bipolar outflow from the black hole is enhancing/triggering star formation in
its vicinity.

Figure 4: (Extended Data Table 1) Summary of observational results regarding
the nature of the nucleus in Henize 2-10. A highly sub-Eddington massive black
hole is consistent with all the available observations including new results
presented here, while a supernova remnant is not.

###  1.2 STIS Observations and Data Reduction:

Spatially resolved spectroscopic observations of the nuclear regions of the
dwarf starburst galaxy Henize 2-10 were obtained using the Space Telescope
Imaging Spectrograph (STIS) instrument on the Hubble Space Telescope (HST). We
obtained observations with two slit orientations. The first is aligned with
the quasi-linear ionized gas structure identified by Reines et
al.reines2011actively_Henize  and covers the central radio/X-ray source and
the bright knot of ionized gas to the east. We refer to this as the East-West
(EW) orientation. The second slit orientation was placed perpendicular to the
EW observation at the location of the central radio/X-ray source. We refer to
this as the North-South (NS) orientation. The candidate AGN itself was too
faint to acquire directly, therefore we used a target acquisition with an
offset from a bright point source 7.9" to the southeast. Spectra were taken
with the G750M and G430M gratings providing medium spectral resolution (R
$\sim$ 5000-6000) coverage of key optical emission lines. The central
wavelengths were set at 6581 Å and 4961 Å for the G750M and G430M gratings,
respectively. At each slit orientation, two orbits were spent in G430M and one
orbit in G750M. The observations were taken with a two-point dither pattern
with CR-SPLITS (multiple exposures taken to aid in cosmic ray rejection) at
each position to help eliminate cosmic rays. The calibrated dithered images
were combined and have a spatial resolution of $\sim$0.1", which corresponds
to a physical scale of $\sim$4 pc at the distance of Henize 2-10 ($\sim$9
Mpc).

Figure 5: (Extended Data Figure 1) Raw 2D spectra showing the [OI]6300
emission line at the location of the nucleus in the EW slit orientation. The
location of the nucleus is indicated by white circles and the two images
correspond to the two dithered sub-exposures.

###  1.3 Emission Line Fitting:

Before we fit emission lines in the spectra, we first modeled and removed the
continuum. The reduced spectra were continuum-subtracted by masking emission
line regions and fitting a low order polynomial to the continuum in each row
in the spatial dimension. A low order polynomial was used given the lack of
absorption features in the spectra. We did, however, consider the potential
impact of stellar absorption lines on our measurements and found that the
absorption line strengths are negligible compared to the emission line
strengths. Scaling a Starburst99 starburst99  model for a 4 Myr stellar
population (see Stellar Ages section below) to our observed spectra, we find
that the flux of H$\alpha$ absorption is smaller by a factor of 61 than the
H$\alpha$ emission and the flux of the H$\beta$ absorption is smaller by a
factor of 7 than the H$\beta$ emission at the location of the central source.
Accounting for this absorption has a negligible effect on the line ratios of
the nuclear source and does not impact the classifications based on the
diagnostic diagrams. Once the spectra were continuum-subtracted, we fit each
emission line of interest with a linear combination of Gaussian profiles to
characterize the flux and estimate kinematic properties along the spatial
dimension of each slit. The fitting was done using lmfit newville2016lmfit , a
non-linear least squares curve fitting package in Python. We fit each emission
line with up to two Gaussian components when needed. To determine if a second
Gaussian component is warranted, we require that the flux of both components
be greater than the 3$\sigma$ error of the flux. This process is performed row
by row in the spatial direction along each slit for emission lines of
interest. During this process we fit the H$\beta$, [OIII]5007, H$\alpha$,
[NII]6548/6583 and [SII]6716/6731 emission lines. We fit the H$\alpha$ and
[NII] lines simultaneously, fixing the spacing between [NII] Gaussian
components to their corresponding component in the H$\alpha$ emission line to
laboratory values. Additionally, we tie the widths of [NII] components and fix
the flux ratio of the [NII] lines to the laboratory value of 1:2.96.
Similarly, the two [SII] emission lines are fit simultaneously with the
spacing between Gaussian components of the two lines held fixed and the widths
of the Gaussian components are tied together. We also fit [OI]6300 in the
spectra of the nuclear source and the eastern star-forming region, but the
line is too weak to be detected all along the slits. Since the [OI]6300 line
has a complex profile at the location of the black hole, with possible double
peaked narrow lines and a much broader component than the other emission
lines, we confirmed that this was not due to an artifact in the data. In
Figure 5, we show close-up views of the raw 2D spectra along the EW slit
position. The two images correspond to the two dithered sub-exposures offset
by 7.5 pixels. The broad [OI] line at the location of the central source is
seen in both images at different positions on the detector (indicated by white
circles). Note that the locations of hot pixels do not change between the two
images. In Figure 6 we also show the final reduced 2D image with the sub-
exposures combined. The broad, double peaked nature of the [OI] line is
clearly visible. We also note that [OI] is similarly broadened in the nuclear
spectrum extracted from the NS slit position, although there is not an obvious
double-peaked narrow line component (see Figure 2 in the main paper). In any
case, broadened [OI] is clearly detected at the location of the nuclear source
in both slit positions indicating an outflow on the order of 500 km s-1.

Figure 6: (Extended Data Figure 2) Combined 2D spectra showing the [OI]6300
emission line at the location of the nucleus in the EW slit orientation. Same
as Figure 5 but showing the reduced 2D image with the dithered sub-exposures
combined.

###  1.4 Gas Density:

We estimate the electron density, ne, using the density sensitive line ratio
[SII]6716/6731 osterbrock2006astrophysics . This ratio is sensitive to
electron densities in the range of $\sim$102-104 cm-3. Along the EW slit
orientation, we find a range of n${}_{e}\sim 10^{2.5}-10^{4}cm^{-3}$,
indicating a relatively high-density gas (see Figure 7). These density
estimates are in general agreement with the gas densities predicted by the
Allen et al. allen2008mappings  shock/shock+precursor models in the central
regions of Henize 2-10 as described in the next section.

Figure 7: (Extended Data Figure 3) The electron density, ne, along the EW slit
orientation. We measure the electron density along the EW slit from the ratio
of [SII]6716/[SII]6731 and find the electron density ranges from $\sim
10^{2.5}-10^{4}cm^{-3}$ , which is within the range the [SII] ratio is
sensitive to density. The high densities are consistent with those predicted
by optical emission line diagnostics derived from the Allen et al.
allen2008mappings  shock models.

###  1.5 Emissin Line Diagnostics - Photoionization and Shock Models:

To understand the ionization mechanisms in the central regions of Henize 2-10
we compare our emission line measurements in various regions to
photoionization and shock models. In addition to the central radio/X-ray
source, we identified 7 regions of interest that are shown in Figure 8 and
serve to provide a spatially resolved picture of the kinematics and ionization
conditions in the central regions of Henize 2-10. The extraction regions taken
along the EW slit orientation were chosen to correspond with emission features
seen in the H$\alpha$ and I-band imaging from HST (young star clusters, knots
of ionized gas) as well as features seen in the STIS spectroscopy (broad
emission, double peaks).

Figure 8: (Extended Data Figure 4) The spatial extraction regions taken along
the EW slit orientation. We place these regions on optical emission line
diagnostic diagrams (Figures 9-11). Top panel: the extraction regions are
shown on the narrow band H$\alpha$ \+ continuum image from HST to highlight
the ionized gas features that several of the spatial extractions probe. Bottom
panel: the extraction regions are shown on the archival 0.8 micron HST image,
showing young star clusters that the EW slit orientation passes through.

We first utilize the standard emission line diagnostic diagrams described by
Baldwin et al. baldwin1981classification  and Veilleux et al.
veilleux1987spectral  that have been expanded upon and summarized in Kewley et
al. kewley2006host . An accreting BH will produce a much harder continuum than
is emitted by hot stars, and these diagrams take advantage of this fact by
comparing strong emission line ratios that are close together in frequency to
mitigate reddening effects. In this study we employ widely used emission line
diagnostic diagrams that take [OIII]/H$\beta$ versus [NII]/H$\alpha$,
[SII]/H$\alpha$, and [OI]/H$\alpha$ (see Figure 9). In the [NII]/H$\alpha$
diagram, line-emitting galaxies separate into a V-shape kewley2006host  with
star forming galaxies occupying the left most plume while AGNs occupy the
right branch of galaxies. These regions are quantified by an empirical
division between HII regions and emission from AGNs developed by Kauffmann et
al. kauffmann2003host . The “composite” region between this empirical division
and the theoretical maximum starburst line from Kewley et al.
kewley2001theoretical  indicates there is likely significant emission form
both HII regions and AGNs. Like the [NII]/H$\alpha$ diagram, the
[SII]/H$\alpha$ and [OI]/H$\alpha$ diagnostics provide diagnostics for
differentiating between emission from HII regions and AGNs. These two diagrams
add a dividing line to distinguish between emission from Seyferts and Low
Ionization Nuclear Emission Regions (LINERs). LINER emission can be generated
both by shocks and very hard AGN spectra and determining the primary
ionization mechanism can be complicated baldwin1981classification . It should
be noted that while these diagnostic diagrams are useful for identifying
regions dominated by luminous AGNs, they have limitations and can yield
ambiguous results for (or completely miss) massive black holes accreting at
very low Eddington ratios such as the one in Henize 2-10. Indeed, non-stellar
ionization is clearly indicated in the [OI]/H$\alpha$ diagram at the location
of the nuclear source, but not so for the other diagnostic diagrams (see
discussion in the main paper). Figure 9 shows the [OIII]/H$\beta$ versus
[NII]/H$\alpha$, [SII]/H$\alpha$, and [OI]/H$\alpha$ diagnostic diagrams for
the various extraction regions along the EW slit.

Figure 9: (Extended Data Figure 5) Narrow emission line diagnostic diagrams
showing various extraction regions along the EW slit orientation (see Figure
8). The nucleus (yellow point) falls in the Seyfert region of the
[OI]/H$\alpha$ diagram. The young star-forming region $\sim$70 pc to the east
of the low-luminosity AGN is depicted with a blue triangle and star for the
primary emission line component and the blue-shifted secondary component,
respectively. [OI] is not detected in all of the regions.

In addition to the diagnostic diagrams discussed above, we also investigate
whether the ionization conditions seen in the central regions of Henize 2-10
can be explained by mechanical excitation from shocks. To investigate this, we
employ ionization models of shock and shock+precursor emission developed by
Allen et al. allen2008mappings . These models provide emission line fluxes for
ionization from a pure shock (possibly driven by an outflow from an AGN or by
regions of intense star formation), where the gas is collisionally ionized, or
a shock+precursor where ionizing photons produced in the shock-heated gas
travel upstream and ionize the gas before the shock reaches it. We explore
models with a variety of electron densities (0.01-1000 cm-3), shock velocities
(100-600 km s-1) and transverse magnetic field strengths (0.01-32 $\mu$G).
These are shown in Figures 10 and 11.

Figure 10: (Extended Data Figure 6) Optical emission line diagnostics from the
shock and shock+precursor models with varying gas density. We place the
spatial extractions from the EW slit orientation shown in Figure 8 on a grid
of shock excitation models (presented in Allen et al. allen2008mappings ) with
varying gas density (n = 0.01-1000 cm-3) and shock velocity (v = 100-600
km/s). We fix the transverse magnetic field to be b = 1$\mu$G and the assume
solar metallicity.

The emission line ratios from the nuclear source are best described by the
shock+precursor models with a low shock velocity (100-250 km s-1), a high-
density gas (n = 1000 cm-3), and a low transverse magnetic field parameter (b
= 0.01 – 1 $\mu$G) (Figures 10 and 11). Shock+precursor models are thought to
be a good description for AGN+outflow emission. Along the filament (extraction
regions 2-5), the line ratios are explained well by a low velocity shock or
shock+precursor model ($\sim$200 km s-1) in a high density (n${}_{e}\sim 1000$
cm-3) gas with a transverse magnetic field parameter in the range of 1-10
$\mu$G. This is consistent with a scenario where the central black hole is
driving a bipolar outflow that shocks the gas and dominates the ionization
conditions along the filament.

Figure 11: (Extended Data Figure 7) Optical emission line diagnostics from the
shock and shock+precursor models with varying magnetic field. The models
(presented in Allen et al. allen2008mappings ) are shown as a grid with dashed
blue lines indicating constant shock velocity and dashed black lines
indicating constant transverse magnetic field. For these models, the density
is fixed to n = 1000 cm-3 and the transverse magnetic field parameter is
allowed to vary from b = 0.01-32 $\mu$G.

At the region of intense star formation located $\sim$70 pc to the east of the
black hole, we observe strong emission lines including a secondary, blue-
shifted kinematic component. To properly fit the emission lines in this
region, an additional Gaussian component is added (see Figure 3 in the main
paper). The separation of this component is determined using the [SII]
emission line region where the secondary peak is most clearly resolved. We
find that the secondary peak is offset from the primary peak by 154 km s-1,
and we use this value when fitting the other emission line regions where the
secondary peak is not as well resolved. When comparing to shock(+ precursor)
models, we find that the [NII]/H$\alpha$ diagram indicates a low shock
velocity (100-250 km s-1) in a high-density gas (n = 1000 cm-3) with a low
transverse magnetic field parameter (b = 0.01 – 1 $\mu$G) (see Figures 10 and
11). The [SII]/H$\alpha$ diagram shows the primary emission peak is
inconsistent with emission from shocks or shocks+precursor models. This
indicates that the primary emission peak at this location is primarily due to
star formation. The secondary emission peak is consistent with shock+precursor
models for the low velocity, high density conditions, indicating that this
kinematically distinct emission component is dominated by a shock+precursor
from the AGN-driven outflow.

Emission from extraction regions 6 and 7 (i.e., the western star-forming
region) is not consistent with any shock or shock+precursor models, which is
in agreement with their location in the HII region of the
Baldwin–Phillips–Terlevich (BPT) diagram. The line ratios are dominated by
star formation in this region.

###  1.6 Star Cluster Ages:

We estimate the ages of the young stellar clusters that fall within the EW
slit from their H$\alpha$ and H$\beta$ equivalent widths. To ascertain ages
from these equivalent widths we employ simple stellar population (SSP) models
from Starburst99 starburst99 . We use models from Version 7.1 with solar
metallicity (appropriate for the central regions of Henize 2-10 martin2006high
), an instantaneous burst of 104 M⊙ with a Kroupa IMF (0.1 – 100 M⊙), the
Geneva evolutionary tracks with high mass loss and the Pauldrach-Hillier
atmospheres.

At the location of the young stellar clusters in the eastern star-forming
region (region 1 in Figure 8) we find an equivalent width of 478 Å and 70 Å
for H$\alpha$ and H$\beta$ respectively. These both give stellar population
age estimates of $\sim$4.3 Myr, which is in good agreement with previous
estimates of the ages of other young star clusters in the region
chandar2003stellar . The ages of these clusters are larger than the crossing
time for the AGN-driven outflow ($\sim$0.3 Myr), based on the minimum outflow
velocity measured from emission line spectra ($\sim$200 km s-1) and the
distance between the AGN and the eastern star-forming knot ($\sim$70 pc).
Therefore, the timescales allow for the AGN-driven outflow to have
triggered/enhanced the formation of star clusters in the Eastern star-forming
knot.

The EW slit orientation also passes through a young star cluster in region 3.
At this location we find equivalent widths of 212 Å and 41 Å for H$\alpha$ and
H$\beta$ respectively, both indicating an age of $\sim$5.2 Myr for the stellar
cluster. Finally, region 6 in the western star-forming region has H$\alpha$
and H$\beta$ equivalent widths of 1092 Å and 196 Å, respectively. These
equivalent widths indicate the stellar clusters have ages $\leq$3 Myr.

###  1.7 Bipolar Outflow Model:

Here we provide a derivation of the model used to describe a precessing
bipolar outflow emanating from the central radio/X-ray source, which can
explain the coherent velocity structure seen in the central $\sim$120 pc of
the EW orientation observations. In this model we align the EW slit
orientation with the $z$-axis and assume the outflow precesses about this axis
with a small angle $\theta$ and an angular precession frequency $\omega$. If
the gas being ejected by the outflow has velocity $v_{0}$, and we orient the
$x$-axis to be in the direction of the observer, then the radial (Doppler
shifted) velocity seen at the location of the AGN as a function of time will
be given by

$v_{r}(t)=v_{x}(t)=v_{0}\sin{\left(\theta\right)}\sin{\left(\omega
t+\gamma\right)},$ (1)

where $\gamma$ represents a phase shift that accounts for small asymmetries in
the outflow profile (see Figure 12 for an illustration of our model).

In order to find the radial velocity as a function of distance ($z$) along the
slit axis, we must consider what angle the outflow made with the (line-of-
sight) $x$-axis when the gas at distance $z$ was emitted. Since this angle is
time dependent as the outflow precesses, the line-of-sight velocity of the gas
will depend on the orientation of the outflow at some time $t_{0}$ in the
past. The time that has passed since the gas at distance $z$ was ejected by
the outflow is determined by the $z$ velocity of the gas. Due to the symmetry
of the model about the $z$-axis, the $z$ component of the gas velocity will be
time independent and only depend on the angle of the outflow with the
$z$-axis,

$v_{z}=v_{0}\cos{(\theta)}.$ (2)

The time, $t_{0}$, for gas to reach a distance $z$ along the slit is then
given by

$t_{0}=\frac{z}{v_{0}\cos(\theta)}.$ (3)

We are then able to find the an expression for $v_{r}(z)$ by evaluating the
expression for $v_{r}(t)$ at the time $-t_{0}$:

$v_{r}(z)=v_{0}\sin{\left(\theta\right)}\sin{\left(\gamma-\frac{\omega}{v_{0}\cos{(\theta)}}z\right)}$
(4)

To fit this model to the data we require a rough estimate of the bulk outflow
velocity, $v_{0}$. We estimate this parameter using $W80$, the velocity
interval containing 80% of the line flux, of the broad emission seen at the
location of the candidate AGN. We find $W80\approx 200-500$ km s-1 based on
measurements of the [OIII]5007 and [OI]6300 lines at the location of the
candidate AGN. This allows us to determine the best-fit angle of precession,
$\theta$, and the frequency of precession, $f=\omega/2\pi$, to be

$\displaystyle\theta$ $\displaystyle=2.4^{\circ}-6.1^{\circ}$ (5)
$\displaystyle f$ $\displaystyle=3.0-7.5\ {\rm revolutions\ Myr}^{-1}$ (6)

where the larger angle of precession and smaller frequency of precession
corresponds to lower outflow velocities. We find consistent results when using
the Doppler shift profile of H$\alpha$, H$\beta$ and [OIII] emission lines to
fit the model derived above (the results using the [OIII] emission line is
shown in Figure 3). The Doppler shift profile can be coherently traced out to
50-60 pc on either side of the candidate AGN, most definitively out to the
bright eastern star forming region after which the Doppler shifts of the
emission lines are influenced by the bright young stars and then shown no
coherent pattern further along the slit. The coherent velocity structure seen
on the scale of 100 pc is not consistent with young supernova remnant as the
compact radio/x-ray source in the central regions of Henize 2-10, which
provides further motivation that a low luminosity AGN is driving the outflow.

These results are roughly consistent with jet parameters derived in other
studies where precessing or reorienting jet models have been applied. The long
precession period we observe ($\sim$200,000 years) is shorter by a factor of a
few than those seen predicted by Dunn et al. dunn2006precession , Nawaz et al.
nawaz2016jet  and Cielo et al. cielo2018feedback  but longer by a factor of a
few than those predicted by Gower et al. gower1982precessing  when jet
precession is invoked to explain the complex bending and knotting seen in
large radio jets.

Figure 12: (Extended Data Figure 8) A diagram of the toy model of the bipolar
outflow generated by the low-luminosity AGN in Henize 2-10. Our simple model
depends on the outflow velocity of the ionized gas (voutflow), the angle the
outflow makes with its precession axis ($\theta$) and the angular frequency
with which the outflow precesses ($\omega$). Similar models have been used to
describe the bending seen in large radio jets gower1982precessing ;
dunn2006precession .

##  2 Acknowledgments

We are grateful to Mallory Molina for useful discussions regarding shocks. We
also thank Mark Whittle and Kelsey Johnson for their assistance with the
HST/STIS proposal while AER was a graduate student at the University of
Virginia, as well as subsequent discussions. Support for Program number HST-
GO-12584.006-A was provided by NASA through a grant from the Space Telescope
Science Institute, which is operated by the Association of Universities for
Research in Astronomy, Incorporated, under NASA contract NAS5-26555. AER also
acknowledges support for this work provided by NASA through EPSCoR grant
number 80NSSC20M0231. ZS acknowledges support for this project from the
Montana Space Grant Consortium.

##  3 Author contributions statement

ZS reduced and analyzed the STIS data and compared the results to models. AER
led the HST/STIS proposal and helped with the data reduction. Both authors
worked on the interpretation of the results and writing of the paper.

##  4 Data Availability Statement

The spectroscopic data analyzed in this study are available from the Mikulski
Archive for Space Telescopes (MAST), https://archive.stsci.edu/

##  5 Competing Interest Statement

The authors declare no competing interests.

## References

  * [1] Mark G Allen, Brent A Groves, Michael A Dopita, Ralph S Sutherland, and Lisa J Kewley.  The mappings iii library of fast radiative shock models.  The Astrophysical Journal Supplement Series, 178(1):20, 2008.
  * [2] Frederick K Baganoff, Y Maeda, M Morris, MW Bautz, WN Brandt, W Cui, JP Doty, ED Feigelson, GP Garmire, SH Pravdo, et al.  Chandra x-ray spectroscopic imaging of sagittarius a* and the central parsec of the galaxy.  The Astrophysical Journal, 591(2):891, 2003.
  * [3] Jack A Baldwin, Mark M Phillips, and Roberto Terlevich.  Classification parameters for the emission-line spectra of extragalactic objects.  Publications of the Astronomical Society of the Pacific, 93(551):5, 1981.
  * [4] Sara C Beck, Jean L Turner, and S Michelle Consiglio.  Dense molecular filaments feeding a starburst: Alma maps of co (3–2) in henize 2–10.  The Astrophysical Journal, 867(2):165, 2018.
  * [5] Kazimierz J Borkowski, Peter Gwynne, Stephen P Reynolds, David A Green, Una Hwang, Robert Petre, and Rebecca Willett.  Asymmetric expansion of the youngest galactic supernova remnant g1. 9+ 0.3.  The Astrophysical Journal Letters, 837(1):L7, 2017.
  * [6] Rupali Chandar, Claus Leitherer, Christy Tremonti, and Daniela Calzetti.  The stellar content of henize 2-10 from space telescope imaging spectrograph ultraviolet spectroscopy.  The Astrophysical Journal, 586(2):939, 2003.
  * [7] S Cielo, A Babul, V Antonuccio-Delogu, J Silk, and M Volonteri.  Feedback from reorienting agn jets-i. jet–icm coupling, cavity properties and global energetics.  Astronomy & Astrophysics, 617:A58, 2018.
  * [8] Anca Constantin, Paul Green, Tom Aldcroft, Dong-Woo Kim, Daryl Haggard, Wayne Barkhouse, and Scott F Anderson.  Probing the balance of agn and star-forming activity in the local universe with champ.  The Astrophysical Journal, 705(2):1336, 2009.
  * [9] Giovanni Cresci, Leonardo Vanzi, Eduardo Telles, Giorgio Lanzuisi, Marcella Brusa, Matilde Mingozzi, Marc Sauvage, and Kelsey Johnson.  The muse view of he 2-10: No agn ionization but a sparkling starburst.  Astronomy & Astrophysics, 604:A101, 2017.
  * [10] RJH Dunn, AC Fabian, and JS Sanders.  Precession of the super-massive black hole in ngc 1275 (3c 84)?  Monthly Notices of the Royal Astronomical Society, 366(3):758–766, 2006.
  * [11] Andrew C Fabian.  Observational evidence of active galactic nuclei feedback.  Annual Review of Astronomy and Astrophysics, 50:455–489, 2012.
  * [12] Volker Gaibler, Sadegh Khochfar, Martin Krause, and Joseph Silk.  Jet-induced star formation in gas-rich galaxies.  Monthly Notices of the Royal Astronomical Society, 425(1):438–449, 2012.
  * [13] Robert Gallagher, Roberto Maiolino, Francesco Belfiore, Niv Drory, Rogério Riffel, and R A Riffel.  Widespread star formation inside galactic outflows.  Monthly Notices of the Royal Astronomical Society, 485(3):3409–3429, 2019.
  * [14] AC Gower and JB Hutchings.  A precessing relativistic jet model for 3c 449.  The Astrophysical Journal, 258:L63–L66, 1982.
  * [15] Jenny E Greene, Jay Strader, and Luis C Ho.  Intermediate-mass black holes.  Annual Review of Astronomy and Astrophysics, 58:257–312, 2020.
  * [16] Pavan R Hebbar, Craig O Heinke, Gregory R Sivakoff, and Aarran W Shaw.  X-ray spectroscopy of the candidate agns in henize 2–10 and ngc 4178: likely supernova remnants.  Monthly Notices of the Royal Astronomical Society, 485(4):5604–5615, 2019.
  * [17] Guinevere Kauffmann, Timothy M Heckman, Christy Tremonti, Jarle Brinchmann, Stéphane Charlot, Simon DM White, Susan E Ridgway, Jon Brinkmann, Masataka Fukugita, Patrick B Hall, et al.  The host galaxies of active galactic nuclei.  Monthly Notices of the Royal Astronomical Society, 346(4):1055–1077, 2003.
  * [18] Lisa J Kewley, MA Dopita, RS Sutherland, CA Heisler, and J Trevena.  Theoretical modeling of starburst galaxies.  The Astrophysical Journal, 556(1):121, 2001.
  * [19] Lisa J Kewley, Brent Groves, Guinevere Kauffmann, and Tim Heckman.  The host galaxies and classification of active galactic nuclei.  Monthly Notices of the Royal Astronomical Society, 372(3):961–976, 2006.
  * [20] P Kharb, S Subramanian, Sravani Vaddi, Mousumi Das, and Z Paragi.  Double-peaked emission lines due to a radio outflow in kissr 1219.  The Astrophysical Journal, 846(1):12, 2017.
  * [21] Henry A Kobulnicky, John M Dickey, Anneila I Sargent, David E Hogg, and Peter S Conti.  Aperture synthesis observations of molecular and atomic gas in the wolf-rayet starburst galaxy.  The Astronomical Journal, 110:116, 1995.
  * [22] Claus Leitherer, Daniel Schaerer, Jeffrey D Goldader, Rosa M González Delgado, Carmelle Robert, Denis Foo Kune, Duília F de Mello, Daniel Devost, and Timothy M Heckman.  Starburst99: synthesis models for galaxies with active star formation.  The Astrophysical Journal Supplement Series, 123(1):3, 1999.
  * [23] Roberto Maiolino, HR Russell, AC Fabian, Stefano Carniani, R Gallagher, S Cazzoli, S Arribas, F Belfiore, E Bellocchi, Luis Colina, et al.  Star formation inside a galactic outflow.  Nature, 544(7649):202–206, 2017.
  * [24] Christina M Manzano-King, Gabriela Canalizo, and Laura V Sales.  Agn-driven outflows in dwarf galaxies.  The Astrophysical Journal, 884(1):54, 2019.
  * [25] NL Martín-Hernández, D Schaerer, E Peeters, AGGM Tielens, and M Sauvage.  High spatial resolution mid-infrared spectroscopy of the starburst galaxies ngc 3256, ii zw 40 and henize 2–10.  Astronomy & Astrophysics, 455(3):853–870, 2006.
  * [26] DS Mathewson, MA Dopita, IR Tuohy, and VL Ford.  A new oxygen-rich supernova remnant in the large magellanic cloud.  The Astrophysical Journal, 242:L73–L76, 1980.
  * [27] Mallory Molina, Amy E Reines, Jenny E Greene, Jeremy Darling, and James J Condon.  Outflows, shocks, and coronal line emission in a radio-selected agn in a dwarf galaxy.  The Astrophysical Journal, 910(1):5, 2021.
  * [28] MA Nawaz, GV Bicknell, AY Wagner, RS Sutherland, and BR McNamara.  Jet–intracluster medium interaction in hydra a–ii. the effect of jet precession.  Monthly Notices of the Royal Astronomical Society, 458(1):802–815, 2016.
  * [29] Matthew Newville, Till Stensitzki, Daniel B Allen, Michal Rawlik, Antonino Ingargiola, and Andrew Nelson.  Lmfit: Non-linear least-square minimization and curve-fitting for python.  Astrophysics Source Code Library, pages ascl–1606, 2016.
  * [30] Dieu D Nguyen, Anil C Seth, Amy E Reines, Mark den Brok, David Sand, and Brian McLeod.  Extended structure and fate of the nucleus in henize 2-10.  The Astrophysical Journal, 794(1):34, 2014.
  * [31] Chris Nixon and Andrew King.  Do jets precess… or even move at all?  The Astrophysical Journal Letters, 765(1):L7, 2013.
  * [32] Donald E Osterbrock and Gary J Ferland.  Astrophysics Of Gas Nebulae and Active Galactic Nuclei.  University science books, 2006.
  * [33] Samantha J Penny, Karen L Masters, Rebecca Smethurst, Robert C Nichol, Coleman M Krawczyk, Dmitry Bizyaev, Olivia Greene, Charles Liu, Mariarosa Marinelli, Sandro B Rembold, et al.  Sdss-iv manga: evidence of the importance of agn feedback in low-mass galaxies.  Monthly Notices of the Royal Astronomical Society, 476(1):979–998, 2018.
  * [34] JE Pringle.  Self-induced warping of accretion discs.  Monthly Notices of the Royal Astronomical Society, 281(1):357–361, 1996.
  * [35] Amy E Reines, James J Condon, Jeremy Darling, and Jenny E Greene.  A new sample of (wandering) massive black holes in dwarf galaxies from high-resolution radio observations.  The Astrophysical Journal, 888(1):36, 2020.
  * [36] Amy E Reines and Adam T Deller.  Parsec-scale radio emission from the low-luminosity active galactic nucleus in the dwarf starburst galaxy henize 2-10.  The Astrophysical Journal Letters, 750(1):L24, 2012.
  * [37] Amy E Reines, Mark T Reynolds, Jon M Miller, Gregory R Sivakoff, Jenny E Greene, Ryan C Hickox, and Kelsey E Johnson.  Deep chandra observations of the compact starburst galaxy henize 2–10: X-rays from the massive black hole.  The Astrophysical Journal Letters, 830(2):L35, 2016.
  * [38] Amy E Reines, Gregory R Sivakoff, Kelsey E Johnson, and Crystal L Brogan.  An actively accreting massive black hole in the dwarf starburst galaxy henize 2-10.  Nature, 470(7332):66, 2011.
  * [39] Rogemar A Riffel.  Evidence for an accreting massive black hole in he 2–10 from adaptive optics integral field spectroscopy.  Monthly Notices of the Royal Astronomical Society, 494(2):2004–2011, 2020.
  * [40] F Santoro, C Tadhunter, D Baron, R Morganti, and J Holt.  Agn-driven outflows and the agn feedback efficiency in young radio galaxies.  Astronomy & Astrophysics, 644:A54, 2020.
  * [41] Joseph Silk.  Unleashing positive feedback: linking the rates of star formation, supermassive black hole accretion, and outflows in distant galaxies.  The Astrophysical Journal, 772(2):112, 2013.
  * [42] Joseph Silk and Colin Norman.  Global star formation revisited.  The Astrophysical Journal, 700(1):262, 2009.
  * [43] Jonathan R Trump, Christopher D Impey, Brandon C Kelly, Francesca Civano, Jared M Gabor, Aleksandar M Diamond-Stanic, Andrea Merloni, C Megan Urry, Heng Hao, Knud Jahnke, et al.  Accretion rate and the physical nature of unobscured active galaxies.  The Astrophysical Journal, 733(1):60, 2011.
  * [44] Jonathan R Trump, Tohru Nagao, Hiro Ikeda, Takashi Murayama, Christopher D Impey, John T Stocke, Francesca Civano, Martin Elvis, Knud Jahnke, Brandon C Kelly, et al.  Spectropolarimetric evidence for radiatively inefficient accretion in an optically dull active galaxy.  The Astrophysical Journal, 732(1):23, 2011.
  * [45] Sylvain Veilleux and Donald E Osterbrock.  Spectral classification of emission-line galaxies.  The Astrophysical Journal Supplement Series, 63:295–310, 1987.


# IMU-based online multi-lidar calibration without lidar odometry

Sandipan Das1,2, Bengt Boberg2 1 KTH EECS, Sweden<EMAIL_ADDRESS>
2 Scania, Sweden. {sandipan.das<EMAIL_ADDRESS>
When deploying autonomous systems that require several sensors for perception,
accurate and reliable extrinsic calibration is required. In this research, we
offer a reliable technique that can extrinsically calibrate numerous lidars in
the base frame of a moving vehicle without the use of odometry estimation or
fiducial markers. Our method is based on comparing the raw IMU signals between
a collocated IMU present with the lidar and the IMU measurements from the GNSS
system in the vehicle base frame. Additionally, based on our observability
criterion, we choose measurements that include the most mutual information
rather than comparing all comparable IMU readings. This enables us to locate
the measurements that are most useful for real-time calibration. Utilizing
data gathered from Scania test vehicles with various sensor setups, we have
successfully validated our methodology.

##  I Introduction

For safe navigation and redundancy, autonomous vehicles need numerous sensors
to generate $360^{\circ}$ sensing coverage. Therefore, it is crucial to
establish the precise mounting location for the vehicle’s sensors, a process
known as extrinsic calibration. Without the proper extrinsic calibration
parameters, it is impossible to fuse the sensing data into a single reference
frame. We concentrated on multi-lidar system extrinsic calibration in our
work. The same principle, though, might also be applied to other modalities.

There are two methods for extrinsic calibration: offline and online. Putting
fiducial markers in their shared field-of-view (FoV) of both sensors allows
the lidars’ characteristics to be matched offline [1, 2]. Online calibration
is carried out by comparing the estimated states of the sensor to those in the
relevant sensor frame [3, 4, 5, 6, 7]. Since the sensors may not have a common
field of view, it is preferable to execute online calibration because it
requires less engineering work than an offline method. In online calibration,
it is widely recognized that not all motion segments give the information
necessary for extrinsic calibration [8, 9, 10]. Hence, identifying the
degenerate motion segments helps to discard information that is not necessary
for extrinsic calibration computation.

We offer an alternative vehicle-run, an observability-aware online calibration
method for simultaneously calibrating multiple lidars using collocated IMU
(Inertial measurement unit) signals without the need for state estimation.
This method can be extended to other modalities seamlessly and does not
require common FoV in between the sensors. Additionally, we compensate for the
biases of the IMU measurements online for improving the calibration accuracy.

Figure 1: Illustration of the 2 lidars with their embedded IMUs positioned
around one data collection vehicle, when the sensors are calibrated. The
vehicle base frame $\mathtt{{B}}$, is located at the center of the rear axle.
The sensor frames of the lidars are:
$\mathtt{L}^{(\mathtt{F}_{\mathtt{L}}-\mathtt{Top})}$ and
$\mathtt{L}^{(\mathtt{F}_{\mathtt{R}}-\mathtt{Top})}$, whereas, the sensor
frames of the IMUs are: $\mathtt{I}^{(\mathtt{F}_{\mathtt{L}}-\mathtt{Top})}$
and $\mathtt{I}^{(\mathtt{F}_{\mathtt{R}}-\mathtt{Top})}$.
${(\mathtt{F}_{\mathtt{L}})}$: front-left and ${(\mathtt{F}_{\mathtt{R}})}$:
front-right.

###  I-A Motivation

As seen in the Fig. 1, the mobile platforms we employed for our studies have a
wide variety of sensors. Therefore, it is essential to confirm that the
sensors are correctly calibrated prior to any autonomous run. As a result, we
must calibrate every sensor online with respect to a common reference frame,
which is often found in the middle of the rear axle of the vehicle. The state
estimates for motion-based online calibration have intrinsic drifts under
aggressive maneuvers which are needed for online calibration to excite the
different degrees of freedom. It could therefore cause inaccuracies in the
calibration process.

This is why we propose calibrating the lidars using a comparison of the raw
IMU signals (with bias compensation) between the collocated IMU(s) with the
lidar(s), and the IMU signals from the GNSS system calibrated to the vehicle
reference frame. Our approach is based on comparing the raw measurements, as
opposed to comparing estimated states, which by their very nature contain
process noise. Furthermore, we solely compare the raw angular velocity data
for observability analysis.

###  I-B Contribution

The extensive literature on motion-based calibration serves as the inspiration
for our work. We would want to contribute the following:

  * •

An observability-aware extrinsic calibration algorithm to calibrate multiple
lidars online using collocated IMUs without the need of lidar odometry
estimation.

  * •

Online bias compensation of the IMU signals to improve the quality of the
extrinsic calibration algorithm.

  * •

Verification of our method based on data collected from Scania autonomous
vehicles with the sensor setup shown in Fig. 2, with FoV schematics similar to
Fig. 1.

##  II Related Work

Extrinsic sensor calibration is a well studied area. In our discussion we
briefly review the relevant literature and motivate our choice of method.

###  II-A Offline calibration

Offline extrinsic calibration is performed by matching known markers placed in
the common FoV of multiple sensors. To calibrate a multi-lidar system,
geometric elements including points, lines, and planes were retrieved and
matched in [1, 2].

###  II-B Online calibration

Online extrinsic calibration based on the Hand-Eye method [11, 12, 3] is an
extensively studied topic. The term “Hand-Eye” refers to early motion-based
calibration research that estimated the motion of the gripper (hand) and the
camera (eye) while constricting their poses with a fixed rigid body
transformation. Multi-camera extrinsic calibration with the Hand-Eye method
has been explored in works like Camodocal [4] and recently extended to multi-
lidar extrinsic calibration [13, 6, 14].

Instead of performing Hand-Eye calibration, in Kalibr [8] the authors
automatically detected sets of measurements from which they could identify an
observable parameter space and then performed a maximum likelihood estimate
(MLE) by minimizing the errors between landmark observations and their known
correspondences. They also discarded parameter updates for numerically
unobservant directions and degenerate scenarios. Kalibr was extended towards
MIMU calibration [15] based on matching poses fitted to a B-spline and further
constraining with known image landmarks. In MIMC-VINS [16], an efficient
multi-state constraint Kalman filter was used to jointly propagate all IMU
states while enforcing rigid body constraints between the IMUs during the
filter update stage, which produced the MIMU calibration.

We claim that because we have integrated IMUs in our setup, calibrating the
multi-IMU (MIMU) configuration will enable us to retrieve extrinsics from the
multi-lidar system, given the lidar-IMU extrinsics are known. We also selected
measurements that have enough signal excitation known as observability-aware
criteria. Nilsson et al. [17] provided housing with noncoplanar orientations
for observability analysis, in which their MIMU setup was placed and the IMU
measurements were matched using MLE. Since we wanted to calibrate the sensors
directly in the vehicle, our observability-aware criteria closely resemble
work by Jiajun et al. [7]. However, unlike prior work, we maximized mutual
information between the angular velocity signals to identify relevant motion
segments.

Unlike prior work, we used the fundamental tenet that a rigid body’s angular
velocity is constant around every point of the body as the foundation for our
MIMU calibration. We used MLE to match the raw angular velocity between the
IMU(s) embedded within the lidar(s) and GNSS system and recovered the
extrinsics in accordance with the signal-to-signal match principle [18]. We
also estimated the IMU biases online by capturing data in a special sequence,
removing the need for any exteroceptive modality.

##  III Problem Statement

###  III-A Sensor platform and reference frames

Figure 2: Reference frames conventions for our vehicle platform. The world
frame $\mathtt{W}$ is a fixed frame, while the base frame $\mathtt{{B}}$, as
shown in Fig. 1, is located at the rear axle center of the vehicle. Each
sensor unit contains the two optical frames $\mathtt{C}$, an IMU frame,
$\mathtt{I}$, and lidar frame $\mathtt{L}$.

The sensor platform with its corresponding reference frames is shown in Fig. 2
along with the illustrative sensor FoV of the bus in Fig. 1. Each of the
sensor housings contains a lidar with an embedded IMU and two cameras. Please
note the cameras are not used for this work. We used logs from a bus and a
truck with similar sensor housings for our experiments.

Now we describe the necessary notation and reference frames used in our
system. The vehicle base frame, $\mathtt{{B}}$ is located in the center of the
rear axle of the vehicle. Sensor readings from lidars and IMUs are represented
in their respective sensor frames as $\mathtt{L}^{(k)}$, and
$\mathtt{I}^{(k)}$ respectively. Here,
$k\in[\mathtt{F}_{\mathtt{L}}-\mathtt{Top},\mathtt{F}_{\mathtt{R}}-\mathtt{Top}]$
denotes the location of the sensor in the vehicle corresponding to front-left-
top and front-right-top respectively. The GNSS measurements are reported in a
world fixed frame, $\mathtt{W}$, and transformed to $\mathtt{{B}}$ frame by
performing a calibration routine outside the scope of this work. The GNSS
system also has an embedded IMU which is transformed to the base frame,
$\mathtt{{B}}$ as well. In our discussions the transformation matrix is
denoted as, $\mathbf{T}=\left[\begin{array}[]{ll}\mathbf{R}_{3\times
3}&\mathbf{t}_{3\times 1}\\\
\mathbf{0}^{\top}&1\end{array}\right]\in\mathrm{SE}(3)$ and
$\mathbf{R}\mathbf{R}^{T}=\mathbf{I}_{3\times 3}$, since the rotation matrix
is orthogonal.

###  III-B Problem formulation

Our primary goal is to estimate the extrinsic calibration of multiple lidar
sensors in base frame, $\mathtt{{B}}$ in real-time without the need of any
fiducial markers. We only use the raw IMU signals from $\mathtt{I}^{(k)}$ and
$\mathtt{{B}}$ frame (from GNSS system) and use the known extrinsics
$\mathbf{T}_{\mathtt{I}^{(k)}\mathtt{L}^{(k)}}$ from sensor supplier to
estimate $\mathbf{T}_{\mathtt{{B}}\mathtt{L}^{(k)}}$. Additionally, we
estimate the IMU biases online to denoise the raw IMU signals which improves
the calibration quality.

##  IV Methodology

###  IV-A Initialization

The IMU measurements are in their corresponding sensor frame. For gravity
alignment, we use the equations [19, Eq. 25, 26] to estimate roll and pitch
and obtain $\mathtt{R}_{\mathtt{W}\mathtt{{B}}}$ after collecting IMU data
when the vehicle is static for a few seconds.

###  IV-B IMU sensor model

We have considered a 6-DoF (degree of freedom) IMU such that it has a 3-axis
accelerometer and 3-axis gyroscope. The IMU sensor data in its corresponding
sensor frame can be represented as,

$\left[\begin{array}[]{c}\boldsymbol{\omega}_{\mathtt{I}_{t}}\\\
\mathbf{a}_{\mathtt{I}_{t}}\end{array}\right]=\left[\begin{array}[]{c}\boldsymbol{\hat{\boldsymbol{\omega}}_{\mathtt{I}_{t}}}\\\
\hat{\mathbf{a}}_{\mathtt{I}_{t}}-\mathbf{R}_{\mathtt{W}\mathtt{I}_{t}}^{T}\mathbf{g}_{\mathtt{W}}\end{array}\right]+\left[\begin{array}[]{l}\mathbf{b}_{\mathtt{I}_{t}}^{\omega}\\\
\mathbf{b}_{\mathtt{I}_{t}}^{a}\end{array}\right]+\left[\begin{array}[]{c}\mathbf{n}_{\mathtt{I}_{t}}^{\omega}\\\
\mathbf{n}_{\mathtt{I}_{t}}^{a}\end{array}\right].$ (1)

Here, $[\boldsymbol{\omega}_{\mathtt{I}_{t}}\in\mathbb{R}^{3\times
1},\mathbf{a}_{\mathtt{I}_{t}}\in\mathbb{R}^{3\times 1}]$ represent the
measured angular velocity and linear acceleration at timestamp $t$ and
$[\hat{\boldsymbol{\omega}}_{\mathtt{I}_{t}},\hat{\mathbf{a}}_{\mathtt{I}_{t}}]$
represent the latent ideal angular velocity and linear acceleration
respectively. $\mathbf{b}_{\mathtt{I}_{t}}^{\omega}\in\mathbb{R}^{3\times 1}$
and $\mathbf{b}_{\mathtt{I}_{t}}^{a}\in\mathbb{R}^{3\times 1}$ represent the
gyro and accelerometer biases which change with time and other factors like
temperature.
$\mathbf{n}_{\mathtt{I}_{t}}^{\omega}\sim\mathcal{N}(0,\boldsymbol{\Sigma}_{\omega})$
and
$\mathbf{n}_{\mathtt{I}_{t}}^{a}\sim\mathcal{N}(0,\boldsymbol{\Sigma}_{a})$
are the additive zero-mean white Gaussian noises for gyroscope and
accelerometer with covariance
$\boldsymbol{\Sigma}_{\omega}\in\mathbb{R}^{3\times 3}$ and
$\boldsymbol{\Sigma}_{a}\in\mathbb{R}^{3\times 3}$ respectively.
$\mathbf{g}_{\mathtt{W}}\in\mathbb{R}^{3\times 1}$, represents the gravity
vector in $\mathtt{W}$ frame and $\mathbf{R}_{\mathtt{{B}}\mathtt{W}}$
represent the gravity alignment rotation matrix.

###  IV-C IMU bias characterization

In our method, since we are not using an external odometry source for bias
compensation we collect the data sequences with periods of rest in between.
This limits the bias covariance growth as they converge faster after rest
detection.

####  IV-C1 IMU state propagation

The IMU dynamical model based on angular velocity in quaternion form is shown
as,

$\mathbf{q}_{\mathtt{{B}}\mathtt{I}_{t}}=\left[\cos\left(\frac{\Delta
t}{2}\left\|\boldsymbol{\omega}_{\mathtt{I}_{t}}\right\|\right)\mathbf{I}_{4}+\frac{\sin\left(\frac{\Delta
t}{2}\left\|\boldsymbol{\omega}_{\mathtt{I}_{t}}\right\|\right)}{\left\|\boldsymbol{\omega}_{\mathtt{I}_{t}}\right\|}\boldsymbol{\Omega}_{\mathtt{I}_{t}}\right]\mathbf{q}_{\mathtt{{B}}\mathtt{I}_{t-1}}$
(2)

and,

$\boldsymbol{\Omega}_{\mathtt{I}_{t}}=\left[\begin{array}[]{cccc}0&{\left[\boldsymbol{\omega}_{\mathtt{I}_{t}}\right]_{z}}&-\left[\boldsymbol{\omega}_{\mathtt{I}_{t}}\right]_{y}&{\left[\boldsymbol{\omega}_{\mathtt{I}_{t}}\right]_{x}}\\\
-\left[\boldsymbol{\omega}_{\mathtt{I}_{t}}\right]_{z}&0&{\left[\boldsymbol{\omega}_{\mathtt{I}_{t}}\right]_{x}}&{\left[\boldsymbol{\omega}_{\mathtt{I}_{t}}\right]_{y}}\\\
{\left[\boldsymbol{\omega}_{\mathtt{I}_{t}}\right]_{y}}&-\left[\boldsymbol{\omega}_{\mathtt{I}_{t}}\right]_{x}&0&{\left[\boldsymbol{\omega}_{\mathtt{I}_{t}}\right]_{z}}\\\
-\left[\boldsymbol{\omega}_{\mathtt{I}_{t}}\right]_{x}&-\left[\boldsymbol{\omega}_{\mathtt{I}_{t}}\right]_{y}&-\left[\boldsymbol{\omega}_{\mathtt{I}_{t}}\right]_{z}&0\end{array}\right].$
(3)

$\Delta t$ denotes the sampling period of the IMU data. We used the Madgwick
filter [20] to estimate the refined rotation which minimizes the difference
between the measured acceleration and the aligned gravity vector using a
gradient descent algorithm as,

$\mathbf{q}_{\mathtt{{B}}\mathtt{I}_{t}}^{*}=\operatorname*{arg\,min}_{\mathbf{q}_{\mathtt{{B}}\mathtt{I}_{t}}\in\mathbb{R}^{4\times
1}}\left(\begin{array}[]{l}\tilde{\mathbf{q}}_{\mathtt{{B}}\mathtt{I}_{t}}\otimes{\mathbf{R}_{\mathtt{{B}}\mathtt{{W}}}^{-1}}\mathbf{g}_{\mathtt{W}}\otimes\mathbf{q}_{\mathtt{{B}}\mathtt{I}_{t}}-\hat{\mathbf{a}}_{\mathtt{{B}}\mathtt{I}_{t}}\\\
\end{array}\right),$ (4)

where, $\tilde{\mathbf{q}}$ represents the quaternion conjugate and $\otimes$
is the quaternion product operator. Note that we use the compensated
accelerometer signals after bias compensation which is discussed in IV-C3.

####  IV-C2 Rest detection

Rest is detected if the difference between the norm of the aligned gravity
vector and the acceleration vector is less than a predefined threshold, $\tau$
for at least more than 2 seconds. Thus,

$\mathtt{Rest\
detected}=\|\begin{array}[]{l}\tilde{\mathbf{q}}^{*}_{\mathtt{{B}}\mathtt{I}_{t}}\otimes\mathbf{g}_{\mathtt{{B}}}\otimes\mathbf{q}^{*}_{\mathtt{{B}}\mathtt{I}_{t}}-\hat{\mathbf{a}}_{\mathtt{{B}}\mathtt{I}_{t}}\end{array}\|_{2}\leq\tau$
(5)

####  IV-C3 Accelerometer bias estimation

The accelerometer bias is estimated by computing the mean of the signal when
stationary and the standard deviation gives us the co-variance of the white
Gaussian noise. We recompute these parameters whenever the vehicle is
stationary and keep it unchanged until the next rest period. If the rest
period duration is $N$ seconds then,

$\displaystyle\mathbf{b}^{a}$ $\displaystyle=\frac{\Delta
t}{N}\sum_{i=1}^{N}(\mathbf{a}_{\mathtt{I}_{i}}-{\mathbf{R}_{\mathtt{{B}}\mathtt{I}_{t}}^{*}}{\mathbf{R}_{\mathtt{{B}}\mathtt{{W}}}^{-1}}\mathbf{g}_{\mathtt{W}}),$
(6) $\displaystyle\boldsymbol{\Sigma}_{a}$
$\displaystyle=\mathtt{diag}\left(\frac{\Delta
t}{N}\sum_{i=1}^{N}|\mathbf{a}_{\mathtt{I}_{i}}-\mathbf{b}^{a}|^{2}\right),$
$\displaystyle\implies\hat{\mathbf{a}}_{\mathtt{I}_{t}}$
$\displaystyle={\mathbf{a}}_{\mathtt{I}_{t}}+{\mathbf{R}_{\mathtt{{B}}\mathtt{I}_{t}}^{*}}{\mathbf{R}_{\mathtt{{B}}\mathtt{{W}}}^{-1}}\mathbf{g}_{\mathtt{W}}-\mathbf{b}^{a}-\mathcal{N}(0,\boldsymbol{\Sigma}_{a}).$

####  IV-C4 Gyro bias estimation

The initialization of the gyro bias is done in a similar fashion as we did for
the accelerometer and recomputed whenever the rest period is detected. Thus,

$\displaystyle\text{System State},\hat{\mathbf{x}}_{0}$
$\displaystyle=\mathbf{b}^{\omega}_{\mathtt{I}_{0}}=\frac{\Delta
t}{N}\sum_{i=1}^{N}(\boldsymbol{\omega}_{\mathtt{I}_{i}})$ (7)
$\displaystyle\text{Measurement Noise},\mathbf{W}_{0}$
$\displaystyle=\mathcal{N}(0,\boldsymbol{\Sigma}_{\omega})$
$\displaystyle\text{Process Noise},\mathbf{Q}$
$\displaystyle=(0.1^{\circ}/\mathtt{sec})^{2}{\mathbf{I}}_{3\times 3}$
$\displaystyle\text{Initial Covariance},\mathbf{P}_{0}$
$\displaystyle=(0.3^{\circ}/\mathtt{sec})^{2}{\mathbf{I}}_{3\times 3}.$

After that we use a Kalman filter [21] to track the gyroscope bias as a state.
The system is modeled by,

$\displaystyle\mathbf{x}_{t}$
$\displaystyle=\hat{\mathbf{x}}_{t-1}+\mathcal{N}(0,\mathbf{Q})$ (8)
$\displaystyle\mathbf{y}_{t}$
$\displaystyle=\mathbf{R}_{\mathtt{{B}}\mathtt{I}_{t}}^{*}\mathbf{x}_{t}+\mathbf{W}_{t}.$

The standard Kalman filter update equations become,

$\displaystyle\mathbf{P}_{t}$ $\displaystyle=\mathbf{P}_{t-1}+\mathbf{Q}$ (9)
$\displaystyle\mathbf{K}_{t}$
$\displaystyle=\mathbf{P}_{t}{\mathbf{R}_{\mathtt{{B}}\mathtt{I}_{t}}^{*T}}(\boldsymbol{\Sigma}_{\omega}+\mathbf{R}_{\mathtt{{W}}\mathtt{I}_{t}}^{*}\mathbf{P}_{t}{\mathbf{R}_{\mathtt{{B}}\mathtt{I}_{t}}^{*T}})^{-1}$
$\displaystyle\hat{\mathbf{x}}_{t}$
$\displaystyle=\hat{\mathbf{x}}_{t-1}+\mathbf{K}_{t}(\mathbf{y}_{t}-\mathbf{R}_{\mathtt{{B}}\mathtt{I}_{t}}^{*}\hat{\mathbf{x}}_{t-1})$
$\displaystyle\mathbf{P}_{t}$
$\displaystyle=\mathbf{P}_{t}-\mathbf{K}_{t}\mathbf{R}_{\mathtt{{B}}\mathtt{I}_{t}}^{*}\mathbf{P}_{t}.$

Thus after estimating the biases we can rearrange eq. 1 to compute
$\hat{\boldsymbol{\omega}}_{\mathtt{I}_{t}}$.

###  IV-D IMU-based lidar calibration

####  IV-D1 Rotation estimation

The extrinsics of $\mathtt{L}^{(k)}$ wrt $\mathtt{{B}}$ frame has a strong
correlation with the extrinsics of $\mathtt{I}^{(k)}$ wrt $\mathtt{{B}}$ frame
as they are collocated with known extrinsics from sensor supplier. A rigid
body has the same angular velocity at all places. The angular velocities of
all the IMU sensors in our sensor arrangement must be equal if their
orientations are the same because they are all securely affixed to the
vehicle. Let, [$\hat{\boldsymbol{\omega}}_{\mathtt{I}_{i}^{(k)}}$,
$\hat{\mathbf{a}}_{\mathtt{I}_{i}^{(k)}}$] and
[$\hat{\boldsymbol{\omega}}_{\mathtt{{B}}_{i}}$,
$\hat{\mathbf{a}}_{\mathtt{{B}}_{i}}$] be the estimated angular velocity and
linear acceleration of the $k^{th}$ IMU and the base IMU at timestamp $t_{i}$
after bias compensation. By removing the $k^{th}$ superscript for brevity, our
optimization problem becomes,

$\displaystyle\mathbf{R}_{\mathtt{{B}}\mathtt{I}}^{\star}$
$\displaystyle=\operatorname*{arg\,min}_{\mathbf{R}_{\mathtt{{B}}\mathtt{{I}}}}\sum_{i=1}^{N}\|\mathbf{R}_{\mathtt{{B}}\mathtt{{I}}}\hat{\boldsymbol{\omega}}_{\mathtt{{B}}_{i}}-\hat{\boldsymbol{\omega}}_{\mathtt{I}_{i}}\|^{2}_{\boldsymbol{\sum_{i}}},$
(10) $\displaystyle\textit{s.t.}\ \
\mathbf{R}_{\mathtt{{B}}\mathtt{I}}\mathbf{R}_{\mathtt{{B}}\mathtt{I}}^{T}=\mathbf{I}_{3},$

which, can be simply solved with Horn alignment [22] or Kabsch alignment [23].

####  IV-D2 Translation estimation

For the acceleration components we compensate for the Coriolis forces as
illustrated in Fig. 3 and equate the translation components as,

$\displaystyle{(\mathbf{R}_{\mathtt{{B}}\mathtt{I}}^{*})}^{-1}\hat{\mathbf{a}}_{\mathtt{I}}$
$\displaystyle=\hat{\mathbf{a}}_{\mathtt{{B}}}+\underbrace{\hat{\boldsymbol{\omega}}_{\mathtt{{B}}}\times\left(\hat{\boldsymbol{\omega}}_{\mathtt{{B}}}\times\mathbf{t}_{\mathtt{{B}}\mathtt{I}}\right)}_{\text{Centrifugal
force
}}+\underbrace{\dot{\hat{\boldsymbol{\omega}}}_{\mathtt{{B}}}\times\mathbf{t}_{\mathtt{{B}}\mathtt{I}}}_{\text{Euler
force }}$ (11)
$\displaystyle=\hat{\mathbf{a}}_{\mathtt{{B}}}+[\hat{\boldsymbol{\omega}}_{\mathtt{{B}}}]^{2}_{\times}\mathbf{t}_{\mathtt{{B}}\mathtt{I}}+[\dot{\hat{\boldsymbol{\omega}}}_{\mathtt{{B}}}]_{\times}\mathbf{t}_{\mathtt{{B}}\mathtt{I}},$

where, $[\boldsymbol{.}]_{\times}$ is a skew-symmetric matrix and
$[\boldsymbol{a}]_{\times}\boldsymbol{b}=\boldsymbol{a}\times\boldsymbol{b}$.

Figure 3: Illustration of IMU signal transformation between $\mathtt{{B}}$ and
$\mathtt{I}^{(\mathtt{F}_{\mathtt{L}}-\mathtt{Top})}$ frame.

Since, we already computed $\mathbf{R}^{\star}_{\mathtt{{B}}\mathtt{I}}$, our
optimization problem for the translation component becomes,

$\displaystyle\mathbf{t}_{\mathtt{{B}}\mathtt{I}}^{\star}=\operatorname*{arg\,min}_{\mathbf{t}_{\mathtt{{B}}\mathtt{{I}}}}\sum_{i=1}^{N}\leavevmode\resizebox{294.85841pt}{}{$\|\underbrace{(\left[\boldsymbol{\hat{\omega}_{\mathtt{{B}}_{i}}}\right]^{2}_{\times}+[\dot{\hat{\boldsymbol{\omega}}}_{\mathtt{{B}}_{i}}]_{\times})}_{\mathbf{A}_{i}}\mathbf{t}_{\mathtt{{B}}\mathtt{I}}-\underbrace{(\mathbf{R}^{-1\star}_{\mathtt{{B}}\mathtt{{I}}}\hat{\boldsymbol{a}}_{\mathtt{{I}}_{i}}-\hat{\boldsymbol{a}}_{\mathtt{{B}}_{i}})}_{\mathbf{B}_{i}}\|^{2},$}$
(12)

where, $\mathbf{A}=[\mathbf{A}_{0}\ \mathbf{A}_{1}\ \dots\
\mathbf{A}_{n}]^{T}$, $\mathbf{B}=[\mathbf{B}_{0}\ \mathbf{B}_{1}\ \dots\
\mathbf{B}_{n}]^{T}$ and $\mathbf{x}=\mathbf{t}_{\mathtt{{B}}\mathtt{I}}$
$\forall i\in[1,N]$. This is a system of linear equations of the form,
$\mathbf{A}{\mathbf{x}}=\mathbf{B}$, which can be solved using a least-square
approach as,

$\mathbf{x}^{\star}=\operatorname*{arg\,min}_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{B}\|^{2}_{\boldsymbol{\Sigma}},$
(13)

where, ${\boldsymbol{\Sigma}}$ denotes the co-variance of the residual. The
main approach [24] used to tackle these problems involves solving a series of
approximations to the original problem repeatedly by linearizing as,
$F(\mathbf{x}+\Delta\mathbf{x})\approx
F(\mathbf{x})+\mathbf{J}(\mathbf{x})\Delta\mathbf{x}$, where, $J$ being the
jacobian of $F(\mathbf{x})$. Thus, $\hat{\mathbf{x}}$ is updated in the
current iteration as,
$\hat{\mathbf{x}}\leftarrow\hat{\mathbf{x}}\boxplus\mathbf{\Delta x}$, where
$\boxplus$ is an addition operator in the manifold and the problem becomes,

$\displaystyle\mathbf{\Delta
x}^{\star}=\operatorname*{arg\,min}_{\mathbf{\Delta
x}}\frac{1}{2}\|(\mathbf{A}\hat{\mathbf{x}}-\mathbf{B})+\mathbf{A}^{T}\mathbf{\Delta
x}\|^{2}_{\boldsymbol{\Sigma}}.$ (14)

The optimal solution is given by,

$\underbrace{\left(\mathbf{J}^{\top}\boldsymbol{\Sigma}^{-1}\mathbf{J}\right)}_{\text{Fisher
information matrix}}\mathbf{\Delta
x}=-\mathbf{J}^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{A}\hat{\mathbf{x}}-\mathbf{B}).$
(15)

In practice, since we receive reliable initial estimate for the translation
component from the CAD parameters, we search for $\mathbf{x}^{\star}$ only in
a local neighborhood within reasonable bounds. Thus, we solve a bounded
variable least squares problem as,

$\mathbf{x}^{\star}=\operatorname*{arg\,min}_{\mathbf{L}\leq\mathbf{x}\leq\mathbf{U}}\|\mathbf{A}\mathbf{x}-\mathbf{B}\|^{2}_{\boldsymbol{\Sigma}},$
(16)

where, $\mathbf{L}$ and $\mathbf{U}$ are the lower and upper bounds of
$\mathbf{x}$ respectively. Thus the solution space in Eq. 14 is modified with
an additional constraint as,

$\displaystyle\mathbf{\Delta
x}^{\star}=\operatorname*{arg\,min}_{\mathbf{\Delta
x}}\frac{1}{2}\|(\mathbf{A}\hat{\mathbf{x}}-\mathbf{B})+\mathbf{A}^{T}\mathbf{\Delta
x}\|^{2}_{\boldsymbol{\Sigma}},$ (17) $\displaystyle\textit{s.t.}\
\mathbf{L}\leq\hat{\mathbf{x}}+\mathbf{\Delta x}\leq\mathbf{U}.$

###  IV-E Observability analysis

All motions do not excite enough degrees of freedom to allow calibration. As a
result, it’s important to identify the motion segments for information-aware
calibration updates. We capture this information by comparing the angular
velocity between $\mathtt{I}^{(k)}$ and $\mathtt{{B}}$ frame based on Eq.10
which can also be solved iteratively with the update step as,

$\displaystyle\underbrace{\left(\sum_{i}\mathbf{J}_{i}^{\top}\Sigma_{i}^{-1}\mathbf{J}_{i}\right)}_{\text{Fisher
information
matrix}}\mathbf{\Delta}{\mathbf{\mathbf{R}}}_{\mathtt{{B}}\mathtt{{I}}}$
$\displaystyle\stackrel{{\scriptstyle
Eq.~{}\ref{eq:bvls_iterative_solution}}}{{=}}-\sum_{i}\mathbf{J}_{i}^{\top}\Sigma_{i}^{-1}(\hat{\mathbf{\mathbf{R}}}_{\mathtt{{B}}\mathtt{{I}}}\boldsymbol{\omega}_{\mathtt{{B}}_{i}}-\boldsymbol{\omega}_{\mathtt{I}_{i}})$
(18)

where,
${\boldsymbol{\Sigma}_{i}}=\operatorname{cov}(\hat{\mathbf{\mathbf{R}}}_{\mathtt{{B}}\mathtt{{I}}}\boldsymbol{\omega}_{\mathtt{{B}}_{i}}-\boldsymbol{\omega}_{\mathtt{I}_{i}})$
and $\mathbf{J}_{i}=\boldsymbol{\omega}_{\mathtt{{B}}_{i}}^{T}$. The Fisher
information matrix, $\mathcal{I}_{N\times N}$ captures all the information
contained in the measurements. We do the processing after arranging the data
in a batch size of $N$. We perform a Singular Value Decomposition of
$\mathcal{I}_{N\times N}$ for each batch as:

$\mathcal{I}_{N\times N}=\textbf{USU}^{T},$ (19)

where, $\textbf{U}=[\textbf{u}_{1},\textbf{u}_{2}\dots,\textbf{u}_{N}]$ and
$\textbf{S}=\text{diag}(\sigma_{1},\sigma_{2}\dots,\sigma_{N})$ is a diagonal
matrix of singular values in decreasing order. Information about the data in
the batch is indicated by the value of the minimal singular value. If the
minimal singular value exceeds a certain threshold (design decision), we may
say that there are sufficient excitations in the batch of data to allow for
extrinsics computation and hence chosen for calibration.

##  V Experimental Results

To demonstrate the online calibration performance, we conducted several
experiments with real-world data collected using two vehicles with the sensor
setup illustrated in Fig. 2.

###  V-A Dataset

We use data from a GNSS receiver for the GT (ground truth) poses. We analyze
our results on 2 different collected test sequences which are described in
detail in Table I, covering different driving scenarios. In Seq-1, we
performed aggressive motion to excite the possible degrees of freedom.
However, in Seq-2 we recorded data in a normal driving scenario in regular
traffic conditions. In Seq-1, there were no periods of rest whereas in Seq-2
there were occasional periods of rest due to traffic signals.

Dataset collection details for the experiments
---
Data | Vehicle | Length ($\mathrm{K}\mathrm{m}$) | Duration ($\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{s}$) | Stop period
Seq-1† | Truck | 0.224 | 81 | No
Seq-2† | Bus | 1.088 | 283 | Yes
†$\mathtt{F}_{\mathtt{L}}-\mathtt{Top},\mathtt{F}_{\mathtt{R}}-\mathtt{Top}$
mounted lidars were used for the experiments.

TABLE I:

(a) Seq-1: Odometry estimation in aggressive driving scenario.

(b) Seq-2: Odometry estimation in non-aggressive driving scenario.

Figure 4: Odometry estimation of $\mathtt{F}_{\mathtt{L}}-\mathtt{Top}$ and
$\mathtt{F}_{\mathtt{R}}-\mathtt{Top}$ lidars using using Fast-lio2 [25] for
both the sequences. The trajectory alignment is done using Kabsch alignment
[23].

We evaluated the performance of a state-of-the-art lidar odometry estimator
(Fast-lio2 [25]) on both the motion sequences as seen in Fig. 4 and observe
that the estimator is sensitive to the aggressive motion sequence, Seq-1. The
RMSE of the absolute pose error (APE) for Seq-1 is 2.2926$\mathrm{m}\text{/}$
and 2.1063$\mathrm{m}\text{/}$ for $\mathtt{F}_{\mathtt{L}}-\mathtt{Top}$ and
$\mathtt{F}_{\mathtt{R}}-\mathtt{Top}$ lidars respectively, which is not
suitable for motion-based calibration using pose alignment strategies. As seen
in Fig. 4(a), the Kabsch alignment [23] algorithm does not lead us to correct
trajectory alignment. For Seq-2, the APE is less than 0.5 $\mathrm{m}\text{/}$
for both lidars when using Fast-lio2. However, this motion sequence doesn’t
have enough excitation in all the different degrees of freedom and thus is not
a great choice for a calibration motion sequence using our proposed method. We
did not run the state estimator for the whole sequence as not enough
excitation for the different degrees of freedom was present in the data.

###  V-B Bias estimation results

To improve our calibration performance, which relies on the signal-to-signal
match policy we compensate for the biases in both the accelerometer and
angular velocity signals. The accelerometer biases are recomputed whenever a
period of rest is detected. The angular velocity biases are tracked online
based on the Kalman filter. The update step is based on the estimated
orientation using the Madgwick filter that uses the bias-compensated
accelerometer signals for gravity alignment. Thus whenever a period of rest is
detected the covariances of the angular velocity biases converge.

Figure 5: Bias estimation for the $\mathtt{F}_{\mathtt{L}}-\mathtt{Top}$
angular velocity signals for Seq-1. Figure 6: Bias estimation for the
$\mathtt{F}_{\mathtt{L}}-\mathtt{Top}$ angular velocity signals for Seq-2.

In Seq-1, we didn’t have any periods of rest. Thus the angular velocity bias
covariances do not converge as seen in Fig. 5. Whereas, in Fig. 6, we see that
the angular velocity bias covariances converge after the detection of the rest
period for a few seconds at the beginning of Seq-2.

###  V-C Online calibration results

GT calibration of the lidars is obtained offline by refining the vehicle CAD
parameters. The refinement process analyzes known static feature positions
around the vehicle in the world frame, $\mathtt{W}$, and matches the
corresponding detected features from the lidar perspective. After that, we
obtain the extrinsics of the lidars to the base frame, $\mathtt{{B}}$ and
compute the transformation matrix,
$\mathbf{T}_{\mathtt{G}\mathtt{T}}=\mathbf{T}^{-1}_{\mathtt{{B}}\mathtt{L}^{(\mathtt{F}_{\mathtt{L}}-\mathtt{Top})}}\mathbf{T}_{\mathtt{{B}}\mathtt{L}^{(\mathtt{F}_{\mathtt{R}}-{\mathtt{Top}})}}$,
between the 2 lidars. We compare our results to this matrix by comparing the
translation and rotation errors as:

$\displaystyle\Delta t$
$\displaystyle=\frac{1}{3}\sqrt{\|\widehat{\mathbf{t}}-\mathbf{t}_{\operatorname{GT}}\|^{2}_{F}},$
(20) $\displaystyle\Delta R$
$\displaystyle=\frac{180}{\pi}\operatorname{cos}^{-1}\left[\frac{1}{2}\left(\operatorname{Tr}(\widehat{\mathbf{R}}^{-1}\mathbf{R}_{\operatorname{GT}})-1\right)\right],$
(21)

where, $\Delta R$ is the rotation along the principle eigenvector of
$(\widehat{\mathbf{R}}^{-1}\mathbf{R}_{\operatorname{GT}})$.

Figure 7: Raw IMU signals from $\mathtt{F}_{\mathtt{L}}-\mathtt{Top}$ and
$\mathtt{F}_{\mathtt{R}}-\mathtt{Top}$ IMUs for Seq-1. Figure 8: Calibrated
IMU signals from $\mathtt{F}_{\mathtt{L}}-\mathtt{Top}$ and
$\mathtt{F}_{\mathtt{R}}-\mathtt{Top}$ IMUs with bias compensation for Seq-1.

For Seq-1, we obtained $\Delta t=0.3387$\mathrm{m}\text{/}$$ and $\Delta
R=0.4577^{\circ}$ and the original IMU signals, as well as the calibrated IMU
signals between the 2 collocated IMUs, can be observed in Fig. 7 and Fig. 8
respectively. For Seq-2, the results were, $\Delta
t=0.4431$\mathrm{m}\text{/}$$ and $\Delta R=4.7938^{\circ}$. For both the
sequences, the translation components were obtained by solving a bound
variable least square optimization problem as defined in Eq.16 with bounds of
$\pm 0.3$\mathrm{m}\text{/}$$. The raw and calibrated signals for Seq-2 are
additionally shown in our complimentary video. As expected, in Seq-1, due to
the presence of rich motion sequence in terms of excitation of different
degrees of freedom our signal matching algorithm produces better results than
Seq-2.

###  V-D Observability analysis

For observability analysis we extract the information matrix by comparing the
angular velocity signals between the collocated IMU, $\mathtt{I}^{(k)}$, and
the $\mathtt{{B}}$ frame by dividing the data into equal segments of 10 sec
each. We analyze the SVD of the Fisher information matrix from Eq.18 and
select the IMU data in the segment for calibration if the minimum singular
value is greater than a threshold. For our experiments, we set the threshold
as $5^{-10}$.

(a) Selected poses with observability criteria for Seq-1 shown with random
colors.

(b) The angular velocity signals in different motion segments in a 10 sec
window.

Figure 9: Seq-1: The observability criteria selects all the motion segments as
there were enough excitation in the angular velocity signals in all the motion
segments.

As seen in Fig. 9(b), we selected all segments (highlighted in green) as there
was enough excitation in the angular velocity during all the time in Seq-1.
The corresponding selected trajectory segments are shown in Fig. 9(a).
Similarly for Seq-2, we can see that in Fig. 10(a), relevant poses for
calibration are selected only when there are significant turns in the
maneuvers as observed in Fig. 10(b), meaning there is enough excitation for
different degrees of freedom.

(a) Selected poses with observability criteria for Seq-2 shown with random
colors.

(b) We observe that in straight motion segments there are not enough
excitation in the angular velocity for Seq-2. Hence, the straight motion
segments are not used from observability perspective and thus ignored for
calibration.

Figure 10: Seq-2: The observability criteria selects motion segments when
there were enough excitation in the angular velocity signals, particularly
during the turns.

##  VI Conclusion

We show that it is possible to estimate the extrinsic calibration of multiple
lidars based on matching collocated IMU signals only. Unlike other methods
which use odometry estimates for matching poses, this is a lightweight method
that relies on raw signal matching. The observability-aware module selects
signals that carry enough excitation for different degrees of freedom required
for calibration. Our method provides comparable performance to GT extrinsics
parameters when sufficient signal excitations are present. This work can be
easily extended to other modalities like cameras and radars. For further
improvement, we can also match features in between the corresponding sensors
if they have a common FoV.

##  VII Acknowledgements

This research has been jointly funded by the Swedish Foundation for Strategic
Research (SSF) and Scania. The research was also affiliated with Wallenberg
AI, Autonomous Systems and Software Program (WASP).

## References

  * [1] M. He, H. Zhao, F. Davoine, J. Cui, and H. Zha, “Pairwise lidar calibration using multi-type 3d geometric features in natural scene,” in _IEEE International Conference on Intelligent Robots and Systems_ , 2013, pp. 1828–1835.
  * [2] D.-G. Choi, Y. Bok, J.-S. Kim, and I. S. Kweon, “Extrinsic calibration of 2-d lidars using two orthogonal planes,” _IEEE Transactions on Robotics_ , vol. 32, no. 1, pp. 83–98, 2016.
  * [3] J. Brookshire and S. Teller, “Extrinsic calibration from per-sensor egomotion,” _Robotics: Science and Systems VIII_ , pp. 504–512, 2013.
  * [4] L. Heng, B. Li, and M. Pollefeys, “Camodocal: Automatic intrinsic and extrinsic calibration of a rig with multiple generic cameras and odometry,” in _IEEE International Conference on Intelligent Robots and Systems_. IEEE, 2013, pp. 1793–1800.
  * [5] P. Furgale, J. Rehder, and R. Siegwart, “Unified temporal and spatial calibration for multi-sensor systems,” in _IEEE International Conference on Intelligent Robots and Systems_ , 2013, pp. 1280–1286.
  * [6] J. Jiao, H. Ye, Y. Zhu, and M. Liu, “Robust odometry and mapping for multi-lidar systems with online extrinsic calibration,” _IEEE Transactions on Robotics_ , vol. 38, no. 1, pp. 351–371, 2021.
  * [7] J. Lv, X. Zuo, K. Hu, J. Xu, G. Huang, and Y. Liu, “Observability-aware intrinsic and extrinsic calibration of lidar-imu systems,” _IEEE Transactions on Robotics_ , 2022.
  * [8] J. Maye, P. Furgale, and R. Siegwart, “Self-supervised calibration for robotic systems,” in _IEEE Intelligent Vehicles Symposium_ , 2013, pp. 473–480.
  * [9] K. Hausman, J. Preiss, G. S. Sukhatme, and S. Weiss, “Observability-aware trajectory optimization for self-calibration with application to uavs,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 3, pp. 1770–1777, 2017\.
  * [10] T. Schneider, M. Li, C. Cadena, J. Nieto, and R. Siegwart, “Observability-aware self-calibration of visual and inertial sensors for ego-motion estimation,” _IEEE Sensors Journal_ , vol. 19, no. 10, pp. 3846–3860, 2019.
  * [11] R. Tsai and R. Lenz, “A new technique for fully autonomous and efficient 3d robotics hand/eye calibration,” _IEEE Transactions on Robotics and Automation_ , vol. 5, no. 3, pp. 345–358, 1989.
  * [12] R. Horaud and F. Dornaika, “Hand-eye calibration,” _The international journal of robotics research_ , vol. 14, no. 3, pp. 195–210, 1995.
  * [13] J. Jiao, Y. Yu, Q. Liao, H. Ye, R. Fan, and M. Liu, “Automatic calibration of multiple 3d lidars in urban environments,” in _IEEE International Conference on Intelligent Robots and Systems_ , 2019, pp. 15–20.
  * [14] S. Das, N. Mahabadi, A. Djikic, C. Nassir, S. Chatterjee, and M. Fallon, “Extrinsic calibration and verification of multiple non-overlapping field of view lidar sensors,” in _IEEE International Conference on Robotics and Automation_ , 2022, pp. 919–925.
  * [15] J. Rehder, J. Nikolic, T. Schneider, T. Hinzmann, and R. Siegwart, “Extending kalibr: Calibrating the extrinsics of multiple imus and of individual axes,” in _IEEE International Conference on Robotics and Automation_. IEEE, 2016, pp. 4304–4311.
  * [16] K. Eckenhoff, P. Geneva, and G. Huang, “Mimc-vins: A versatile and resilient multi-imu multi-camera visual-inertial navigation system,” _IEEE Transactions on Robotics_ , vol. 37, no. 5, pp. 1360–1380, 2021.
  * [17] J.-O. Nilsson, I. Skog, and P. Händel, “Aligning the forces—eliminating the misalignments in imu arrays,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 63, no. 10, pp. 2498–2500, 2014.
  * [18] M. E. Pittelkau, “Calibration and attitude determination with redundant inertial measurement units,” _Journal of Guidance, Control, and Dynamics_ , vol. 28, no. 4, pp. 743–752, 2005.
  * [19] M. Pedley, “Tilt sensing using a three-axis accelerometer,” _Freescale semiconductor application note_ , vol. 1, pp. 2012–2013, 2013.
  * [20] S. O. H. Madgwick, A. J. L. Harrison, and R. Vaidyanathan, “Estimation of imu and marg orientation using a gradient descent algorithm,” in _2011 IEEE International Conference on Rehabilitation Robotics_ , 2011, pp. 1–7.
  * [21] R. E. Kalman, “A new approach to linear filtering and prediction problems,” _Transactions of the ASME–Journal of Basic Engineering_ , vol. 82, no. Series D, pp. 35–45, 1960.
  * [22] B. K. Horn, “Closed-form solution of absolute orientation using unit quaternions,” _Journal of the Optical Society of America_ , vol. 4, no. 4, pp. 629–642, 1987.
  * [23] W. Kabsch, “A discussion of the solution for the best rotation to relate two sets of vectors,” _Crystallographica: Crystal Physics, Diffraction, Theoretical and General Crystallography_ , pp. 827–828, 1978.
  * [24] J. Nocedal and S. J. Wright, _Numerical optimization_. Springer, 1999.
  * [25] W. Xu, Y. Cai, D. He, J. Lin, and F. Zhang, “FAST-LIO2: Fast direct lidar-inertial odometry,” _IEEE Transactions on Robotics_ , vol. 38, no. 4, pp. 2053–2073, 2022.