# Field linkage and magnetic helicity density

K. Lund1, M. Jardine1, A. J. B. Russell2, J.-F. Donati3, R. Fares4, C. P.
Folsom3, S. V. Jeffers5, S. C. Marsden6, J. Morin7, P. Petit4 and V. See8
1SUPA, School of Physics and Astronomy, University of St Andrews, North Haugh,
St Andrews, KY16 9SS, UK
2School of Science & Engineering, University of Dundee, Nethergate, Dundee DD1
4HN, UK
3IRAP, Université de Toulouse, CNRS, UPS, CNES, 14 Avenue Edouard Belin,
31400, Toulouse, France
4Physics Department, United Arab Emirates University, P.O. Box 15551, Al-Ain,
United Arab Emirates
5Institut für Astrophysik, Universität Göttingen, Friedrich-Hund-Platz 1,
D-37077 Göttingen, Germany
6University of Southern Queensland, Centre for Astrophysics, Toowoomba, QLD,
4350, Australia
7LUPM, Université de Montpellier, CNRS, Place Eugène Bataillon, F-34095
Montpellier, France
8University of Exeter, Department of Physics & Astronomy, Stocker Road, Devon,
Exeter, EX4 4QL, UK

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The helicity of a magnetic field is a fundamental property that is conserved
in ideal MHD. It can be explored in the stellar context by mapping large-scale
magnetic fields across stellar surfaces using Zeeman-Doppler imaging. A recent
study of 51 stars in the mass range 0.1-1.34 M⊙ showed that the photospheric
magnetic helicity density follows a single power law when plotted against the
toroidal field energy, but splits into two branches when plotted against the
poloidal field energy. These two branches divide stars above and below $\sim$
0.5 M⊙. We present here a novel method of visualising the helicity density in
terms of the linkage of the toroidal and poloidal fields that are mapped
across the stellar surface. This approach allows us to classify the field
linkages that provide the helicity density for stars of different masses and
rotation rates. We find that stars on the lower-mass branch tend to have
toroidal fields that are non-axisymmetric and so link through regions of
positive and negative poloidal field. A lower-mass star may have the same
helicity density as a higher-mass star, despite having a stronger poloidal
field. Lower-mass stars are therefore less efficient at generating large-scale
helicity.

###### keywords:

stars: magnetic field – methods: analytical

††pubyear: 2020††pagerange: Field linkage and magnetic helicity density–Field
linkage and magnetic helicity density

##  1 Introduction

Magnetic helicity is a fundamental property of magnetic fields that measures
the amount of linkage and twist of field lines within a given volume. Since it
is exactly conserved in ideal MHD and highly conserved for high magnetic
Reynolds numbers in general (Woltjer, 1958; Taylor, 1974), helicity is an
important factor when attempting to understand how magnetic fields are
generated and evolve (e.g. Brandenburg & Subramanian, 2005; Chatterjee et al.,
2011; Pipin et al., 2019). Until recently, this could only be measured for the
Sun (e.g. reviews by Démoulin, 2007; Démoulin & Pariat, 2009). We can,
however, now map all three components of the large-scale magnetic field at the
surfaces of stars using the spectropolarimetric technique of Zeeman-Doppler
imaging (Semel, 1989).

These magnetic field maps now exist for a large enough sample of stars that
trends with stellar mass and rotation period have become apparent (Donati &
Landstreet, 2009). In particular, it appears that magnetic fields show
different strengths and topologies in the mass ranges above and below $\sim$
0.5 M⊙, which is believed to correspond to the onset of the transition from
partially to fully convective interiors. Rapidly-rotating stars in the mass
range above $\sim$ 0.5 M⊙ tend to have fields that are predominantly toroidal
(Donati et al., 2008a). The stronger the toroidal field, the more likely it is
to be axisymmetric (See et al., 2015). In the mass range below $\sim$ 0.5 M⊙,
stars show predominantly axisymmetric poloidal fields. For the lowest masses,
however, a bimodal behaviour is found, such that stars may have strong,
predominantly axisymmetric poloidal fields, or much weaker, non-axisymmetric
poloidal fields (Donati et al., 2008a; Morin et al., 2008b; Donati &
Landstreet, 2009; Morin et al., 2010).

This difference in magnetic fields in stars that are partially or fully
convective is also apparent in their photospheric helicity densities. Using
observations of 51 stars, Lund et al. (2020) found that the helicity density
scales with the toroidal energy according to $|\langle{h\,}\rangle|$ $\propto$
$\langle{\rm{B_{tor}}^{2}\,\rangle}^{0.86\,\pm\,0.04}$. The scaling with the
poloidal energy is more complex, however, revealing two groups with different
behaviours. Specifically, stars less massive than $\sim$ 0.5 M⊙ appear to have
an excess of poloidal energy when compared to more massive stars with similar
helicity densities. It appears that stars with different internal structures
and different total magnetic energies may nonetheless generate magnetic fields
with the same helicity density at their surfaces. The aim of this paper is to
explore the nature of this division and the types of flux linkage that support
the measured helicity densities. In order to do that, we have developed a
novel method of visualising the linkages of different field components across
the surfaces of stars.

##  2 Methods

###  2.1 Poloidal and toroidal magnetic field components

For the purposes of this paper the stellar magnetic fields discussed will be
decomposed into their poloidal and toroidal components:
$\boldsymbol{B}=\boldsymbol{B}_{\rm{pol}}+\boldsymbol{B}_{\rm{tor}}$. The
poloidal and toroidal fields can be expressed in a general form as (see
Appendix III of Chandrasekhar, 1961):

$\boldsymbol{B}_{\rm{pol}}=\curl{[\curl{[\Phi\hat{r}]}]},$ (1)
$\boldsymbol{B}_{\rm{tor}}=\curl{[{\Psi\hat{r}}]}.$ (2)

In a spherical coordinate system111We use a right-handed spherical coordinate
system where a positive radial field component points out of the star, the
$\theta$ component is positive pointing from North to South and the $\phi$
component is positive in the clockwise direction as viewed from the South
pole. the scalars $\Phi$ and $\Psi$ take the form:

$\Phi=S(r)c_{lm}P_{lm}e^{im\phi},$ (3) $\Psi=T(r)c_{lm}P_{lm}e^{im\phi}.$ (4)

$P_{lm}\equiv P_{lm}(\cos{\theta})$ is the associated Legendre polynomial of
mode $l$ and order $m$ and

$c_{lm}\equiv\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}$ (5)

is a normalisation constant. $S(r)$ and $T(r)$ are functions describing the
radial behaviour of the magnetic field components. Determining the complete
form of these functions from observations of stellar magnetic fields is
impossible, however values can be obtained at stellar surfaces
($r=R_{\star}$).

Figure 1: The cartoon shows four different combinations of symmetries
(axisymmetric and non-axisymmetric relative to the rotation axis) of the
poloidal (red arrows) and toroidal (blue arrows) fields. Poloidal field lines
that link with the toroidal field are represented by solid lines, the ones
that do not are dashed.

The Zeeman-Doppler imaging technique (Semel, 1989) describes the large-scale
(low $l$ modes) magnetic fields at the surfaces of stars in terms of
$\alpha_{lm}$, $\beta_{lm}$, and $\gamma_{lm}$ coefficients (e.g. Donati et
al., 2006; Vidotto, 2016):

$\begin{split}\boldsymbol{B}_{\rm{pol}}(\theta,\phi)&=\sum_{lm}\alpha_{lm}c_{lm}P_{lm}e^{im\phi}\hat{r}\\\
&+\sum_{lm}\frac{\beta_{lm}}{(l+1)}c_{lm}\derivative{P_{lm}}{\theta}e^{im\phi}\hat{\theta}\\\
&+\sum_{lm}\frac{\beta_{lm}im}{(l+1)\sin\theta}c_{lm}P_{lm}e^{im\phi}\hat{\phi},\end{split}$
(6)
$\begin{split}\boldsymbol{B}_{\rm{tor}}(\theta,\phi)&=\sum_{lm}\frac{\gamma_{lm}im}{(l+1)\sin\theta}c_{lm}P_{lm}e^{im\phi}\hat{\theta}\\\
&-\sum_{lm}\frac{\gamma_{lm}}{(l+1)}c_{lm}\derivative{P_{lm}}{\theta}e^{im\phi}\hat{\phi}.\end{split}$
(7)

These expressions are consistent with the general form of the poloidal and
toroidal fields evaluated at the stellar surface when

$S(R_{\star})=\frac{\alpha_{lm}R_{\star}^{2}}{l(l+1)},\quad{\derivative{S(r)}{r}}\Big{\rvert}_{r=R_{\star}}=\frac{\beta_{lm}R_{\star}}{(l+1)}$
(8)

and

$T(R_{\star})=\frac{\gamma_{lm}R_{\star}}{(l+1)}.$ (9)

Given the surface magnetic field, magnetic energies are estimated by
calculating the mean squared magnetic flux density $\langle{B^{2}}\rangle$.
For instance, in the case of the poloidal energy222When calculating the mean
squared magnetic flux density we integrate over a full sphere, dividing by the
solid angle of $\Omega=$ 4$\pi$.:
$\langle{B^{2}}_{\rm{pol}}\rangle=\frac{1}{\Omega}\int\boldsymbol{B}_{\rm{pol}}\cdot\boldsymbol{B}_{\rm{pol}}\mathrm{d}\Omega$.
Accordingly, the fraction of axisymmetric poloidal magnetic field energy is
given by333We define axisymmetric as $m=0$.:
$\langle{B^{2}_{\rm{pol},m=0}}\rangle/\langle{B^{2}}_{\rm{pol}}\rangle$. The
toroidal energy and axisymmetry fraction are calculated analogously.

###  2.2 Magnetic helicity density

Magnetic helicity can be defined as
$H=\int\boldsymbol{A}\cdot\boldsymbol{B}\,\rm{dV}$ (Woltjer, 1958), where
$\boldsymbol{A}$ is a vector potential corresponding to the magnetic field
$\boldsymbol{B}$. As magnetic helicity is a quantity measuring the linkage of
fields within a volume, our surface magnetic fields limit us to evaluating the
magnetic helicity density $h=\boldsymbol{A}\cdot\boldsymbol{B}$. The
separation of the magnetic field into its poloidal and toroidal components is
particularly useful in this regard. It dispenses with the need to invoke a
gauge (Berger & Hornig, 2018) since the usual gauge field (the corresponding
potential field with the same boundary flux) has zero helicity. In addition,
in a spherical coordinate system, toroidal field lines lie purely on spherical
surfaces while poloidal field lines pass through these surfaces. This makes
visualising the linkage of field lines straightforward. Fig. 1 illustrates how
poloidal fields lines (shown red) that pass through the stellar surface may
link through loops of toroidal field (shown blue) that lie on the stellar
surface. It is notable that in these examples, only some fraction of the
poloidal field links with the toroidal field line shown.

Interpreting magnetic helicity as the linking of poloidal and toroidal fields
(Berger, 1985; Berger & Hornig, 2018) allows the helicity density to be
calculated for any stellar magnetic map given only its $\alpha_{lm}$ and
$\gamma_{lm}$ coefficients and the stellar radius $R_{\star}$ (Lund et al.,
2020):

$\begin{split}h(\theta,\phi)&=\real(\sum_{lm}\sum_{l^{\prime}m^{\prime}}\frac{\alpha_{lm}\gamma_{l^{\prime}m^{\prime}}R_{\star}}{(l^{\prime}+1)l(l+1)}c_{lm}c_{l^{\prime}m^{\prime}}e^{i\phi(m+m^{\prime})}\\\
&\bigg{(}P_{lm}P_{l^{\prime}m^{\prime}}\bigg{(}l(l+1)-\frac{mm^{\prime}}{\sin^{2}\theta}\bigg{)}+\derivative{P_{lm}}{\theta}\derivative{P_{l^{\prime}m^{\prime}}}{\theta}\bigg{)}\bigg{)}.\end{split}$
(10)

The magnetic helicity density, as expressed in Eq. 10, depends on the
$\alpha_{lm}$ and $\gamma_{lm}$ coefficients, but not the $\beta_{lm}$
coefficients found in the $\theta$ and $\phi$ components of the poloidal
field. This is because only the radial component of the poloidal field
($\boldsymbol{B}_{\rm{pol},r}$) passes through the spherical surfaces
containing the toroidal field.

When comparing different magnetic maps, e.g. for different stars, it is often
useful to summarise the overall helicity with a single number. For this
purpose, we consider an average value across the hemisphere facing the
observer. We note that typically only one hemisphere is fully visible as part
of the star never comes into view. Furthermore, we take the absolute value of
the averaged helicity density as we are interested in comparing magnitudes,
not signs. This absolute average helicity density (|$\langle{h}\rangle$|) will
for simplicity be referred to as the “helicity density” for the remainder of
this paper.

Figure 2: (a) Contours of the function $\Psi$ lying on the surface of the star
are shown in blue. (b) This shows an enlarged version, illustrating that
$\boldsymbol{B}_{\rm{tor}}=\gradient{\Psi}\crossproduct\hat{r}$.

It is possible to visualise the linkage of poloidal and toroidal fields that
results in helicity density at the stellar surface through maps showing the
strength of $\boldsymbol{B}_{\rm{pol},r}$ (calculated from the $\alpha_{lm}$
according to Eq. 6) with the field lines of $\boldsymbol{B}_{\rm{tor}}$
superimposed. Expanding Eq. 2 as

$\boldsymbol{B}_{\rm{tor}}=\gradient{\Psi}\crossproduct\hat{r}$ (11)

shows that the contours of $\Psi$ correspond to the field lines of
$\boldsymbol{B}_{\rm{tor}}$ (see Fig. 2). In particular, at the stellar
surface,

$\Psi=\frac{\gamma_{lm}R_{\star}}{(l+1)}c_{lm}P_{lm}e^{im\phi},$ (12)

from Eqs. 4 and 9.

##  3 Stellar sample

Figure 3: Absolute helicity density averaged across a single hemisphere versus
the mean squared poloidal magnetic flux density ($l\leq$ 4). The colour of the
symbols correspond to Rossby number and the shape of the symbols splits the
sample into two mass groups; circles represent $M_{\star}$ > 0.5 M⊙ and
diamonds represent $M_{\star}\leq$ 0.5 M⊙. When there are multiple
measurements for the same star these are connected by lines. The thick grey
lines show the best fit of |$\langle{h\,}\rangle$| =
$\langle{B^{2}_{\rm{pol}}\,}\rangle^{\alpha}10^{\beta}$ for stars in the two
mass groups; $\alpha=1.04\,\pm\,0.05$, $\beta=10.18\,\pm\,0.13$ for
$M_{\star}$ > 0.5 M⊙ (circles) and $\alpha=0.88\,\pm\,0.15$,
$\beta=9.57\,\pm\,0.74$ for $M_{\star}\leq$ 0.5 M⊙ (diamonds). The stars
outlined are shown in Fig. 4.

Our sample of stellar magnetic maps are all created using Zeeman-Doppler
imaging. They describe the magnetic fields of 51 different stars, 15 of which
are represented by multiple maps. The stars range in spectral type from F to
M, and in mass from 0.1-1.34 M⊙. Details of each star/map are provided in
Table 1, along with the calculated helicity densities and magnetic energy
components. For the sake of a fair comparison every magnetic map is evaluated
to the same resolution, which means every calculation is performed up to the
same $l$-mode ($l\leq 4$) even when higher modes are available (Lund et al.,
2020).

##  4 The role of axisymmetry in helicity density

Figure 4: A visualisation of the linkage of the dipole ($l=1$) poloidal and
toroidal field components of GJ 182, WX UMa (2008), GJ 49 and GJ 1245B (2006).
The colour shows the strength of the radial magnetic (poloidal) field, and the
black contours represent the toroidal magnetic field lines. The heavy black
contour separates regions of positive (solid) and negative (dashed) toroidal
field. These examples correspond roughly to the four classes shown in Fig. 1.

The very well-defined dependence of helicity density on the toroidal field of
$|\langle{h\,}\rangle|$ $\propto$
$\langle{\rm{B_{tor}}^{2}\,\rangle}^{0.86\,\pm\,0.04}$ was shown in Lund et
al. (2020). One enduring puzzle, however, is that the dependence on the
poloidal field revealed two branches, as shown in Fig. 3. The higher-mass
branch ($M_{\star}>$ 0.5 M⊙) follows |$\langle{h\,}\rangle$|
$=\langle{B^{2}_{\rm{pol}}\,}\rangle^{1.04\,\pm\,0.05}10^{10.18\,\pm\,0.13}$
and the lower-mass branch ($M_{\star}\leq$ 0.5 M⊙) follows
|$\langle{h\,}\rangle$|
$=\langle{B^{2}_{\rm{pol}}\,}\rangle^{0.88\,\pm\,0.15}10^{9.57\,\pm\,0.74}$.
When fitting these power laws the sample is split specifically at 0.5 M⊙
because a number of magnetic properties, including helicity density, have been
shown to change behaviour across this value (Donati et al., 2008b; Morin et
al., 2008b; Morin et al., 2010; See et al., 2015; Lund et al., 2020). Fig. 3
shows that the lowest mass stars have higher poloidal energies than higher-
mass stars with the same helicity density. The lowest mass stars in our sample
also typically have the lowest Rossby numbers, as indicated by the colours in
the plot.

It appears from Fig. 3 that the lower-mass fully-convective stars have excess
poloidal field that does not contribute to their helicity density. In order to
explore the distribution of the poloidal and toroidal fields on these two
branches and to determine their contribution to the helicity density, we plot
maps showing how their poloidal and toroidal fields link. As an example, Fig.
4 presents maps for GJ 182, WX UMa, GJ 49 and GJ 1245B. These stars are
highlighted in Fig. 3 and represent two pairs of stars with approximately the
same helicity density, and thus similar toroidal energies. Each pair consists
of one star from the higher-mass branch and one star from the lower-mass
branch, and the two pairs are distinguished by the magnitude of their helicity
densities (GJ 182 and WX UMa (top row) have higher helicity densities than GJ
49 and GJ 1245B). To illustrate field linkages at the largest scale the maps
show the dipole ($l=1$) mode. We note that these four stars fall into the four
categories shown in Fig. 1.

Figure 5: Axisymmetric toroidal energy as a fraction of total toroidal energy
versus the ratio of the mean squared toroidal to poloidal magnetic flux
densities ($l\leq$ 4). The symbols are the same as in Fig. 3. Figure 6: The
helicity energy fraction ($l\leq 4$), defined as $\tilde{h}\equiv$
|$\langle{h\,}\rangle$|/$\langle{R_{\star}B^{2}\,}\rangle$, versus stellar
mass. The symbols are the same as in Fig. 3.

Figure 7: Helicity density of the large-scale magnetic fields ($l\leq 4$) of
our stellar sample shown according to stellar mass and rotation period. The
dashed lines split the sample at 0.5 M⊙. The size of the symbol indicates
relative strength of the helicity density. The colour corresponds to the
fractional axisymmetry of the toroidal (top) and poloidal (bottom) magnetic
energies. The stars outlined are shown in Fig. 4.

The clearest trend to emerge is in the toroidal field. It is notable that the
stars on the higher-mass branch, GJ 182 and GJ 49 (left column of Fig. 4),
both have fairly axisymmetric toroidal fields, whilst the stars on the lower
mass branch, WX UMa and GJ 1245B, have non-axisymmetric toroidal fields. We
can quantify this trend by plotting the fraction of the toroidal field energy
that is held in axisymmetric modes as a function of the ratio of toroidal to
poloidal energy (see Fig. 5). Lower-mass stars tend to have toroidal fields
that are non-axisymmetric and magnetic energy budgets that are dominated by
the poloidal field.

Lund et al. (2020) showed that stars with a similar helicity density have
toroidal fields of similar strengths. As can be seen from Fig. 4, however,
these toroidal fields may have different symmetries if the stars lie on two
different branches. The fact that these different symmetries are able to
produce the same helicity density is because the poloidal field, although
typically stronger on the lower-mass branch, can also be of different symmetry
to the toroidal field. In both rows of Fig. 4, the transition from the higher-
mass branch to the lower-mass branch is accompanied by an increase in the
strength of the poloidal field. At the same time, for the lower mass stars,
toroidal field lines enclose regions of both positive and negative poloidal
field such that the total linked flux is relatively small. The difference in
magnetic field topology roughly cancels the difference in field strength, such
that the helicity is very similar within each pair, although it is arrived at
very differently.

To separate the effects of field strength and geometry Fig. 6 plots stellar
mass against the “helicity energy fraction”; $\tilde{h}\equiv$
|$\langle{h\,}\rangle$|/$\langle{R_{\star}B^{2}\,}\rangle$. Dividing the
helicity density by the mean squared magnetic flux density and stellar radius
results in a dimensionless helicity density normalised by magnetic field
strength. Fig. 6 shows an even spread of low helicity energy fractions across
the entire range of stellar masses, however only the higher mass stars
($M_{\star}>$ 0.5 M⊙) exceed a fraction of 0.35. This illustrates that even
though the lower-mass stars are among those whose dynamos are most efficient
at injecting magnetic energy into the largest spatial scales that ZDI is able
to detect and map (e.g. Morin et al., 2008b), they are apparently less
efficient at generating helicity at these largest scales. This is most likely
because of the inefficient linking between their poloidal and toroidal
magnetic field components.

Inefficient linking can be arrived at in more than one way. For the higher-
helicity pair in Fig. 4 (GJ 182, WX UMa), the lower-mass star has a strongly
axisymmetric poloidal field, and since the toroidal field is non-axisymmetric,
this combination produces inefficient linking. In the case of the lower-
helicity pair (GJ 49 and GJ 1245B), the lower-mass star’s poloidal and
toroidal fields are both strongly non-axisymmetric; nonetheless, they still
offset each another by approximately 90 degrees, which again gives inefficient
linking.

We can place these trends in a broader context by showing how the axisymmetry
and helicity density varies across the stellar mass-rotation period plane.
This is shown in Fig. 7 which also shows separately the variation of the
toroidal and poloidal axisymmetry fractions. High-helicity stars exist in both
mass ranges. The decline in the axisymmetry of the toroidal fields with
decreasing mass is apparent, but the trends in poloidal axisymmetry are more
complex. The reasons for this, and the potential role of magnetic cycles, are
not yet clear.

Fig. 7 also demonstrates very clearly the differences in the strength and
structure of the magnetic field that is possible for stars in the bimodal
regime, which includes the stars of lowest mass and shortest period. Two of
our example stars (WX Uma and GJ 1245B) lie in this regime. They have similar
masses (0.1 and 0.12 M⊙ respectively) and rotation periods (0.78 and 0.71
days, respectively) but their magnetic fields and the helicity densities they
support are quite different.

##  5 Discussion and conclusions

Helicity measures the linkage within a field. By studying the linkage of the
poloidal and toroidal components of stellar magnetic fields, we can learn
about the underlying dynamo processes generating the field, and thereby the
form of the magnetic cycles that might take place and the evolution of stellar
fields as stars spin down over their main sequence lifetimes. This study
reveals that stars in different mass ranges, which may be either fully or
partially convective, generate their helicity through different forms of
toroidal/poloidal field linkage.

For partially convective stars (those that lie on the higher-mass branch in
Fig. 3) the toroidal fields are mainly axisymmetric. An increase in rotation
rate (or a decrease in Rossby number) generally leads to increased helicity
density. For fully-convective stars (those that lie on the lower-mass branch
in Fig. 3) the toroidal fields are mainly non-axisymmetric. For the lowest-
mass stars, a bimodal behaviour is present in the form and strength of the
magnetic field (Morin et al., 2011) which may be weak with a non-axisymmetric
poloidal field, or strong with an axisymmetric poloidal field. These two types
correspond to low and high helicity densities respectively. This may have
implications for the possibility that this represents a bimodality in the
dynamo operating in this regime (Morin et al., 2011).

Stars can evolve from one mass-branch to the other if their field structure
changes, for instance as a result of their internal structure changing from
mainly convective to mainly radiative at a very young age, or if they
transition from one bimodal dynamo mode to the other. Furthermore, stars can
evolve along each branch as their rotation rates decay with age. Their field
linkages can also evolve on much shorter timescales due to magnetic cycles. It
is notable in Fig. 5 that where there are multiple observations of a star,
taken at different times, these typically follow the trend that an increase in
the ratio of $B^{2}_{\rm tor}$ to $B^{2}_{\rm pol}$ leads to an increase in
the axisymmetry of the toroidal field. A similar behaviour is seen in Fig. 3
where for each star with multiple observations, these all lie within the
scatter about the best-fit line.

It is not clear what causes the non-axisymmetry of the toroidal fields in the
lowest-mass stars. Their deep-seated convection may produce bipoles that
emerge through the stellar surface with randomised axial tilts, leading to a
lack of axisymmetry in the toroidal field. Their low surface differential
rotation may also reduce the shearing of bipoles and hence the diffusive
cancellation of poloidal field that results. Both of these processes, however,
occur at length scales well below what can be resolved by these Zeeman-Doppler
field measurements. The field characteristics that are most robustly recovered
by Zeeman-Doppler imaging (field axisymmetry and the ratio of poloidal to
toroidal field (Lehmann et al., 2019)) are nonetheless the very ones that
underpin the helicity density.

In summary, we find that lowest-mass stars tend to be inefficient at
generating helicity on the largest scales, given their magnetic energy. The
helicity density at a stellar surface depends not only on the stellar radius
and the strengths of the individual poloidal and toroidal field components
(see Eq. 10) but also on their spatial distribution relative to each other.
The fraction of the poloidal flux that links with the toroidal flux is
maximised when the axes of symmetry of the two fields align perfectly. The
bottom row of Figure 1 illustrates two ways in which such an alignment is
possible: where both fields are axisymmetric (left) and where both are non-
axisymmetric (right). Conversely, if the poloidal and toroidal axes of
symmetry are orthogonal to each other there is no linkage, and consequently no
helicity. Different orientations with different amounts of field linking can
nonetheless result in the same helicity density due to the dependence on field
strengths. In short, to achieve a full understanding of the source of the
helicity density at a stellar surface, it is not enough simply to look at, for
example, the radius and field strength, it is also necessary to produce a
surface map showing the field linkage. It is only through a combination of all
these components that a clear picture can be formed.

Table 1: Our stellar sample, with the four stars we focus on in this paper highlighted in bold (GJ 182, GJ 49, GJ 1245B and WX UMa). From left to right the columns show: star name, mass, radius, rotation period,Rossby number, absolute helicity density averaged across the visible hemisphere, $\langle{B^{2}_{\rm{pol}}\,}\rangle$, poloidal axisymmetric magnetic energy as a fraction of poloidal energy, $\langle{B^{2}_{\rm{tor}}\,}\rangle$, toroidal axisymmetric magnetic energy as a fraction of toroidal energy, $l_{\rm{max}}$ and observation epoch. The helicity density and the energies are all calculated for $l\leq 4$. References for the stellar parameters are given in the last column, where references to the papers where the magnetic maps were published are in italic. A more comprehensive table of parameters for these stars can be found in Vidotto et al. (2014). Star ID |  $M_{\star}$ |  $R_{\star}$ |  $P_{\rm{rot}}$ |  $R_{o}$ |  |$\langle{h\,}\rangle$| |  $\langle{B^{2}_{\rm{pol}}\,}\rangle$ |  Pol |  $\langle{B^{2}_{\rm{tor}}\,}\rangle$ |  Tor |  $l_{\rm{max}}$ |  Obs. | Ref.
---|---|---|---|---|---|---|---|---|---|---|---|---
|  (M⊙) |  (R⊙) |  (d) |  |  (Mx${}^{2}\rm{cm}^{-3}$) |  (G2) |  Axi |  (G2) |  Axi |  |  epoch |
Solar-like stars
HD 3651 |  0.88 |  0.88 |  44.0 |  1.916 |  1.82E+11 |  1.49E+01 |  0.87 |  4.96E-01 |  0.98 |  10 |  – | 1, 2, 3
HD 9986 |  1.02 |  1.04 |  22.4 |  1.621 |  7.16E+09 |  4.71E-01 |  0.50 |  3.43E-02 |  0.94 |  10 |  – | 1, 2, 3
HD 10476 |  0.82 |  0.82 |  35.2 |  0.576 |  6.77E+09 |  6.23E+00 |  0.00 |  5.69E-01 |  0.40 |  10 |  – | 1, 4, 2, 3
HD 20630 |  1.03 |  0.95 |  9.00 |  0.593 |  2.09E+13 |  2.61E+02 |  0.32 |  4.56E+02 |  0.90 |  10 |  Oct 2012 | 1, 5, 2, 6
HD 22049 |  0.86 |  0.77 |  11.76 |  0.366 |  7.19E+11 |  1.16E+02 |  0.70 |  3.35E+00 |  0.23 |  10 |  - | 1, 5, 2, 7
HD 39587 |  1.03 |  1.05 |  5.136 |  0.295 |  1.95E+12 |  1.45E+02 |  0.07 |  2.04E+02 |  0.84 |  10 |  – | 1, 5, 2, 3
HD 56124 |  1.03 |  1.01 |  20.7 |  1.307 |  2.43E+11 |  5.22E+00 |  0.90 |  9.32E-01 |  0.91 |  10 |  – | 1, 2, 3
HD 72905 |  1 |  1 |  5.227 |  0.272 |  1.39E+13 |  1.43E+02 |  0.15 |  8.70E+02 |  0.97 |  10 |  – | 1, 5, 2, 3
HD 73350 |  1.04 |  0.98 |  12.3 |  0.777 |  3.54E+11 |  8.99E+01 |  0.00 |  8.91E+01 |  0.90 |  10 |  – | 1, 8, 2, 3
HD 75332 |  1.21 |  1.24 |  3.870 |  >1.105 |  9.88E+11 |  4.58E+01 |  0.80 |  3.96E+00 |  0.39 |  15 |  – | 1, 5, 2, 3
HD 78366 |  1.34 |  1.03 |  11.4 |  >2.781 |  2.39E+12 |  2.11E+02 |  0.94 |  8.36E+00 |  0.29 |  10 |  2008 | 1, 2, 9
… |  … |  … |  … |  … |  4.27E+11 |  4.83E+01 |  0.06 |  7.91E+00 |  0.54 |  … |  2010 | …
… |  … |  … |  … |  … |  4.29E+11 |  2.18E+01 |  0.77 |  2.25E+00 |  0.30 |  … |  2011 | …
HD 101501 |  0.85 |  0.9 |  17.04 |  0.663 |  3.02E+12 |  1.28E+02 |  0.26 |  5.68E+01 |  0.80 |  10 |  – | 1, 5, 2, 3
HD 131156A |  0.90 |  0.80 |  6.25 |  0.256 |  2.11E+13 |  1.48E+03 |  0.30 |  4.12E+03 |  0.93 |  10 |  Aug 2007 | 10, 2, 11
… |  … |  … |  … |  … |  5.24E+12 |  5.58E+02 |  0.59 |  3.00E+02 |  0.62 |  … |  Feb 2008 | …
… |  … |  … |  … |  … |  6.64E+12 |  4.91E+02 |  0.07 |  4.54E+02 |  0.81 |  … |  Jun 2009 | …
… |  … |  … |  … |  … |  4.05E+12 |  3.89E+02 |  0.08 |  1.91E+02 |  0.37 |  … |  Jan 2010 | …
… |  … |  … |  … |  … |  1.40E+13 |  2.73E+02 |  0.34 |  2.03E+02 |  0.95 |  … |  Jun 2010 | …
… |  … |  … |  … |  … |  6.60E+13 |  7.18E+02 |  0.73 |  1.12E+03 |  0.96 |  … |  Aug 2010 | …
… |  … |  … |  … |  … |  2.91E+13 |  5.25E+02 |  0.27 |  1.77E+03 |  0.97 |  … |  Jan 2011 | …
HD 131156B |  0.66 |  0.55 |  11.1 |  0.611 |  5.59E+12 |  2.60E+02 |  0.25 |  1.25E+02 |  0.81 |  10 |  – | 10, 2, 3
HD 146233 |  0.98 |  1.02 |  22.7 |  1.324 |  8.87E+08 |  1.90E+00 |  0.09 |  1.48E-02 |  0.05 |  10 |  Aug 2007 | 12, 2, 8
HD 166435 |  1.04 |  0.99 |  4.2 |  0.259 |  5.12E+12 |  2.53E+02 |  0.50 |  1.70E+02 |  0.79 |  10 |  – | 1, 2, 3
HD 175726 |  1.06 |  1.06 |  4.0 |  0.272 |  1.71E+12 |  6.65E+01 |  0.18 |  2.65E+01 |  0.80 |  10 |  – | 1, 13, 2, 3
HD 190771 |  0.96 |  0.98 |  8.80 |  0.453 |  5.02E+12 |  8.50E+01 |  0.35 |  1.53E+02 |  0.98 |  10 |  2007 | 9, 2, 8
HD 201091A |  0.66 |  0.62 |  34.1 |  0.786 |  8.29E+10 |  3.02E+01 |  0.03 |  3.57E+00 |  0.42 |  10 |  – | 1, 5, 2, 14
HD 206860 |  1.1 |  1.04 |  4.6 |  0.388 |  2.77E+13 |  3.03E+02 |  0.49 |  3.71E+02 |  0.94 |  10 |  – | 1, 2, 15
Young suns
AB Dor |  0.76 |  1.00 |  0.5 |  0.026 |  2.87E+14 |  3.34E+04 |  0.14 |  1.21E+04 |  0.56 |  25 |  Dec 2001 | 16, 17, 2, 18
… |  … |  … |  … |  … |  6.80E+13 |  2.97E+04 |  0.09 |  5.92E+03 |  0.50 |  … |  Dec 2002 | ….
BD-16351 |  0.9 |  0.88 |  3.21 |  0.14${}^{+0.01}_{-0.02}$ |  5.04E+13 |  2.94E+03 |  0.04 |  2.45E+03 |  0.88 |  15 |  Sep 2012 | 19
HII 296 |  0.9 |  0.93 |  2.61 |  0.13${}^{+0.01}_{-0.01}$ |  4.74E+13 |  3.36E+03 |  0.38 |  3.96E+02 |  0.53 |  15 |  Oct 2009 | 19
HII 739 |  1.15 |  1.07 |  1.58 |  0.25${}^{+0.01}_{-0.08}$ |  1.28E+12 |  1.78E+02 |  0.28 |  1.09E+02 |  0.65 |  15 |  Oct 2009 | 19
HIP 12545 |  0.95 |  1.07 |  4.83 |  0.14${}^{+0.02}_{-0.02}$ |  8.28E+14 |  1.36E+04 |  0.49 |  7.63E+03 |  0.85 |  15 |  Sep 2012 | 19
HIP 76768 |  0.80 |  0.85 |  3.70 |  0.09${}^{+0.03}_{-0.02}$ |  4.63E+14 |  5.23E+03 |  0.75 |  9.10E+03 |  0.95 |  15 |  May 2013 | 19
TYC 0486-4943-1 |  0.75 |  0.69 |  3.75 |  0.13${}^{+0.03}_{-0.03}$ |  4.42E+12 |  7.19E+02 |  0.13 |  2.13E+02 |  0.70 |  15 |  Jun 2013 | 19
TYC 5164-567-1 |  0.90 |  0.89 |  4.68 |  0.19${}^{+0.04}_{-0.05}$ |  3.23E+13 |  3.21E+03 |  0.59 |  4.54E+02 |  0.36 |  15 |  Jun 2013 | 19
TYC 6349-0200-1 |  0.85 |  0.96 |  3.41 |  0.07${}^{+0.01}_{-0.02}$ |  4.58E+13 |  3.27E+03 |  0.26 |  8.35E+02 |  0.58 |  15 |  Jun 2013 | 19
TYC 6878-0195-1 |  1.17 |  1.37 |  5.70 |  0.10${}^{+0.04}_{-0.03}$ |  9.02E+13 |  3.17E+03 |  0.19 |  1.40E+03 |  0.85 |  15 |  Jun 2013 | 19
Hot Jupiter Hosts |
$\tau$ Boo |  1.34 |  1.42 |  3 |  >0.732 |  1.52E+11 |  1.61E+00 |  0.48 |  1.54E+00 |  0.86 |  5 |  Jun 2006 | 20, 21, 2, 22
… |  … |  … |  … |  … |  3.11E+11 |  9.42E+00 |  0.59 |  1.80E+00 |  0.65 |  8 |  Jun 2007 | 20, 21, 2, 23
… |  … |  … |  … |  … |  7.83E+10 |  3.10E+00 |  0.13 |  6.16E+00 |  0.91 |  … |  Jan 2008 | 20, 2, 21
… |  … |  … |  … |  … |  7.76E+10 |  3.12E+00 |  0.30 |  5.48E-01 |  0.45 |  … |  Jun 2008 | 20, 2, 21
… |  … |  … |  … |  … |  2.53E+10 |  2.31E+00 |  0.62 |  2.12E-01 |  0.36 |  … |  Jul 2008 | 20, 2, 21
… |  … |  … |  … |  … |  1.71E+11 |  3.75E+00 |  0.56 |  1.18E+00 |  0.63 |  … |  Jun 2009 | 20, 21, 2, 24
HD 73256 |  1.05 |  0.89 |  14 |  0.962 |  3.28E+11 |  4.30E+01 |  0.03 |  1.18E+01 |  0.78 |  4 |  Jan 2008 | 25, 2, 24
HD 102195 |  0.87 |  0.82 |  12.3 |  0.473 |  2.43E+12 |  7.09E+01 |  0.23 |  9.30E+01 |  0.88 |  4 |  Jan 2008 | 26, 27, 2, 24
Table 1: -continued

Star ID |  $M_{\star}$ |  $R_{\star}$ |  $P_{\rm{rot}}$ |  $R_{o}$ |  |$\langle{h\,}\rangle$| |  $\langle{B^{2}_{\rm{pol}}\,}\rangle$ |  Pol |  $\langle{B^{2}_{\rm{tor}}\,}\rangle$ |  Tor |  $l_{\rm{max}}$ |  Obs. | Ref.
---|---|---|---|---|---|---|---|---|---|---|---|---
|  (M⊙) |  (R⊙) |  (d) |  |  (Mx${}^{2}\rm{cm}^{-3}$) |  (G2) |  Axi |  (G2) |  Axi |  |  epoch |
HD 130322 |  0.79 |  0.83 |  26.1 |  0.782 |  1.20E+11 |  5.40E+00 |  0.58 |  1.03E+00 |  0.96 |  4 |  Jan 2008 | 20, 28, 29, 2, 24
HD 179949 |  1.21 |  1.19 |  7.6 |  >1.726 |  5.82E+10 |  5.15E+00 |  0.57 |  1.17E+00 |  0.81 |  6 |  Jun 2007 | 12, 2, 30
HD 189733 |  0.82 |  0.76 |  12.5 |  0.403 |  5.26E+12 |  1.73E+02 |  0.30 |  2.73E+02 |  0.91 |  5 |  Jun 2007 | 31, 2, 32
… |  … |  … |  … |  … |  1.05E+13 |  2.72E+02 |  0.17 |  9.50E+02 |  0.96 |  5 |  Jul 2008 | …
M dwarf stars
GJ 569A |  0.48 |  0.43 |  14.7 |  <0.288 |  1.57E+14 |  1.38E+04 |  0.96 |  6.75E+02 |  1.00 |  5 |  Jan 2008 | 2, 33
DS Leo |  0.58 |  0.52 |  14 |  <0.267 |  1.76E+14 |  2.22E+03 |  0.58 |  1.01E+04 |  0.99 |  5 |  Jan 2007 | 2, 33
… |  … |  … |  … |  … |  8.76E+13 |  2.09E+03 |  0.15 |  8.31E+03 |  0.94 |  … |  Dec 2007 | …
GJ 182 |  0.75 |  0.82 |  4.35 |  0.054 |  5.62E+14 |  1.09E+04 |  0.17 |  2.65E+04 |  0.90 |  8 |  Jan 2007 | 2, 33
GJ 49 |  0.57 |  0.51 |  18.6 |  <0.352 |  2.15E+13 |  4.19E+02 |  0.67 |  4.51E+02 |  1.00 |  5 |  Jul 2007 | 2, 33
GJ 494A |  0.59 |  0.53 |  2.85 |  0.092 |  6.79E+13 |  9.99E+03 |  0.12 |  1.70E+04 |  0.91 |  8 |  2007 | 2, 33
… |  … |  … |  … |  … |  1.95E+14 |  1.09E+04 |  0.27 |  1.44E+04 |  0.89 |  … |  2008 | …
GJ 388 |  0.42 |  0.38 |  2.24 |  0.047 |  9.77E+13 |  4.45E+04 |  0.97 |  3.97E+02 |  0.25 |  8 |  2007 | 2, 34
… |  … |  … |  … |  … |  1.11E+14 |  4.33E+04 |  0.92 |  1.78E+03 |  0.08 |  … |  2008 | …
EQ Peg A |  0.39 |  0.35 |  1.06 |  0.02 |  4.38E+14 |  1.81E+05 |  0.71 |  2.41E+04 |  0.29 |  4 |  Aug 2006 | 2, 34
EQ Peg B |  0.25 |  0.25 |  0.4 |  0.005 |  4.25E+14 |  2.14E+05 |  0.95 |  2.70E+03 |  0.42 |  8 |  Aug 2006 | 2, 34
GJ 873 |  0.32 |  0.3 |  4.37 |  0.068 |  1.33E+15 |  3.45E+05 |  0.28 |  3.05E+04 |  0.61 |  8 |  2006 | 2, 34
… |  … |  … |  … |  … |  2.12E+14 |  2.82E+05 |  0.30 |  4.12E+03 |  0.20 |  … |  2007 | …
GJ 9520 |  0.55 |  0.49 |  3.4 |  0.097 |  1.76E+14 |  1.63E+04 |  0.85 |  4.44E+03 |  0.87 |  8 |  2007 | 2, 33
… |  … |  … |  … |  … |  2.38E+14 |  1.33E+04 |  0.63 |  7.56E+03 |  0.91 |  … |  2008 | …
V374 Peg |  0.28 |  0.28 |  0.45 |  0.006 |  1.29E+14 |  5.30E+05 |  0.82 |  1.43E+04 |  0.01 |  10 |  2005 | 34, 2, 35
… |  … |  … |  … |  … |  8.11E+13 |  3.86E+05 |  0.81 |  1.12E+04 |  0.01 |  … |  2006 | …
GJ 1111 |  0.1 |  0.11 |  0.46 |  0.0059 |  1.84E+13 |  1.29E+04 |  0.79 |  6.54E+02 |  0.68 |  6 |  2007 | 2, 36
… |  … |  … |  … |  … |  9.22E+12 |  5.43E+03 |  0.31 |  1.26E+03 |  0.88 |  … |  2008 | …
… |  … |  … |  … |  … |  1.57E+13 |  4.47E+03 |  0.66 |  1.81E+03 |  0.79 |  … |  2009 | …
GJ 1156 |  0.14 |  0.16 |  0.49 |  0.0081 |  1.39E+12 |  4.38E+03 |  0.02 |  3.74E+02 |  0.19 |  6 |  2007 | 2, 36
… |  … |  … |  … |  … |  2.55E+13 |  1.44E+04 |  0.12 |  1.69E+03 |  0.58 |  … |  2008 | …
… |  … |  … |  … |  … |  3.86E+11 |  1.04E+04 |  0.01 |  4.88E+02 |  0.05 |  … |  2009 | …
GJ 1245B |  0.12 |  0.14 |  0.71 |  0.011 |  2.57E+13 |  3.34E+04 |  0.06 |  5.30E+03 |  0.37 |  4 |  2006 | 2, 36
… |  … |  … |  … |  … |  4.89E+12 |  4.68E+03 |  0.13 |  4.86E+02 |  0.30 |  … |  2008 | …
WX UMa |  0.1 |  0.12 |  0.78 |  0.01 |  2.28E+15 |  2.19E+06 |  0.93 |  4.84E+04 |  0.38 |  4 |  2007 | 2, 36
… |  … |  … |  … |  … |  3.89E+14 |  2.15E+06 |  0.85 |  3.33E+04 |  0.03 |  … |  2008 | …

1: Marsden et al. (2014); 2: Vidotto et al. (2014); 3: Petit et al. (in prep);
4: Saar & Brandenburg (1999); 5: Hempelmann et al. (2016); 6: do Nascimento et
al. (2016); 7: Jeffers et al. (2014); 8: Petit et al. (2008); 9: Morgenthaler
et al. (2011); 10: Fernandes et al. (1998); 11: Jeffers et al. (in prep); 12:
Valenti & Fischer (2005); 13: Mosser et al. (2009); 14: Boro Saikia et al.
(2016); 15: Boro Saikia et al. (2015); 16: Maggio et al. (2000); 17: Innis et
al. (1988); 18: Donati et al. (2003); 19: Folsom et al. (2016); 20: Takeda et
al. (2007); 21: Fares et al. (2009); 22: Catala et al. (2007); 23: Donati et
al. (2008a); 24: Fares et al. (2013); 25: Udry et al. (2003); 26: Melo et al.
(2007); 27: Ge et al. (2006); 28: Udry et al. (2000); 29: Simpson et al.
(2010); 30: Fares et al. (2012); 31: Bouchy et al. (2005); 32: Fares et al.
(2010); 33: Donati et al. (2008b); 34: Morin et al. (2008b); 35: Morin et al.
(2008a); 36: Morin et al. (2010)

## Acknowledgements

The authors would like to thank the referee for thoughtful and constructive
comments. MJ and KL acknowledge support from Science and Technology Facilities
Council (STFC) consolidated grant number ST/R000824/1. VS acknowledges funding
from the European Research Council (ERC) under the European Unions Horizon
2020 research and innovation programme (grant agreement no. 682393,
AWESoMeStars). JFD acknowledges funding from the European Research Council
(ERC) under the H2020 research and innovation programme (grant agreement
740651 NewWorlds). RF acknowledges funding from United Arab Emirates
University (UAEU) startup grant number G00003269.

## Data Availability

The data used to make Figures 3, 5 and 7 are in Table 1 and can be accessed at
https://doi.org/10.17630/4dea9f0c-9626-4655-ba61-8af5c05fbe38. Archival data
underpinning the plots is available at polarbase
(http://polarbase.irap.omp.eu).

## References

  * Berger (1985) Berger M. A., 1985, The Astrophysical Journal Supplement Series, 59, 433
  * Berger & Hornig (2018) Berger M. A., Hornig G., 2018, Journal of Physics A Mathematical General, 51, 495501
  * Boro Saikia et al. (2015) Boro Saikia S., Jeffers S. V., Petit P., Marsden S., Morin J., Folsom C. P., 2015, A&A, 573, A17
  * Boro Saikia et al. (2016) Boro Saikia S., et al., 2016, A&A, 594, A29
  * Bouchy et al. (2005) Bouchy F., et al., 2005, A&A, 444, L15
  * Brandenburg & Subramanian (2005) Brandenburg A., Subramanian K., 2005, Physics Reports, 417, 1
  * Catala et al. (2007) Catala C., Donati J.-F., Shkolnik E., Bohlender D., Alecian E., 2007, MNRAS, 374, L42
  * Chandrasekhar (1961) Chandrasekhar S., 1961, Hydrodynamic and Hydromagnetic Stability. Clarendon Press, https://books.google.co.uk/books?id=XvZ0tgEACAAJ
  * Chatterjee et al. (2011) Chatterjee P., Guerrero G., Brandenburg A., 2011, A&A, 525, A5
  * Démoulin (2007) Démoulin P., 2007, Advances in Space Research, 39, 1674
  * Démoulin & Pariat (2009) Démoulin P., Pariat E., 2009, Advances in Space Research, 43, 1013
  * Donati & Landstreet (2009) Donati J. F., Landstreet J. D., 2009, ARA&A, 47, 333
  * Donati et al. (2003) Donati J.-F., et al., 2003, MNRAS, 345, 1145
  * Donati et al. (2006) Donati J.-F., et al., 2006, MNRAS, 370, 629
  * Donati et al. (2008a) Donati J.-F., et al., 2008a, MNRAS, 385, 1179
  * Donati et al. (2008b) Donati J.-F., et al., 2008b, MNRAS, 390, 545
  * Fares et al. (2009) Fares R., et al., 2009, MNRAS, 398, 1383
  * Fares et al. (2010) Fares R., et al., 2010, MNRAS, 406, 409
  * Fares et al. (2012) Fares R., et al., 2012, MNRAS, 423, 1006
  * Fares et al. (2013) Fares R., Moutou C., Donati J.-F., Catala C., Shkolnik E. L., Jardine M. M., Cameron A. C., Deleuil M., 2013, MNRAS, 435, 1451
  * Fernandes et al. (1998) Fernandes J., Lebreton Y., Baglin A., Morel P., 1998, A&A, 338, 455
  * Folsom et al. (2016) Folsom C. P., et al., 2016, MNRAS, 457, 580
  * Ge et al. (2006) Ge J., et al., 2006, ApJ, 648, 683
  * Hempelmann et al. (2016) Hempelmann A., Mittag M., Gonzalez-Perez J. N., Schmitt J. H. M. M., Schröder K. P., Rauw G., 2016, A&A, 586, A14
  * Innis et al. (1988) Innis J. L., Thompson K., Coates D. W., Evans T. L., 1988, MNRAS, 235, 1411
  * Jeffers et al. (2014) Jeffers S. V., Petit P., Marsden S. C., Morin J., Donati J.-F., Folsom C. P., 2014, A&A, 569, A79
  * Lehmann et al. (2019) Lehmann L. T., Hussain G. A. J., Jardine M. M., Mackay D. H., Vidotto A. A., 2019, MNRAS, 483, 5246
  * Lund et al. (2020) Lund K., et al., 2020, MNRAS, 493, 1003
  * Maggio et al. (2000) Maggio A., Pallavicini R., Reale F., Tagliaferri G., 2000, A&A, 356, 627
  * Marsden et al. (2014) Marsden S. C., et al., 2014, MNRAS, 444, 3517
  * Melo et al. (2007) Melo C., et al., 2007, A&A, 467, 721
  * Morgenthaler et al. (2011) Morgenthaler A., Petit P., Morin J., Aurière M., Dintrans B., Konstantinova-Antova R., Marsden S., 2011, Astronomische Nachrichten, 332, 866
  * Morin et al. (2008a) Morin J., et al., 2008a, MNRAS, 384, 77
  * Morin et al. (2008b) Morin J., et al., 2008b, MNRAS, 390, 567
  * Morin et al. (2010) Morin J., Donati J. F., Petit P., Delfosse X., Forveille T., Jardine M. M., 2010, MNRAS, 407, 2269
  * Morin et al. (2011) Morin J., Dormy E., Schrinner M., Donati J. F., 2011, MNRAS, 418, L133
  * Mosser et al. (2009) Mosser B., et al., 2009, A&A, 506, 33
  * Petit et al. (2008) Petit P., et al., 2008, MNRAS, 388, 80
  * Pipin et al. (2019) Pipin V. V., Pevtsov A. A., Liu Y., Kosovichev A. G., 2019, The Astrophysical Journal, 877, L36
  * Saar & Brandenburg (1999) Saar S. H., Brandenburg A., 1999, ApJ, 524, 295
  * See et al. (2015) See V., et al., 2015, MNRAS, 453, 4301
  * Semel (1989) Semel M., 1989, A&A, 225, 456
  * Simpson et al. (2010) Simpson E. K., Baliunas S. L., Henry G. W., Watson C. A., 2010, MNRAS, 408, 1666
  * Takeda et al. (2007) Takeda G., Ford E. B., Sills A., Rasio F. A., Fischer D. A., Valenti J. A., 2007, ApJS, 168, 297
  * Taylor (1974) Taylor J. B., 1974, Physical Review Letters, 33, 1139
  * Udry et al. (2000) Udry S., et al., 2000, A&A, 356, 590
  * Udry et al. (2003) Udry S., et al., 2003, A&A, 407, 679
  * Valenti & Fischer (2005) Valenti J. A., Fischer D. A., 2005, ApJS, 159, 141
  * Vidotto (2016) Vidotto A. A., 2016, MNRAS, 459, 1533
  * Vidotto et al. (2014) Vidotto A. A., et al., 2014, MNRAS, 441, 2361
  * Woltjer (1958) Woltjer L., 1958, Proceedings of the National Academy of Science, 44, 489
  * do Nascimento et al. (2016) do Nascimento J. D. J., et al., 2016, ApJ, 820, L15


# Visibility-Aware RRT* for Safety-Critical Navigation of Perception-Limited
Robots in Unknown Environments

Taekyung Kim and Dimitra Panagou Taekyung Kim is with the Department of
Robotics, University of Michigan, Ann Arbor, MI, 48109, USA
<EMAIL_ADDRESS>Panagou is with the Department of Robotics and
Department of Aerospace Engineering, University of Michigan, Ann Arbor, MI,
48109, USA<EMAIL_ADDRESS>work has been partially supported by Amazon.

###### Abstract

Safe autonomous navigation in unknown environments remains a critical
challenge for robots with limited sensing capabilities. While safety-critical
control techniques, such as Control Barrier Functions (CBFs), have been
proposed to ensure safety, their effectiveness relies on the assumption that
the robot has complete knowledge of its surroundings. In reality, robots often
operate with restricted field-of-view and finite sensing range, which can lead
to collisions with unknown obstacles if the planning algorithm is agnostic to
these limitations. To address this issue, we introduce the visibility-aware
RRT* algorithm that combines sampling-based planning with CBFs to generate
safe and efficient global reference paths in partially unknown environments.
The algorithm incorporates a collision avoidance CBF and a novel visibility
CBF, which guarantees that the robot remains within locally collision-free
regions, enabling timely detection and avoidance of unknown obstacles. We
conduct extensive experiments interfacing the path planners with two different
safety-critical controllers, wherein our method outperforms all other compared
baselines across both safety and efficiency aspects. 111Our project page can
be found at: https://www.taekyung.me/visibility-rrt

##  I INTRODUCTION

Autonomous navigation in unknown environments is a fundamental challenge in
robotics, with applications ranging from exploration and search-and-rescue to
autonomous driving. A general navigation stack comprises four key components:
localization, mapping, planning, and control. Localization and mapping enable
robots to build a map of their surroundings using onboard sensors, such as
Visual Simultaneous Localization and Mapping (VSLAM) systems. However, these
sensors often have limitations, including restricted field-of-view (FOV) and
finite sensing range. For instance, popular off-the-shelf RGB-D cameras’
horizontal FOV are often less than 70∘ [1].

To ensure safe navigation, there has been growing interest in safety-critical
controllers, such as those based on Control Barrier Functions (CBFs) [2],
which enforce safety objectives on controlled systems. These controllers
theoretically guarantee safety by keeping the robot within a safe set,
assuming that nearby obstacles are fully known a priori. However, in practice,
the effectiveness of these controllers relies heavily on the robot’s ability
to perceive potential obstacles in the local environment. If the global
planner does not account for the robot’s limited sensing capabilities under
the presence of unknown obstacles, the theoretical safety guarantees may not
hold, leading to potentially hazardous situations.

Several approaches have been proposed to improve the safety of navigation
under sensing limitations. Some methods use nonlinear optimization to generate
[3] or refine [4] perception-aware paths, but these approaches can be
computationally intensive and the solution convergence is not guaranteed.
Other approaches first plan a path and then optimize the heading angle to
enhance safety based on sensing objectives [5, 6, 7], although they do not
always guarantee safety. Another method computes local safe regions during the
exploration of candidate paths to account for local sensing information and
ensure safety [8], but it requires storing all combinations of such regions,
making the planning module overly computationally expensive.

In this paper, we propose a visibility-aware RRT* algorithm that explicitly
reasons about the robot’s limited sensing capabilities and combines the
benefits of sampling-based planning with CBFs to enable safe and efficient
navigation in partially unknown environments. Our approach incorporates two
types of CBFs: (i) a collision avoidance CBF to ensure that the planned path
is collision-free with respect to known obstacles, and (ii) a novel visibility
CBF, which we introduce in this work, that guarantees the local tracking
controller following the resulting path will always keep the robot within
locally collision-free regions. This ensures that the robot can detect unknown
obstacles in a timely manner, enabling the robot to avoid them by virtue of
the local controller. We demonstrate the effectiveness of our method through
extensive simulations and comparisons with other planning methods. The results
indicate that our method can be seamlessly incorporated with safety-critical
controllers, outperforming existing baselines in both safety and efficiency
metrics across multiple scenarios.

##  II PRELIMINARIES

####  II-1 System Description

We consider a robot that is modeled as a continuous-time, control-affine
system:

$\dot{{\bm{x}}}=F({\bm{x}},{\bm{u}})=f({\bm{x}})+g({\bm{x}}){\bm{u}},$ (1)

where ${\bm{x}}\in\mathcal{X}\subset\mathbb{R}^{n}$ is the state, and
${\bm{u}}\in\mathcal{U}\subset\mathbb{R}^{m}$ is the control input, with
$\mathcal{U}$ being a set of admissible controls for System (1). The functions
$f:\mathcal{X}\to\mathbb{R}^{n}$ and $g:\mathcal{U}\to\mathbb{R}^{n\times m}$
are assumed to be locally Lipschitz continuous.

We assume the robot has limited sensing capabilities, modeled via a field-of-
view angle $\theta_{\text{fov}}<360^{\circ}$ and a maximum sensing range
$l_{\text{range}}<\infty$ of the onboard sensor.

####  II-2 Linear Quadratic Regulator

Linear Quadratic Regulator (LQR) computes optimal control policies for linear
systems, and can be applied to nonlinear systems by linearizing the system
around an operating point $({\bm{x}}_{o},{\bm{u}}_{o})$, i.e., taking the
first-order Taylor expansion of a nonlinear dynamics $F$ around this point:

$\dot{{\bm{x}}}\approx
F({\bm{x}}_{o},{\bm{u}}_{o})+A\bar{{\bm{x}}}+B\bar{{\bm{u}}},$ (2)

where $\bar{{\bm{x}}}={\bm{x}}-{\bm{x}}_{o}$,
$\bar{{\bm{u}}}={\bm{u}}-{\bm{u}}_{o}$, $A=\frac{\partial
F}{\partial{\bm{x}}}|_{({\bm{x}}_{o},{\bm{u}}_{o})}$, $B=\frac{\partial
F}{\partial{\bm{u}}}|_{({\bm{x}}_{o},{\bm{u}}_{o})}$. Given an infinite time
horizon cost function:

$J=\int_{0}^{\infty}(\,\bar{{\bm{x}}}^{\top}Q\bar{{\bm{x}}}+\bar{{\bm{u}}}^{\top}R\bar{{\bm{u}}}\,)\,\mathrm{d}t,$
(3)

where $Q=Q^{\top}\succeq 0$ and $R=R^{\top}\succ 0$ are weight matrices, one
can compute the locally optimal control policy
$\mathcal{K}_{\text{lqr}}^{*}({\bm{x}})=-K_{\text{lqr}}{\bm{x}}$ by solving
the algebraic Riccati equation:

$A^{\top}P+PA-PBR^{-1}B^{\top}P+Q=0.$ (4)

Solving (4) for $P$, the LQR gain matrix is:

$K_{\text{lqr}}=R^{-1}B^{\top}P.$ (5)

####  II-3 Control Barrier Functions

######  Definition 1 (Control Barrier Function (CBF)).

Let $\mathcal{C}=\\{{\bm{x}}\in\mathcal{X}:h({\bm{x}})\geq 0\\}$, where
$h:\mathcal{X}\rightarrow\mathbb{R}$ is a continuously differentiable
function. The function $h$ is a CBF on $\mathcal{C}$ for System (1) if there
exists a class $\mathcal{K}$ function $\alpha$ such that

$\sup_{{\bm{u}}\in\mathcal{U}}\left[L_{f}h({\bm{x}})+L_{g}h({\bm{x}}){\bm{u}}\right]\geq-\alpha(h({\bm{x}})),$
(6)

for all ${\bm{x}}\in\mathcal{X}$ and $L_{f}h({\bm{x}})$ and $L_{g}h({\bm{x}})$
denote the Lie derivatives of $h$ along $f$ and $g$, respectively [9].

######  Theorem 1.

Given a CBF $h$ with the associated set $\mathcal{C}$, any Lipschitz
continuous controller ${\bm{u}}\in\mathcal{K}_{\textup{cbf}}({\bm{x}})$, with
$\mathcal{K}_{\textup{cbf}}({\bm{x}})=\\{{\bm{u}}\in\mathcal{U}:L_{f}h({\bm{x}})+L_{g}h({\bm{x}}){\bm{u}}+\alpha(h({\bm{x}}))\geq
0\\}$, renders $\mathcal{C}$ forward invariant for System (1). See [9].

####  II-4 High Order CBF

The concept of CBFs has been extended to High Order CBFs (HOCBFs) [10], which
provide a more general formulation that can be used for constraints of high
relative degree. In this work, we focus on time-invariant HOCBFs, which
generalize the concept of exponential CBFs proposed in [11].

For a $r_{b}^{\textup{th}}$ order differentiable function
$h:\mathcal{X}\rightarrow\mathbb{R}$, we define a series of functions
$\psi_{i}:\mathcal{X}\rightarrow\mathbb{R}$, $i=0,\ldots,r_{b}$, as:

$\displaystyle\psi_{0}({\bm{x}})$ $\displaystyle:=h({\bm{x}}),$ (7a)
$\displaystyle\psi_{1}({\bm{x}})$
$\displaystyle:=\dot{\psi}_{0}({\bm{x}})+\alpha_{1}(\psi_{0}({\bm{x}})),$ (7b)
$\displaystyle\vdots$ $\displaystyle\psi_{r_{b}}({\bm{x}})$
$\displaystyle:=\dot{\psi}_{r_{b}-1}({\bm{x}})+\alpha_{r_{b}}(\psi_{r_{b}-1}({\bm{x}})),$
(7c)

where $\alpha_{i}:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$,
$i=1,\ldots,r_{b}$, denote class $\mathcal{K}$ functions, and
$\dot{\psi}_{i}({\bm{x}})=L_{f}\psi_{i}({\bm{x}})+L_{g}\psi_{i}({\bm{x}}){\bm{u}}$.
We further define a series of sets $\mathcal{C}_{i}$ associated with these
functions as:

$\displaystyle\mathcal{C}_{1}$
$\displaystyle:=\\{{\bm{x}}\in\mathcal{X}:\psi_{0}({\bm{x}})\geq 0\\},$ (8a)
$\displaystyle\mathcal{C}_{2}$
$\displaystyle:=\\{{\bm{x}}\in\mathcal{X}:\psi_{1}({\bm{x}})\geq 0\\},$ (8b)
$\displaystyle\vdots$ $\displaystyle\mathcal{C}_{r_{b}}$
$\displaystyle:=\\{{\bm{x}}\in\mathcal{X}:\psi_{r_{b}-1}({\bm{x}})\geq 0\\}.$
(8c)

######  Definition 2 (HOCBF).

The function $h:\mathcal{X}\rightarrow\mathbb{R}$ is a HOCBF of relative
degree $r_{b}$ for System (1) if there exist differentiable class
$\mathcal{K}$ functions $\alpha_{1},\alpha_{2},\ldots,\alpha_{r_{b}}$ such
that

$\psi_{r_{b}}({\bm{x}})\geq
0,\quad\forall{\bm{x}}\in\mathcal{C}_{1}\cap\mathcal{C}_{2}\cap\ldots\cap\mathcal{C}_{r_{b}}.$
(9)

######  Theorem 2.

Given a HOCBF $h$ with the associated sets
$\mathcal{C}_{1},\mathcal{C}_{2},\ldots,\mathcal{C}_{r_{b}}$, if
${\bm{x}}(t_{0})\in\mathcal{C}_{1}\cap\mathcal{C}_{2}\cap\ldots\cap\mathcal{C}_{r_{b}}$,
then any Lipschitz continuous controller
${\bm{u}}(t)\in\mathcal{K}_{\textup{hocbf}}({\bm{x}})$, with
$\mathcal{K}_{\textup{hocbf}}({\bm{x}})=\\{{\bm{u}}\in\mathcal{U}:\psi_{r_{b}}({\bm{x}})\geq
0\\}$, renders the set
$\mathcal{C}_{1}\cap\mathcal{C}_{2}\cap\ldots\cap\mathcal{C}_{r_{b}}$ forward
invariant for system (1) [10].

##  III PROBLEM FORMULATION

Consider a robot modeled as the nonlinear system (1) navigating in an unknown
environment $\mathcal{W}$. The environment $\mathcal{W}$ consists of two types
of obstacles: a set $\mathcal{O}_{0}\subset\mathcal{W}$ of known obstacles at
the initial time $t_{0}$, and a set $\mathcal{H}\subset\mathcal{W}$ of unknown
(or “hidden”) obstacles that can be sensed by the robot on-the-fly as it
navigates through the environment. The true obstacle-free space $\mathcal{S}$
is then defined as
$\mathcal{S}=\mathcal{W}\setminus(\mathcal{O}_{0}\cup\mathcal{H})$. A local
tracking control policy
$\pi:[t_{0},\infty)\times\mathcal{X}\times\mathcal{X}\to\mathcal{U}$ takes the
current time $t$, the current state ${\bm{x}}_{t}:=\bm{x}(t)$, and a reference
state ${\bm{x}}_{t}^{\textup{ref}}\in\mathcal{X}$ as inputs, and generates as
outputs the control inputs that track a nominal path from the current to
reference state. Given the initial condition ${\bm{x}}(t_{0})={\bm{x}}_{0}$,
$t\geq t_{0}$, the robot dynamics obey:

$\dot{{\bm{x}}}(t)=f({\bm{x}}(t))+g({\bm{x}}(t))\pi(t,{\bm{x}}_{t},{\bm{x}}_{t}^{\textup{ref}}),\quad{\bm{x}}(t_{0})={\bm{x}}_{0}.$
(10)

The navigation problem is divided into two hierarchical steps: 1) finding a
global reference path $p^{\textup{ref}}$, and 2) tracking $p^{\textup{ref}}$
using a local controller.

Global Path Planning: Given a goal state
${\bm{x}}_{\textup{goal}}\in\mathcal{X}$, the objective of the global path
planner is to find a global reference path
$p^{\textup{ref}}:\mathcal{X}\to\mathcal{P}$, represented by a set of
waypoints
$\mathcal{P}=\\{{\bm{x}}_{i}^{\textup{ref}}\in\mathcal{X}\\}_{i=0}^{N}$, where
${\bm{x}}_{0}^{\textup{ref}}={\bm{x}}_{0}$ and
${\bm{x}}_{N}^{\textup{ref}}={\bm{x}}_{\textup{goal}}$ are the initial and
final waypoints in $\mathcal{P}$, respectively. For any given state
${\bm{x}}_{t}$, the function yields the next waypoint that the robot should
navigate towards as:

${\bm{x}}_{i}^{\textup{ref}}=p^{\textup{ref}}({\bm{x}}_{t}),$ (11)

where ${\bm{x}}_{i}^{\textup{ref}}\in\mathcal{P}\subset\mathcal{X}$ indicates
the immediate waypoint for navigation. The reference path should satisfy:

$p^{\textup{ref}}({\bm{x}}_{t})\in\mathcal{S},\quad\forall t\geq
t_{0},\quad\forall{\bm{x}}_{t}\in\mathcal{X}.$ (12)

However, the true collision-free set $\mathcal{S}$ is often unknown in
advance.

######  Definition 3 (Collision-Free Path).

We define the known collision-free set at $t_{0}$ as
$\mathcal{S}_{0}=\mathcal{W}\setminus\mathcal{O}_{0}$. A reference path
$p^{\textup{ref}}$ is considered collision-free if it satisfies the condition:

$p^{\textup{ref}}({\bm{x}}_{t})\in\mathcal{S}_{0},\quad\forall t\geq
t_{0},\quad\forall{\bm{x}}_{t}\in\mathcal{X}.$ (13)

Let ${\bm{r}}_{t}=[x_{t},y_{t}]^{\top}\in\mathbb{R}^{2}$ and
$\theta_{t}\in\mathbb{R}$ denote the position and orientation of the robot
with respect to a global frame at time $t$, respectively. The sensing
footprint $\mathcal{F}_{t}$, which represents the currently sensed area at
time $t$, is defined as:

$\displaystyle\mathcal{F}_{t}=\bigg{\\{}{\bm{p}}\in\mathcal{S}\mid$
$\displaystyle\,\|{\bm{p}}-{\bm{r}}_{t}\|\leq l_{\textup{range}}$
$\displaystyle\cap\,|\,\angle({\bm{p}}-{\bm{r}}_{t},\theta_{t})\,|\leq\frac{\theta_{\textup{fov}}}{2}\bigg{\\}},$
(14)

where $\|{\bm{p}}-{\bm{r}}_{t}\|$ represents the Euclidean distance between a
point ${\bm{p}}\in\mathcal{W}$ and the robot’s position ${\bm{r}}_{t}$,
$\angle({\bm{p}}-{\bm{r}}_{t},\theta_{t})$ denotes the angle between the
vector ${\bm{p}}-{\bm{r}}_{t}$ and the robot’s orientation $\theta_{t}$. Given
the accumulated sensory information collected from time $t_{0}$ to $t$, we
define the known collision-free space $\mathcal{B}_{t}$ at time $t$ as:

$\mathcal{B}_{t}:=\bigcup_{\tau=t_{0}}^{t}\mathcal{F}_{\tau}\subseteq\mathcal{S}.$
(15)

Consequently, the part of the environment that has not been sensed yet is
$\mathcal{S}\setminus\mathcal{B}_{t}$, and is called the unmeasured space at
time $t$.

It is important to note that the collision-free path criterion (13) treats all
unmeasured space $\mathcal{S}\setminus\mathcal{B}_{t_{0}}$ at time $t_{0}$ as
free space. However, while following the path $p^{\textup{ref}}$ using the
local controller $\pi$ after the planning cycle, unknown obstacles
$\mathcal{H}$ might exist within such unmeasured spaces.

Given this context, we introduce a stronger criterion that underlines the
path’s visibility and safety by requiring the robot to stay within the known
local collision-free space $\mathcal{B}_{t}$, thereby mitigating the risk
posed by unknown obstacles [8].

######  Definition 4 (Visibility-Aware Path).

A reference path $p^{\textup{ref}}$ is called visibility-aware w.r.t. the
dynamics (1) if the updated state ${\bm{x}}(t+\Delta t)$, for a small time
increment $\Delta t>0$, under the local tracking controller
$\pi(t,{\bm{x}}_{t},p^{\textup{ref}}({\bm{x}}_{t}))$ satisfies:

${\bm{x}}(t+\Delta t)\in\mathcal{B}_{t},\quad\forall t\geq
t_{0},\quad\forall{\bm{x}}_{t}\in\mathcal{X}.$ (16)

In the following sections, we present the construction of a global path
planner that imposes (16) (see Sec. IV), and demonstrate how (16) ensures
hidden obstacles are avoidable, in the sense that the local tracking
controller is able to generate collision-free trajectories (see Sec. V).

Local Tracking Control: A local control policy $\pi$ is used to stabilize the
systems around the subsequent waypoint ${\bm{x}}_{i}^{\textup{ref}}$ (11) and
to ensure that the robot remains within the true collision-free set
$\mathcal{S}$. The local controller is assumed to be capable of tracking the
path between subsequent waypoints with a tracking error bounded by a maximum
displacement $\epsilon$. As $p^{\textup{ref}}$ continuously updates the
waypoint, it drives the robot towards the goal state
${\bm{x}}_{\textup{goal}}$. There exists a plethora of methodologies to design
such controllers [9, 12]. Since the primary focus of this paper is on
designing a visibility-aware path planner, we will discuss two case studies of
such controllers in Sec. V and Sec. VI to highlight the advantages of our
proposed planner when utilized alongside various local controllers.

##  IV GLOBAL PATH PLANNER

In this section, we introduce the Visibility-Aware RRT* algorithm, which aims
to solve the motion planning problem under limited sensing capabilities while
satisfying the visibility condition given in (16). The proposed algorithm
builds upon the LQR-RRT* framework [13] by incorporating two CBFs that guide
the search process to ensure satisfaction of both the collision-free (13) and
visibility-aware (16) properties.

###  IV-A Visibility-Aware RRT*

Initialization: $\mathcal{V}\leftarrow$ $\\{{\bm{x}}_{0}\\}$;
$\mathcal{E}\leftarrow\emptyset$;

for  _$i=1,\ldots,\texttt{{maxIter}}$_ do

${\bm{x}}_{\textup{rand}}\leftarrow\texttt{Sample}()$;

${\bm{x}}_{\textup{nearest}}\leftarrow\texttt{Nearest}(\mathcal{V},{\bm{x}}_{\textup{rand}})$;

${\bm{x}}_{\textup{rand}}\leftarrow\texttt{SetAngle}({\bm{x}}_{\textup{rand}},{\bm{x}}_{\textup{nearest}})$;

${\bm{x}}_{\textup{new}},\\_\leftarrow\texttt{LQR-CBF-
Steer}({\bm{x}}_{\textup{nearest}},{\bm{x}}_{\textup{rand}})$;

$\mathcal{X}_{\textup{near}}\leftarrow\texttt{NearbyNode}(\mathcal{V},{\bm{x}}_{\textup{new}})$;

${\bm{x}}_{\textup{min}}\leftarrow\texttt{ChooseParent}(\mathcal{X}_{\textup{near}},{\bm{x}}_{\textup{new}})$;

if  _${\bm{x}}_{\textup{min}}\neq\texttt{{None}}$_ then

$\mathcal{V}\leftarrow\mathcal{V}\cup\\{{\bm{x}}_{\textup{new}}\\}$;

$\mathcal{E}\leftarrow\mathcal{E}\cup\\{({\bm{x}}_{\textup{min}},{\bm{x}}_{\textup{new}})\\}$;

$\texttt{Rewire}(\mathcal{X}_{\textup{near}},{\bm{x}}_{\textup{new}})$;

$\mathcal{P}\leftarrow\texttt{ExtractPath}(\mathcal{V},{\bm{x}}_{\textup{0}},{\bm{x}}_{\textup{goal}})$;

return $\mathcal{P}$;

Algorithm 1 Visibility-Aware RRT*

The LQR-RRT* [13] is a sampling-based motion planning algorithm that
incorporates the LQR controller as a steering function for nonlinear systems
with non-holonomic constraints. It guarantees probabilistic completeness and
inherits the asymptotic optimality property from the RRT* algorithm through
its rewiring process [14]. Recently, Yang et al. [15] proposed an extension to
the LQR-RRT* that incorporates CBF constraints during the planning process.
Instead of formulating and solving Quadratic Programs (QPs) iteratively, as
done in existing works [16], it checks whether the CBF constraints are
satisfied and uses them as the termination condition for the steering process.

In our work, we modify the LQR-RRT* baseline to include the robot’s
orientation angle $\theta$ at each node, in addition to the position
coordinates $(x,y)$. This modification allows it to consider the robot’s
rotational motion and field of view during the planning process, essential for
accommodating the robot’s sensing limitations.

The overall algorithm, as shown in Alg. 1, iteratively expands a tree
$(\mathcal{V},\mathcal{E})$ in the configuration space $\mathcal{W}$ until a
maximum number of iterations (denoted as maxIter) is reached. The components
of the algorithm are as follows:

  * •

Sample: Randomly samples a state
${\bm{x}}_{\text{rand}}=[x_{\text{rand}},y_{\text{rand}},\theta_{\text{rand}}]^{\top}$
within the configuration space to expand the tree towards.

  * •

Nearest: Finds the nearest node ${\bm{x}}_{\text{nearest}}$ in the current
vertices in the tree $\mathcal{V}$ to the sampled state
${\bm{x}}_{\text{rand}}$.

  * •

SetAngle: Assigns the orientation angle $\theta_{\text{rand}}$ to the sampled
state ${\bm{x}}_{\text{rand}}$ based on the direction from the nearest node
${\bm{x}}_{\text{nearest}}$ to ${\bm{x}}_{\text{rand}}$:
$\theta_{\text{rand}}=\arctan\left(\frac{y_{\text{rand}}-y_{\text{nearest}}}{x_{\text{rand}}-x_{\text{nearest}}}\right)$.

  * •

LQR-CBF-Steer: Generates a new node ${\bm{x}}_{\text{new}}$ by steering from
${\bm{x}}_{\text{nearest}}$ towards ${\bm{x}}_{\text{rand}}$. The details are
described below.

  * •

NearbyNode: Identifies a set of nearby nodes $\mathcal{X}_{\text{near}}$
within a pre-defined euclidean distance from ${\bm{x}}_{\text{new}}$.

  * •

ChooseParent: Selects the best parent node ${\bm{x}}_{\text{min}}$ for the new
node ${\bm{x}}_{\text{new}}$ from the set of nearby nodes
$\mathcal{X}_{\text{near}}$ based on the cost and the feasibility of the
connecting trajectory generated by the steering function. Given a trajectory
$\bm{\sigma}=({\bm{x}}(t),{\bm{u}}(t))$, where $t\in[0,T]$ and $T$ is the
duration of the trajectory, the evaluated cost is defined as:

$\texttt{cost}(\bm{\sigma})=\int_{0}^{T}\left({\bm{x}}(t)^{\top}Q{\bm{x}}(t)+{\bm{u}}(t)^{\top}R{\bm{u}}(t)\right)\,\mathrm{d}t,$
(17)

where the $Q$ and $R$ are the same weight matrices as in (3) (see Alg. 2).

  * •

Rewire: Performs the rewiring process by considering the new node
${\bm{x}}_{\text{new}}$ as a potential parent for the nearby nodes in
$\mathcal{X}_{\text{near}}$. If a trajectory from ${\bm{x}}_{\text{new}}$ to
any nearby node ${\bm{x}}_{\text{near}}\in\mathcal{X}_{\text{near}}$ is found
with a smaller cost (17), the tree is optimized by rewiring
${\bm{x}}_{\text{near}}$ to ${\bm{x}}_{\text{new}}$ (see Alg. 3).

  * •

ExtractPath: Searches for a path from ${\bm{x}}_{\text{goal}}$ to
${\bm{x}}_{\text{start}}$ within the tree $(\mathcal{V},\mathcal{E})$. If a
feasible path is found, it returns the solution path as a set of waypoints
$\mathcal{P}$.

minCost $\leftarrow\infty$;
${\bm{x}}_{\textup{min}}\leftarrow\texttt{{None}}$;

for  _${\bm{x}}_{\textup{near}}\in\mathcal{X}_{\textup{near}}$_ do

$\\_,\bm{\sigma}\leftarrow\texttt{LQR-CBF-
Steer}({\bm{x}}_{\textup{near}},{\bm{x}}_{\textup{new}})$;

if  _$\bm{\sigma}\neq\texttt{{None}}$_ then

if  _${\bm{x}}_{\textup{near}}.\texttt{{cost}}+\texttt{{Cost}}(\bm{\sigma})
<\texttt{{minCost}}$_ then

minCost
$\leftarrow{\bm{x}}_{\textup{near}}.\texttt{cost}+\texttt{Cost}(\bm{\sigma})$;

${\bm{x}}_{\textup{min}}\leftarrow{\bm{x}}_{\textup{near}}$;

${\bm{x}}_{\text{new}}.\texttt{parent}\leftarrow{\bm{x}}_{\text{min}}$;

${\bm{x}}_{\text{new}}.\texttt{cost}\leftarrow\texttt{minCost}$;

return ${\bm{x}}_{\textup{min}}$;

Algorithm 2 ChooseParent($\mathcal{X}_{\textup{near}}$,
${\bm{x}}_{\textup{new}}$)

for  _${\bm{x}}_{\textup{near}}\in\mathcal{X}_{\textup{near}}$_ do

$\\_,\bm{\sigma}\leftarrow\texttt{LQR-CBF-
Steer}({\bm{x}}_{\textup{new}},{\bm{x}}_{\textup{near}})$;

if  _$\bm{\sigma}\neq\texttt{{None}}$_ then

if  _${\bm{x}}_{\textup{new}}.\texttt{{cost}}+\texttt{{Cost}}(\bm{\sigma})
<{\bm{x}}_{\textup{near}}.\texttt{{cost}}$_ then

${\bm{x}}_{\textup{near}}.\texttt{parent}\leftarrow{\bm{x}}_{\textup{new}}$;

${\bm{x}}_{\textup{new}}.\texttt{cost}\leftarrow\texttt{{Cost}}(\bm{\sigma})$;

Algorithm 3 Rewire($\mathcal{X}_{\textup{near}}$, ${\bm{x}}_{\textup{new}}$)

$\mathcal{X}^{\prime}\leftarrow[{\bm{x}}_{\textup{start}}]$;
$\mathcal{U}^{\prime}\leftarrow[]$;
${\bm{x}}^{\prime}_{0}\leftarrow{\bm{x}}_{\textup{start}}$;

$K_{\textup{lqr}}\leftarrow\texttt{LQRsolver}({\bm{x}}_{\textup{next}})$

for  _$t=1,\ldots,T$_ do

${\bm{x}}^{\prime}_{t},{\bm{u}}^{\prime}_{t}\leftarrow\texttt{Integrator}({\bm{x}}^{\prime}_{t-1},K_{\textup{lqr}})$;

if  _
$\forall\,\texttt{{CBFConstraint}}({\bm{x}}^{\prime}_{t},{\bm{u}}^{\prime}_{t})$
are satisfied_ then

$\mathcal{X}^{\prime}.\texttt{append}({\bm{x}}^{\prime}_{t})$;
$\mathcal{U}^{\prime}.\texttt{append}({\bm{u}}^{\prime}_{t})$;

else

break;

${\bm{x}}_{\textup{new}}\leftarrow{\bm{x}}^{\prime}_{t}$;

return
${\bm{x}}_{\textup{new}},\bm{\sigma}=(\mathcal{X}^{\prime},\mathcal{U}^{\prime})$;

Algorithm 4 LQR-CBF-
Steer(${\bm{x}}_{\textup{start}},{\bm{x}}_{\textup{next}}$)

The LQR-CBF-Steer is the key function that allows the algorithm to use LQR (5)
and CBF (6) as an extension heuristic between two nodes
${\bm{x}}_{\text{start}}$ and ${\bm{x}}_{\text{next}}$ (see Alg. 4). It first
linearizes the dynamics around ${\bm{x}}_{\text{next}}$ (2) and computes the
optimal gain $K_{\text{lqr}}$ using the Riccatti equation (5). The function
then employs the LQR controller $\mathcal{K}_{\text{lqr}}^{*}$ to generate a
sequence of intermediate states
$\mathcal{X}^{\prime}=\\{{\bm{x}}_{\text{start}},{\bm{x}}^{\prime}_{1},\ldots,{\bm{x}}^{\prime}_{T}\\}$
and control inputs
$\mathcal{U}^{\prime}=\\{{\bm{u}}^{\prime}_{1},\ldots,{\bm{u}}^{\prime}_{T}\\}$,
steering the state from ${\bm{x}}_{\text{start}}$ towards
${\bm{x}}_{\text{next}}$ based on the linearized dynamics. Throughout this
process, it verifies that all CBF constraints are satisfied for each state-
action pair $({\bm{x}}^{\prime}_{t},{\bm{u}}^{\prime}_{t})$ at each
discretized time step $t$. If any constraint is violated, the steering
extension is immediately terminated [15].

One of the main advantages of this method is its efficiency in handling
multiple CBF constraints. By avoiding the need to solve QPs iteratively, the
algorithm significantly reduces computational overhead, and also circumvents
the issue of recursive feasibility that arises from the parameters of CBF and
QP [17]. In the following sections, we introduce two types of CBFs to ensure
that the resulting path is both collision-free (13) and visibility-aware (16).

###  IV-B HOCBF for Collision Avoidance

In this work, we employ the unicycle model to demonstrate the proposed
algorithm; however, the algorithm can be derived in the same way for more
complex system dynamics that include rotational motion. The unicycle model is
as follows:

$\displaystyle\dot{{\bm{x}}}$ $\displaystyle=\begin{bmatrix}\cos\theta&0\\\
\sin\theta&0\\\ 0&1\end{bmatrix}\begin{bmatrix}v\\\
\omega\end{bmatrix}=g({\bm{x}}){\bm{u}},$ (18)

where ${\bm{x}}=[x,y,\theta]^{\top}\in\mathcal{X}$ is the state vector and the
control input ${\bm{u}}=[v,\omega]^{\top}\in\mathcal{U}$ represents the
translational velocity $v$ and angular velocity $\omega$.

Hereafter, we abuse the notation and use ${\bm{x}}_{t}$ and ${\bm{u}}_{t}$ to
refer to ${\bm{x}}^{\prime}_{t}$ and ${\bm{u}}^{\prime}_{t}$ in the context of
the LQR-CBF-Steer function (Alg. 4). Let one of the obstacles with radius
$r_{\text{obs}}\in\mathbb{R}^{+}$ be located at
$(x_{\text{obs}},y_{\text{obs}})\in\mathbb{R}^{2}$. The constraint function
(7a) for collision avoidance, considering the robot radius $r_{\text{r}}$ and
the maximum tracking error $\epsilon$ mentioned in Sec. III, is defined as:

$h_{\textup{col}}({\bm{x}}_{t})=(x_{t}-x_{\text{obs}})^{2}+(y_{t}-y_{\text{obs}})^{2}-(r_{\text{obs}}+r_{\text{r}}+\epsilon)^{2}.$
(19)

Given the dynamics (18), the function $h_{\textup{col}}({\bm{x}}_{t})$ (19)
has relative degree $r_{b}=1$ for the control input $v$ and $r_{b}=2$ for the
control input $\omega$. To simplify the construction of the CBF constraint, we
fix the translational velocity $v$ to a constant value and only control the
angular velocity $\omega$ [16].

Then, the HOCBF constraint (9) with $r_{b}=2$ is [16]:

$\displaystyle\psi_{\textup{col}}({\bm{x}}_{t})$
$\displaystyle=2(x_{t}-x_{\text{obs}})v^{2}\cos^{2}{\theta_{t}}+2(y_{t}-y_{\text{obs}})v^{2}\sin^{2}{\theta_{t}}$
$\displaystyle+[2(y_{t}-y_{\text{obs}})v\cos{\theta_{t}}-2(x_{t}-x_{\text{obs}})v\sin{\theta_{t}}]\omega_{t}$
$\displaystyle+k_{1}L_{f}h_{\textup{col}}({\bm{x}}_{t})+k_{2}h_{\textup{col}}({\bm{x}}_{t})\geq
0,$ (20)

where $k_{1}$ and $k_{2}$ are positive constants that are designed to satisfy
the HOCBF condition (9), and

$L_{f}h_{\textup{col}}({\bm{x}}_{t})=2v(x_{t}-x_{\text{obs}})\cos{\theta_{t}}+2v(y_{t}-y_{\text{obs}})\sin{\theta_{t}}.$
(21)

Constraint (20) is then incorporated into the LQR-CBF-Steer function, ensuring
that every node added to the tree is collision-free, as illustrated in Fig.
1a.

Figure 1: Illustration of the CBF constraints in the LQR-CBF-Steer function.
The Integrator generates a sequence of intermediate states using the LQR gain
$K_{\text{lqr}}$. At each intermediate state, the CBFConstraints are checked
to ensure that the state satisfies the given constraints. If any constraint is
violated, the steering process is terminated, and the last satisfied state is
returned as the new node. (a) An example demonstrating the collision avoidance
HOCBF constraint (20). (b) An example demonstrating the visibility CBF
constraint (28). ${\bm{x}}_{\text{next},1}$ is already within the FOV at
${\bm{x}}_{t}$, satisfying the visibility constraint. For
${\bm{x}}_{\text{next},2}$, although it is outside the current FOV, it
satisfies the visibility constraint. ${\bm{x}}_{\text{next},3}$ might satisfy
the visibility constraint $h_{\text{vis}}({\bm{x}}_{t})\geq 0$ (27), but it
violates the visibility CBF constraint $\psi_{\textup{vis}}({\bm{x}}_{t})\geq
0$ (28), causing the steering to terminate.

###  IV-C Incorporating Visibility Constraint into CBF

The collision-free path planning discussed in Section IV-B assumes the
environment is fully known, i.e., $\mathcal{S}_{0}=\mathcal{S}$. However,
partially-known environments and limited perception, such as restricted FOV
and sensing range, can pose significant challenges in ensuring safe
navigation. To mitigate the risks posed by these limitations, we introduce a
visibility constraint formulated as a CBF to ensure that the generated path
allows for any newly-detected obstacles to be avoided by a local tracking
controller.

######  Definition 5 (Critical Point).

A critical point, denoted as ${\bm{x}}_{c}=[x_{c},y_{c},\theta_{c}]^{\top}$,
is the intersecting point between the frontier $\partial\mathcal{B}_{t}$,
which represents the boundary of the explored collision-free space
$\mathcal{B}_{t}$ up to the current time $t$, and the path from the start node
${\bm{x}}_{t}$ to the next node ${\bm{x}}_{\textup{next}}$. The angle
$\theta_{c}$ represents the direction from ${\bm{x}}_{t}$ to ${\bm{x}}_{c}$.

The critical point ${\bm{x}}_{c}$ represents the last sensed point along the
linearized path (as illustrated in Fig. 1b), beyond which the robot has no
updated information about potential unknown obstacles. The goal is to enforce
the robot to have seen the critical point within its FOV before reaching it.
To achieve this, we formulate a visibility constraint based on the time
required for the robot to rotate towards and observe the critical point,
denoted as $t_{\text{rot}}(\theta_{t})$, and the time needed to reach it,
denoted as $t_{\text{reach}}(x_{t},y_{t})$. The visibility constraint requires
that $t_{\text{reach}}(x_{t},y_{t})\geq t_{\text{rot}}(\theta_{t})$,
guaranteeing that the robot will have sufficient time to detect any potential
hidden obstacles at the critical point and allow the local tracking controller
$\pi$ to take appropriate actions if necessary.

####  IV-C1 Time-to-Reach

Given the robot’s current position ${\bm{r}}_{t}=[x_{t},y_{t}]^{\top}$ and the
critical point ${\bm{x}}_{c}$, let
$\Delta{\bm{r}}=\sqrt{(x_{t}-x_{c})^{2}+(y_{t}-y_{c})^{2}}$ denote the
Euclidean distance between them. To account conservatively for the robot
radius $r_{\text{r}}$ and the maximum tracking error $\epsilon$ of the local
controller $\pi$, we define $\Delta d=\Delta{\bm{r}}-r_{\text{r}}-\epsilon$.
The time-to-reach is:

$t_{\text{reach}}(x_{t},y_{t})=\frac{\Delta d}{v},$ (22)

where $v$ is the constant velocity assumed during the planning cycle as
mentioned in Sec. IV-B.

####  IV-C2 Time-to-Rotate

Given the current robot orientation $\theta_{t}$, the amount of angle that the
robot should rotate to observe the critical point using the onboard sensor is:

$\Delta\theta=|\theta_{t}-\theta_{c}|-\frac{\theta_{\text{fov}}}{2}.$ (23)

Using the cosine difference, we can rewrite $\Delta\theta$ as:

$\Delta\theta=\arccos{(\cos{\theta_{t}}\cos{\theta_{c}}+\sin{\theta_{t}}\sin{\theta_{c}})}-\frac{\theta_{\text{fov}}}{2}.$
(24)

Then, time-to-rotate $t_{\text{rot}}(\theta_{t})$ can be obtained as:

$t_{\text{rot}}(\theta_{t})=\frac{\Delta\theta}{\bar{\omega}_{t}},$ (25)

where $\bar{\omega}_{t}$ is the average angular velocity computed by
simulating the LQR controller $\mathcal{K}_{\text{lqr}}^{*}$ for a rotation of
$\Delta\theta$ starting from time $t$:

$\bar{\omega}_{t}=\frac{1}{T_{\omega}}\sum_{\tau=t}^{t+T_{\omega}}\omega_{\tau}^{\text{lqr}},$
(26)

with $\omega_{\tau}^{\text{lqr}}$ denoting the angular velocity generated by
the LQR controller $\mathcal{K}_{\text{lqr}}^{*}$ at time step $\tau$, and
$T_{\omega}$ being the total number of time steps required for the rotation.

####  IV-C3 Visibility CBF

Based on the visibility constraint $t_{\text{reach}}(x_{t},y_{t})\geq
t_{\text{rot}}(\theta_{t})$, we define the visibility CBF candidate
$h_{\textup{vis}}:\mathbb{R}^{3}\rightarrow\mathbb{R}$ as:

$h_{\textup{vis}}({\bm{x}}_{t})=t_{\text{reach}}(x_{t},y_{t})-t_{\text{rot}}(\theta_{t}).$
(27)

Then, the CBF constraint with $r_{b}=1$ (6) is defined as:

$\psi_{\textup{vis}}({\bm{x}}_{t})=\dot{h}_{\textup{vis}}({\bm{x}}_{t})+k_{3}\,h_{\textup{vis}}({\bm{x}}_{t})\geq
0,$ (28)

where $k_{3}$ is positive constant that is designed to satisfy the CBF
condition. The gradient $\dot{h}_{\textup{vis}}({\bm{x}}_{t})$ can be derived
as:

$\displaystyle\dot{h}_{\textup{vis}}({\bm{x}}_{t})$
$\displaystyle=L_{f}h_{\textup{vis}}({\bm{x}}_{t})+L_{g}h_{\textup{vis}}({\bm{x}}_{t}){\bm{u}}_{t}$
(29)
$\displaystyle=\frac{x_{t}-x_{c}}{\Delta{\bm{r}}}\cos\theta_{t}+\frac{y_{t}-y_{c}}{\Delta{\bm{r}}}\sin\theta_{t}$
$\displaystyle\quad-\frac{1}{\bar{\omega}_{t}}\left(\frac{\sin\theta_{t}\cos\theta_{c}-\cos\theta_{t}\sin\theta_{c}}{\sqrt{1-z_{t}^{2}}}\right)\omega_{t},$
(30)

where
$z_{t}=\cos{\theta_{t}}\cos{\theta_{c}}+\sin{\theta_{t}}\sin{\theta_{c}}$.

To demonstrate that this is a valid CBF, let us verify that (6) holds for all
${\bm{x}}\in\mathcal{X}$. $\frac{x_{t}-x_{c}}{\Delta{\bm{r}}}$,
$\frac{y_{t}-y_{c}}{\Delta{\bm{r}}}$, $\cos\theta_{t}$ and $\sin\theta_{t}$
range between $-1$ and $1$. Plus,
$\left|\frac{\sin\theta_{t}\cos\theta_{c}-\cos\theta_{t}\sin\theta_{c}}{\sqrt{1-z_{t}^{2}}}\right|\leq
1$. $\bar{\omega}_{t}$ is bounded since it is computed from the average LQR
control input, which converges to zero over time, and the largest control
input is bounded in the control space $\mathcal{X}$. Provided that all terms
in (28) are bounded and $h_{\textup{vis}}({\bm{x}}_{t})\geq 0$, we can choose
$k_{3}$ such that the condition (6) is satisfied for all
${\bm{x}}_{t}\in\mathcal{X}$. As a result, $h_{\textup{vis}}$ is a valid CBF.

The visibility CBF constraint $\psi_{\textup{vis}}({\bm{x}}_{t})\geq 0$ is
incorporated into the LQR-CBF-Steer function, as illustrated in Fig. 1b,
alongside the collision avoidance CBF constraint (20). By enforcing both
constraints during the steering process, the algorithm generates paths that
satisfy the feasibility of both safety conditions.

One issue in this formulation is that computing the location of the critical
point ${\bm{x}}_{c}$ whenever the LQR-CBF-Steer function is called, either
during tree expansion or rewiring, is computationally intensive since it
requires calculating the entire local collision-free space $\mathcal{B}_{t}$
[8]. To mitigate this issue in the implementation, we assume that
$\mathcal{B}_{t}$ forms a linear tube along the LQR trajectory between
${\bm{x}}_{t}$ and the parent node of ${\bm{x}}_{t}$. This allows us to
calculate the critical point ${\bm{x}}_{c}$ geometrically without the need to
track and store the entire $\mathcal{B}$ at each iteration, while still
providing a reasonable approximation of the critical point.

##  V LOCAL TRACKING CONTROLLER

(a)

(b)

(c)

(d)

Figure 2: Visualization of the global planning results generated by Baseline 3
and the proposed method in two environments. maxIter is set to 2000 and 3000
for Env. 1 and Env. 2, respectively. The blue and yellow squares represent the
start and goal position. The black circles represent the known obstacles. The
green lines depict the edges of the tree $\mathcal{E}$ appended during the
planning process. The red line depict the final reference path. The shaded
areas in gray represent the local collision-free set $\mathcal{B}_{t}$ that
the robot will sense while following the reference path.

After the global planning cycle described in Sec. IV, the local controller
$\pi$ tracks the waypoints ${\bm{x}}_{i}^{\text{ref}}\in\mathcal{P}$ enabling
the robot to navigate towards the goal. During navigation, the robot utilizes
its onboard sensor to create a map of the environment and detect obstacles on-
the-fly. The environment may contain unknown obstacles that were not
considered during the global planning phase. The local controller must ensure
that the robot remains within the true collision-free set $\mathcal{S}$ during
navigation, avoiding collisions with both known and unknown obstacles.

Two main strategies exist for designing safety-critical controllers in such
cases. The first strategy treats unmeasured space as free space and focuses on
avoiding previously known and newly detected obstacles. A common approach to
this problem is the use of CBFs to formally guarantee collision avoidance.
CBF-based quadratic programs (CBF-QPs) have emerged as a popular and
computationally efficient method for enforcing safety-critical control
objectives [2, 18]. However, a notable challenge emerges under limited sensing
capabilities. Specifically, if an obstacle is detected in close proximity to
the robot, particularly during rotational maneuvers, the CBF-QP may become
infeasible as there might not be a feasible solution to avoid the obstacle due
to the proximity and the robot’s dynamics.

The second strategy treats the unmeasured space as an unsafe set and only
controls the robot within the collision-free space sensed by the onboard
sensor. The GateKeeper algorithm [12] is a recent advance in this strategy,
acting as a safety filter between the planner and the low-level controller. If
the nominal trajectory enters the unmeasured space, GateKeeper executes a
backup controller to ensure the robot remains within the known collision-free
space. While this strategy prevents collisions with hidden obstacles if a
valid backup controller exists, the main challenge arises when the nominal
global planner is agnostic to the robot’s local sensing capabilities. In such
cases, GateKeeper frequently executes the backup controller and requires
global replanning based on the updated environmental information, which can
significantly degrade the system’s operational efficiency.

Assuming a safety-critical local controller $\pi$ interfacing with the global
planner, we can now prove the properties of the resulting path from our
proposed planner:

######  Theorem 3.

Let $p^{\textup{ref}}:\mathcal{X}\to\mathcal{P}$ be the global reference path
generated by the Visibility-Aware RRT* (Alg. 1), represented by a set of
waypoints $\mathcal{P}=\\{{\bm{x}}^{\textup{ref}}_{i}\\}_{i=0}^{N}$. If there
exists a local tracking controller $\pi$ with a maximum tracking error
$\epsilon$ that can track the trajectory between consecutive waypoints, then
the reference path $p^{\textup{ref}}$, when tracked by such controller $\pi$,
is both collision-free and visibility-aware.

###### Proof.

The trajectory $\bm{\sigma}_{i}$ connecting consecutive waypoints
${\bm{x}}^{\textup{ref}}_{i}$ and ${\bm{x}}^{\textup{ref}}_{i+1}$ is generated
by the LQR-CBF-Steer function. By construction, each trajectory satisfies both
collision avoidance (20) and visibility (28) CBF constraints, considering the
robot radius $r_{\text{r}}$ and the maximum tracking error $\epsilon$. Since
(20) is satisfied, the trajectories remain within the set
$\mathcal{C}_{\textup{col}}=\\{{\bm{x}}\in\mathcal{X}:h_{\textup{col}}({\bm{x}})\geq
0\\}$, where $h_{\textup{col}}$ is given by (19). This guarantees that the
local controller $\pi$ follows the trajectory $\bm{\sigma}_{i}$ with a maximum
tracking error of $\epsilon$ to remain collision-free (13). Moreover, since
(28) is met, the local controller $\pi$ forces the robot to move along
position trajectories $\bm{\sigma}_{i}$ from which the critical point
(according to Definition 5) is visible before it is reached. This guarantees
that the local controller $\pi$ can detect and avoid hidden obstacles that, in
the worst-case scenario, are placed at the boundary of the local safe set
$\mathcal{B}_{t}$. Consequently, the local controller $\pi$ can keep the robot
within the known local collision-free sets $\mathcal{B}_{t}$ while tracking
the reference path $p^{\textup{ref}}$, satisfying the visibility-aware
condition (16). ∎

##  VI RESULTS

In this section, we evaluate the proposed visibility-aware RRT* through
simulation results in two distinct environments. We demonstrate that our
approach can address the challenges faced by both CBF-QP and GateKeeper
methods, providing a safe and efficient solution for navigating in
environments with unknown obstacles, which visibility-agnostic planners fail
to achieve.

The first environment (Env. 1) spans an area of 15 m $\times$ 15 m, while the
second one (Env. 2) covers a larger area of 35 m $\times$ 30 m. Both
environments contain known obstacles that prevent the robot from taking a
straight-line path to the goal. The onboard sensor’s FOV $\theta_{\text{fov}}$
and the sensing range $l_{\text{range}}$ are set to $70^{\circ}$ and 3 m,
respectively, unless otherwise specified.

(a)

(b)

Figure 3:  Visualization of the CBF-QP experiments for Env. 1 with a
$45^{\circ}$ FOV. The reference waypoints (green dots) correspond to the paths
in Fig. 2(a) and Fig. 2(b). The blue shaded areas represent the actual local
collision-free set $\mathcal{B}_{t}$ collected from the onboard sensor. The
orange circles are the hidden obstacles $\mathcal{H}$ and the red dots
indicate the detection points of these hidden obstacles. (a) The CBF-QP
becomes infeasible when it detects the hidden obstacle, as the robot is
already within an unsafe distance to perform collision avoidance. (b) Tracking
the reference path from our method, the CBF-QP successfully avoids the hidden
obstacle and reaches the goal.

(a)

(b)

Figure 4: Visualization of the GateKeeper experiments for Env. 2. The reference waypoints correspond to the paths in Fig. 2(c) and Fig. 2(d). Red shaded areas depict the subsequent nominal trajectories at the robot’s position. (a) The GateKeeper executes a stop command as the nominal trajectory is deemed unsafe, where the trajectory lies outside of the local collision-free set $\mathcal{B}_{t}$. Consequently, global replanning is required from the current position. (b) By tracking our reference path, the GateKeeper successfully navigates to the goal without violating visibility constraint. TABLE I: (left) Collision rates of the CBF-QP controller tracking reference paths generated by different planners in two environments with varying FOV angles. (right) Percentage of reference paths that triggered the execution of the backup controller in the GateKeeper algorithm, indicating the need for global replanning from scratch. Note that Baseline 1 is a deterministic planner, there is only one reference path tested. Controller | CBF-QP | GateKeeper
---|---|---
Env. | Env. 1 | Env. 2 | Env. 1 | Env. 2
FOV | $45^{\circ}$ | $70^{\circ}$ | $45^{\circ}$ | $70^{\circ}$ | $70^{\circ}$ | $70^{\circ}$
Baseline 1 | 100% | 0% | 0% | 0% | 100% | 100%
Baseline 2 | 5% | 3% | 13% | 9% | 23% | 32%
Baseline 3 | 4% | 0% | 1% | 1% | 18% | 13%
Ours | 0% | 0% | 0% | 0% | 1% | 0%

We compare four planning methods in our experiments: (i) The A* algorithm with
an obstacle distance cost [19] (Baseline 1), (ii) the standard LQR-RRT*
algorithm [13] (Baseline 2), (iii) our proposed method without the visibility
constraint (Baseline 3), and (iv) the propsed visibility-aware RRT*. To
highlight the effect of the visibility constraint, we illustrate the reference
paths generated by Baseline 3 and our method in Fig. 2. The results
demonstrate that our method maintains fewer nodes in the tree compared to
Baseline 3, as it only extends the nodes that satisfy the visibility
constraint. This selective expansion leads to a more computationally efficient
exploration of the configuration space [20].222https://github.com/tkkim-
robot/visibility-rrt Furthermore, the final reference path generated by our
method tends to take shallower turns around the known obstacles. This behavior
complements the local tracking controller by providing a path that allows it
to detect and avoid potential hidden obstacles.

To demonstrate the effectiveness of the generated paths when interfacing with
the local tracking controllers, we employ two types of controllers to track
the paths as discussed in Sec. V: a CBF-QP controller and the GateKeeper
algorithm. We employ the dynamic unicycle model, which takes translational
acceleration as input instead of velocity, in our tracking controllers to
effectively capture the realistic dynamics and constraints of the robot [21].

First, we evaluate the safety perspective of the reference path when it
interfaces with the CBF-QP controller. We formulate the CBF-QP as in [21],
using the dynamic unicycle model, incorporating a CBF constraint for collision
avoidance. In this experiment, we test two different FOV angles: $45^{\circ}$
and $70^{\circ}$. For Baseline 1, which is a deterministic planner, we run the
algorithm once and for Baselines 2-3 and our proposed method, which are
sampling-based planners, we generate 100 reference paths for each setting.
After the planning phase, we place unknown obstacles near the known obstacles.
Then, we evaluate whether the CBF-QP controller can successfully navigate the
robot to the goal without collisions. One example of this experiment is shown
in Fig. 3. The results, summarized in Table I, demonstrate that even with a
safety-critical controller using a valid CBF for collision avoidance, safety
cannot be guaranteed if the planner does not consider the robot’s perception
limitations. This is because the robot detects the hidden obstacle too late,
leaving no feasible solution for the CBF-QP to avoid the obstacle. This aspect
becomes more prominent as the sensing capability becomes more limited. In
contrast, the visibility-aware paths generated by our method can effectively
complement the local tracking controller, ensuring safe navigation even in the
presence of limited perception.

Next, we investigate the efficiency of the reference paths when interfacing
with the GateKeeper algorithm [12]. GateKeeper assesses the safety of the
nominal trajectory and decides whether to follow it or execute a backup
controller, which, in this case, is a stop command. If the next nominal
trajectory is deemed unsafe, meaning the robot might not be able to completely
stop within the known collision-free space given its control authority,
GateKeeper opts not to follow the nominal trajectory. Instead, it executes the
backup controller and subsequently plans for global path again. We evaluate
the percentage of the same set of reference paths that trigger the execution
of the backup controller, indicating that the path was inefficient for
navigation. One example of this experiment is shown in Fig. 4. The results,
presented in Table I, show that visibility-aware RRT* significantly reduces
the need for replanning compared to the baseline planners. This demonstrates
that the proposed method generates reference paths that are more efficient for
navigation by minimizing the frequency of replanning required when used with
controllers that treat unmeasured space as an unsafe set.

##  VII CONCLUSION

In this paper, we considered that safety may not be guaranteed for robots with
limited perception capabilities navigating in unknown environments, even when
using safety-critical controllers. To address this challenge, we introduced
the Visibility-Aware RRT* algorithm, which generates a global reference path
that is both collision-free and visibility-aware, enabling the tracking
controller to avoid hidden obstacles at all times. We introduced the
visibility CBF to encode the robot’s visibility limitations, utilizing its
constraint as a termination criterion during the steering process, alongside a
standard collision avoidance CBF constraint. We also proved that our algorithm
ensures the satisfaction of visibility constraints by the local controller,
and we evaluated the performance of the integrated path planner with two
different controllers via simulations. The results indicated that our method
outperformed all other compared baselines in both safety and efficiency
aspects.

## References

  * [1] J. G. da Silva Neto, P. J. da Lima Silva, F. Figueredo, J. M. X. N. Teixeira, and V. Teichrieb, “Comparison of RGB-D sensors for 3D reconstruction,” in _Symposium on Virtual and Augmented Reality (SVR)_ , 2020, pp. 252–261.
  * [2] A. D. Ames, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs with application to adaptive cruise control,” in _IEEE Conference on Decision and Control_ , 2014, pp. 6271–6278.
  * [3] I. Spasojevic, V. Murali, and S. Karaman, “Perception-aware time optimal path parameterization for quadrotors,” in _IEEE International Conference on Robotics and Automation_ , 2020, pp. 3213–3219.
  * [4] B. Zhou, J. Pan, F. Gao, and S. Shen, “RAPTOR: Robust and Perception-Aware Trajectory Replanning for Quadrotor Fast Flight,” _IEEE Transactions on Robotics_ , vol. 37, no. 6, pp. 1992–2009, 2021.
  * [5] V. Murali, I. Spasojevic, W. Guerra, and S. Karaman, “Perception-aware trajectory generation for aggressive quadrotor flight using differential flatness,” in _American Control Conference (ACC)_ , 2019, pp. 3936–3943.
  * [6] R. M. Bena, C. Zhao, and Q. Nguyen, “Safety-Aware Perception for Autonomous Collision Avoidance in Dynamic Environments,” _IEEE Robotics and Automation Letters_ , vol. 8, no. 12, pp. 7962–7969, 2023.
  * [7] X. Chen, Y. Zhang, B. Zhou, and S. Shen, “APACE: Agile and Perception-Aware Trajectory Generation for Quadrotor Flights,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2024.
  * [8] G. Oriolo, M. Vendittelli, L. Freda, and G. Troso, “The SRT method: randomized strategies for exploration,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2004, pp. 4688–4694.
  * [9] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada, “Control Barrier Functions: Theory and Applications,” in _European Control Conference (ECC)_ , 2019, pp. 3420–3431.
  * [10] W. Xiao and C. Belta, “Control Barrier Functions for Systems with High Relative Degree,” in _IEEE Conference on Decision and Control (CDC)_ , 2019, pp. 474–479.
  * [11] Q. Nguyen and K. Sreenath, “Exponential Control Barrier Functions for enforcing high relative-degree safety-critical constraints,” in _American Control Conference (ACC)_ , 2016, pp. 322–328.
  * [12] D. Agrawal, R. Chen, and D. Panagou, “gatekeeper: Online Safety Verification and Control for Nonlinear Systems in Dynamic Environments,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2023, pp. 259–266.
  * [13] A. Perez, R. Platt, G. Konidaris, L. Kaelbling, and T. Lozano-Perez, “LQR-RRT*: Optimal sampling-based motion planning with automatically derived extension heuristics,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2012, pp. 2537–2542.
  * [14] S. Karaman and E. Frazzoli, “Sampling-based algorithms for optimal motion planning,” _The International Journal of Robotics Research_ , vol. 30, no. 7, pp. 846–894, 2011.
  * [15] G. Yang, M. Cai, A. Ahmad, A. Prorok, R. Tron, and C. Belta, “LQR-CBF-RRT*: Safe and Optimal Motion Planning,” in _arXiv preprint arXiv:2304.00790_ , 2023.
  * [16] G. Yang, B. Vang, Z. Serlin, C. Belta, and R. Tron, “Sampling-based Motion Planning via Control Barrier Functions,” in _International Conference on Automation, Control and Robots_ , 2019, pp. 22–29.
  * [17] H. Parwana and D. Panagou, “Recursive Feasibility Guided Optimal Parameter Adaptation of Differential Convex Optimization Policies for Safety-Critical Systems,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2022, pp. 6807–6813.
  * [18] M. Jankovic and M. Santillo, “Collision Avoidance and Liveness of Multi-agent Systems with CBF-based Controllers,” in _IEEE Conference on Decision and Control (CDC)_ , 2021, pp. 6822–6828.
  * [19] T. Kim, S. Lim, G. Shin, G. Sim, and D. Yun, “An Open-Source Low-Cost Mobile Robot System With an RGB-D Camera and Efficient Real-Time Navigation Algorithm,” _IEEE Access_ , vol. 10, pp. 127 871–127 881, 2022.
  * [20] C. Tao, H. Kim, H. Yoon, N. Hovakimyan, and P. Voulgaris, “Control Barrier Function Augmentation in Sampling-based Control Algorithm for Sample Efficiency,” in _American Control Conference (ACC)_ , 2022, pp. 3488–3493.
  * [21] H. Parwana, M. Black, B. Hoxha, H. Okamoto, G. Fainekos, D. Prokhorov, and D. Panagou, “Feasible Space Monitoring for Multiple Control Barrier Functions with application to Large Scale Indoor Navigation,” in _arXiv preprint arXiv.2312.07803_ , 2023.