# Stars on the edge: Galactic tides and the outskirts of the Sculptor dwarf
spheroidal

Federico Sestito1, Joel Roediger2, Julio F. Navarro1, Jaclyn Jensen1, Kim A.
Venn1, Simon E. T. Smith1, Christian Hayes2, and Alan W. McConnachie2,1
1 Department of Physics and Astronomy, University of Victoria, PO Box 3055,
STN CSC, Victoria BC V8W 3P6, Canada
2 NRC Herzberg Astronomy & Astrophysics, 5071 West Saanich Road, Victoria, BC
V9E 2E7, Canada
Email<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The formation of "stellar halos" in dwarf galaxies have been discussed in
terms of early mergers or Galactic tides, although fluctuations in the
gravitational potential due to stellar feedback is also a possible candidate
mechanism. A Bayesian algorithm is used to find new candidate members in the
extreme outskirts of the Sculptor dwarf galaxy. Precise metallicities and
radial velocities for two distant stars are measured from their spectra taken
with the Gemini South GMOS spectrograph. The radial velocity, proper motion
and metallicity of these targets are consistent with Sculptor membership. As a
result, the known boundary of the Sculptor dwarf extends now out to an
elliptical distance of $\sim 10$ half-light radii, which corresponds to a
projected physical distance of $\sim 3$ kpc. As reported in earlier work, the
overall distribution of radial velocities and metallicities indicate the
presence of a more spatially and kinematically dispersed metal-poor population
that surrounds the more concentrated and colder metal-rich stars. Sculptor’s
density profile shows a "kink" in its logarithmic slope at a projected
distance of $\sim 25$ arcmin (620 pc), which we interpret as evidence that
Galactic tides have helped to populate the distant outskirts of the dwarf. We
discuss further ways to test and validate this tidal interpretation for the
origin of these distant stars.

###### keywords:

stars: abundances – stars: Population II – galaxies : formation – galaxies:
dwarf – galaxies: individual: Sculptor – galaxies: evolution

††pubyear: 2023††pagerange: Stars on the edge: Galactic tides and the
outskirts of the Sculptor dwarf spheroidal–References

##  1 Introduction

The dwarf satellites orbiting the Milky Way are among the oldest and most
metal-poor galaxies that are known to exist (e.g., Tolstoy, Hill & Tosi,
2009). These galaxies have a broad range of masses, sizes, and luminosities. A
luminosity threshold of L$=10^{5}$ L⊙ was proposed by Simon (2019) to
differentiate between classical dwarf galaxies (such as Fornax, Sculptor, and
Ursa Minor) and the least massive ultra-faint dwarf galaxies (UFDs, also see
Bullock & Boylan-Kolchin, 2017). The UFDs are situated at the low-mass end of
the hierarchical formation process, with just enough mass to produce very
metal-poor stars (VMP, with ${\rm[Fe/H]}{}\leq-2.0$, Simon, 2019). The dwarf
satellites are a prime target to comprehend the characteristics of dark matter
and the early stages of our galaxy’s formation (e.g., Bullock & Boylan-
Kolchin, 2017).

The presence of extended stellar halos around dSphs is of particular interest.
In the context of the hierarchical formation of $\Lambda-$Cold Dark Matter
($\Lambda-$CDM) cosmology (e.g., White & Rees, 1978; Frenk et al., 1988;
Navarro, Frenk & White, 1997), stellar haloes are built primarily through the
accretion and disruption of smaller systems. While this process is clearly at
work in large galaxies like the Milky Way (MW), the possible presence of
stellar haloes around dSphs remains inconclusive (e.g., Deason et al., 2022,
and references therein). In $\Lambda-$CDM, the fraction of stars accreted
through mergers is thought to become negligible in the regime of dwarfs,
suggesting, at face value, that dwarfs should generally lack stellar halos
(e.g., Genel et al., 2010; Moster, Naab & White, 2013).

Faint dwarf galaxies have shallow gravitational potentials, and are therefore
extremely susceptible to a number of processes that may act to influence their
individual morphologies and perhaps lead to the formation of an extended
envelope of stars (e.g., Higgs et al., 2021, and references therein). Internal
processes include star formation and the subsequent stellar feedback (e.g.,
El-Badry et al., 2018); external processes include mergers (e.g., Deason,
Wetzel & Garrison-Kimmel, 2014), ram pressure stripping (e.g., Grebel,
Gallagher & Harbeck, 2003) and stirring (e.g., Kazantzidis et al., 2011),
tidal interaction (e.g., Fattahi et al., 2018), and reionization (e.g.,
Wheeler et al., 2019). The outskirts of a dwarf satellite galaxy is fertile
ground to look for signatures of these internal and external processes.

Very recently, new candidate member stars have been reported in the extreme
outskirts of several dwarf galaxies (e.g., McConnachie & Venn, 2020a; Chiti et
al., 2021; Filion & Wyse, 2021; Longeard et al., 2022; Yang et al., 2022;
Jensen et al., 2023, in prep.; Sestito et al., 2023b; Waller et al., 2023).
Chiti et al. (2021) spectroscopically identified members as far out as $\sim$9
half-light radii ($r_{h}$), or physical distances up to 1 kpc, away from the
centre of Tucana II. Their chemo-dynamical investigation showed that the outer
envelope of Tucana II is more metal-poor than the central regions, a result
that suggests that Tucana II was probably formed by an early merger which
dispersed the old stellar component and brought gas to the inner regions to
fuel the formation of subsequent generations of stars (Benítez-Llambay et al.,
2016).

Recently, Jensen et al. (2023, in prep.) improved a Bayesian algorithm,
developed by McConnachie & Venn (2020a), to estimate the probability that a
star in the vicinity of a dwarf galaxy is a member of the dwarf. Their method
combines the full astrometric and photometric data from Gaia Early Data
Release 3 (EDR3, Gaia Collaboration et al., 2021). The algorithm is able to
identify as members the distant stars found by Chiti et al. (2021) around
Tucana II (Jensen et al., 2023, in prep.), and it has been used to hunt for an
extended stellar halo in Coma Berenices, Ursa Major I, Boötes I ($\sim 6\ \
r_{h}$, Waller et al., 2023), and Ursa Minor ($\sim 12\ \ r_{h}$, Sestito et
al., 2023b). Both Waller et al. (2023) and Sestito et al. (2023b) discussed
that the outermost stars in these systems may have formed in the central
regions, then moved out through tidal stripping and/or supernova feedback.

Other algorithms have also been applied to search for extended stellar haloes
in faint systems. Yang et al. (2022) report a new halo component in Fornax,
representing about 10 percent of the total stellar mass; Longeard et al.
(2022) combined Gaia EDR3 data with photometric metallicites from the Pristine
survey (Starkenburg et al., 2017) to discover 27 new members in the outskirt
of Boötes I. In the same system, Filion & Wyse (2021) found several blue
horizontal branch stars beyond its tidal radius.

In this work, candidates identified by this extremely efficient Bayesian
algorithm (Jensen et al., 2023, in prep.) are used to confirm the presence of
stars in the outermost region of the Sculptor dwarf galaxy (Scl). The two
outermost candidate member stars at $\sim 10$ half-light radii are observed
with low-resolution spectroscopy. Their radial velocities and metallicities
confirm their membership to Sculptor.

Sculptor is the first faint MW satellite discovered (Shapley, 1938) after the
Magellanic Clouds. Since then, the number of known faint dwarf galaxies has
grown tremendously thanks to new observational techniques and survey, e.g.,
the Sloan Digital Sky Survey (SDSS, York et al., 2000), the Panoramic Survey
Telescope And Rapid Response System (Pan-STARRS, Chambers et al., 2016), the
Dark Energy Survey (DES, The Dark Energy Survey Collaboration, 2005), the DESI
Legacy Imaging Surveys (Dey et al., 2019), the Pan-Andromeda Archaeological
Survey (PAndAS, e.g., Martin et al., 2009; McConnachie et al., 2009), and the
DECam Local Volume Exploration Survey (DELVE, Drlica-Wagner et al., 2021).
Many spectroscopic campaigns have targeted stars in Sculptor. In particular,
the Dwarf Abundance and Radial velocity Team (DART, Tolstoy et al., 2006;
Battaglia et al., 2008) observed more than 600 members measuring precise
metallicities and radial velocities, showing that the system spans 3 dex in
metallicity.

High-resolution spectroscopic observations probing the chemical properties of
Sculptor are more scarce. There are only 12 stars with ${\rm[Fe/H]}<-2.5$
observed so far (Frebel, Kirby & Simon, 2010; Tafelmeyer et al., 2010;
Starkenburg et al., 2013; Jablonka et al., 2015; Simon et al., 2015;
Skúladóttir et al., 2021). The first ultra-metal-poor star (UMP,
${\rm[Fe/H]}\leq-4.0$) observed at high-resolution in a dwarf galaxy is a
Sculptor member, as recently discovered by Skúladóttir et al. (2021) (vs. 42
detected UMPs in the Milky Way, Sestito et al., 2019). The authors discuss
that the atmosphere of this ancient star carries the imprints of a hypernova
generated from a metal-free star, also known as First star or Population III.
This finding is in agreement with Ishigaki et al. (2018), in which they
discuss that the the majority of the primordial supernovae would have masses
$<40{\rm\,M_{\odot}}$. However, Hartwig et al. (2023) showed that only $\sim
40$ percent of the stars at that metallicity (${\rm[Fe/H]}=-4.11$, Skúladóttir
et al., 2021) would be mono-enriched by primordial supernovae. The chemical
evolution of Sculptor is similar to other dwarf galaxies. At very low-
metallicities, Sculptor is dominated by contributions from supernovae type II
(e.g., Jablonka et al., 2015; Skúladóttir et al., 2021), while at
${\rm[Fe/H]}\geq-1.9$ indications of supernovae type Ia are observable in the
distribution of chemical abundances (e.g., Venn et al., 2004; Tolstoy, Hill &
Tosi, 2009; Hill et al., 2019). Star formation in Sculptor has been modelled
to have lasted between 6 to 8 Gyr (e.g., Fenner et al., 2006; Revaz &
Jablonka, 2018).

This paper is divided as follow. Section 2 reports the target selection, the
observational setup, and the spectra reduction of the two targets. Section 3
describes the inference of the stellar parameters. Section 4 outlines the
measurement of the radial velocity and the metallicity from the near infrared
Ca Triplet lines. The discussion on the membership of the targets and the
presence of tidal perturbations in Sculptor are reported in Section 5.
Conclusions are summarised in Section 6.

##  2 Data

###  2.1 Target selection

Candidate member stars for spectroscopic follow-up are firstly selected using
the algorithm described in Jensen et al. (2023, in prep., hereafter J23). The
algorithm, much like its predecessor described in McConnachie & Venn (2020a,
hereafter MV20), is designed to search for member stars in a given dwarf
galaxy by determining the probability of membership to the satellite. The
probability of being a satellite member, Psat, is determined using three
likelihoods based on the dwarf’s (1) colour-magnitude diagram (CMD), (2)
systemic proper motion, and (3) the radial distance from the center of the
satellite, using relevant data from Gaia EDR3 (Gaia Collaboration et al.,
2021). The stellar density of a dwarf can often sufficiently be summarized by
a single exponential profile (see MV20). However, motivated to search for
dwarfs which may have evidence of tidal features or extended stellar haloes,
the spatial likelihood in J23 assumes that each dwarf may host a secondary,
extended, and lower density, outer profile. Only a handful of systems are
found to host a secondary outer profile (see J23), a few of which are already
known to be actively tidally disrupting (e.g., Bootes III and Tucana III).
Sculptor is a system for which a secondary outer profile is observed,
indicating either an extended stellar halo or tidal features.

The most radially distant (9.7 and 10.1 $r_{h}$), relatively bright (G $\sim
17.6-17.9$ mag) candidates of Sculptor are selected. J23 demonstrated that the
purity of candidate members is extremely high for P${}_{\rm{sat}}>0.40$ and
the majority of foreground stars has P${}_{\rm{sat}}<0.10$. Therefore, the
algorithm is excellent at removing foreground contaminants even at such low
probabilities. As a matter of fact, these stars are also listed as likely Scl
members by Qi et al. (2022), with a probability $>98$ percent. Thus far, this
algorithm has been proved to be very efficient in finding new members in the
extreme outskirts of Coma Berenices, Ursa Major 1, Bootes 1 (Waller et al.,
2023), and in Ursa Minor (Sestito et al., 2023b).

For the most distant candidate, Target 2, a measurement of its radial velocity
is already present in the literature. Westfall et al. (2006) measured radial
velocities of stars in the direction of Scl. They found a radial velocity (RV)
compatible with Scl ($\sim 119{\rm\,km\ s^{-1}}$). Given its large distance
from Scl centre, the authors had no reason to suspect Target 2 was a member,
therefore, no metallicity was estimated. The main properties of Sculptor and
our two targets are reported in Tables 1 and 2.

Figure 1 shows the position of our two targets (black diamond and star
markers), $\sim 7000$ candidate members selected from J23 (olive diamonds),
and other known Sculptor members (lime green and hot pink markers) in
projected sky coordinates, on colour-magnitude diagram, and in proper motion
space. The selection algorithm is also applied to stars within the Dwarf
Abundance and Radial velocity Team (DART, Tolstoy et al., 2006; Battaglia et
al., 2008) and to the APOGEE data release 17 (DR17, Abdurro’uf et al., 2022).
From the initial sample of DART (656 stars), stars with low probability
membership ($<40$ percent, 118 stars) have been removed. Then, for the cleaned
DART sample in common with APOGEE DR17, radial velocities and metallicities
from the former are used, since the metallicity grid in APOGEE stops at
[Fe/H]$=-2.4$ and several stars have metallicity below that threshold
according to DART. We also note that RV in APOGEE DR17 are in agreement with
DART even at ${\rm[Fe/H]}<-2.4$. The total sample with reliable radial
velocities and metallicities is composed of 617 stars, 538 from DART and 79
from APOGEE DR17.

Figure 1: Sculptor as seen in Gaia EDR3. All panels: Targets 1 and 2 are
marked with a diamond and a star marker, respectively. Olive diamonds are
candidate members from Jensen et al. (2023, in prep.). Lime green circles and
hot pink squares are Sculptor stars from the cleaned DART sample and APOGEE
DR17, respectively. MW foreground stars are shown with small grey dots. These
are selected from Gaia EDR3 in the direction of Sculptor and within the field
of view of the $\eta-\xi$ panel. Left panel: Projected sky coordinates and
projected distance from Sculptor’s centre. The orange ellipses denote
elliptical distances of 3, 5, 7, and 9 $r_{h}$. The arrow points in the
direction of Sculptor’s proper motion. Central panel: Colour-magnitude
diagram. Dark green dashed line is a Padova isochrone at ${\rm[Fe/H]}=-2.0$
and age of 12 Gyr (Bressan et al., 2012). Right panel: Proper motion space.
Table 1: Galactic parameters of Sculptor. The coordinates $\alpha,\delta$, the
mean metallicity, the mean radial velocity, the dispersion in velocity, the
heliocentric distance D⊙, the distance modulus, the ellipticity, the position
angle $\phi$, and the half-light radius rh in arcmin and pc, the Gaia EDR3
mean proper motion, and the dynamical mass are reported with the respective
references. (a) refers to McConnachie (2012), (b) to McConnachie & Venn
(2020a), and (c) to McConnachie & Venn (2020b).

Property | Value | Reference
---|---|---
$\alpha$ | 15.03917 deg | (b)
$\delta$ |  $-33.70917$ deg | (b)
$\overline{{\rm[Fe/H]}}$ | $-1.68\pm 0.01$ | (b)
$\overline{\rm{RV}}$ |  $111.4$ km s-1 | (a)
$\sigma_{v}$ |  $9.2\pm 1.4$ km s-1 | (b)
D⊙ |  $86\pm 6$ kpc | (a)
Dmod |  $19.67\pm 0.14$ mag | (a)
ellipticity | $0.37\pm 0.01$ | (b)
$\phi$ |  $94\pm 1$ deg | (b)
rh |  $12.33\pm 0.05$ arcmin | (b)
rh |  $308\pm 20$ pc | (b)
$\mu_{\alpha}\rm{cos}\delta$ |  $0.099\pm 0.002$ mas yr-1 | (c)
$\mu_{\delta}$ |  $-0.160\pm 0.002$ mas yr-1 | (c)
Mdyn($\leq$rh) | $14\times 10^{6}{\rm\,M_{\odot}}$ | (a)

Table 2: The Gaia EDR3 source ID, the coordinates $(\alpha,\delta)$, the
projected coordinates $(\xi,\eta)$, the elliptical radius distance
$r_{\rm{ell}}$ in $r_{h}$ unit, the probability to be a member from Jensen et
al. (2023, in prep.), the Gaia EDR3 photometry G and BP$-$RP, and the
reddening AV from Schlafly & Finkbeiner (2011) are reported for each target.

Target | source id | $\alpha$ | $\delta$ | $\xi$ | $\eta$ | $r_{\rm{ell}}$ | $P_{\rm{sat}}$ | G | BP$-$RP | AV
---|---|---|---|---|---|---|---|---|---|---
|  | (deg) | (deg) | (deg) | (deg) | ($r_{h}$) |  | (mag) | (mag) | (mag)
Target 1 | 5006419626331394048 | $12.67539$ | $-33.46838$ | $-1.97246$ | $0.21835$ | $9.67$ | $0.46$ | $17.62$ | $1.31$ | 0.0419
Target 2 | 5026130884816022016 | $16.84351$ | $-34.66460$ | $1.48453$ | $-0.96882$ | $10.06$ | $0.45$ | $17.86$ | $1.29$ | 0.0439

###  2.2 GMOS observations and reduction

Targets are observed with the Gemini Multi-Object Spectrograph (GMOS, Hook et
al., 2004; Gimeno et al., 2016) at the Gemini South with a Fast Turnaround
program, GS-2022B-FT-205 (P.I. F. Sestito). The R831 grating centred at the
855 nm Ca ii Triplet line, with a slit width of 1.0 arcsec, is used for these
observations. Each target is observed with 4 exposures of 900s each to obtain
a signal-to-noise ratio per spectral pixel of $\sim$50 in the Ca ii Triplet
region, in order to measure the RV with an uncertainty around $\sim 10$
${\rm\,km\ s^{-1}}$ and the equivalent width (EW) of the Ca ii Triplet with a
sufficient precision for reaching an uncertainty of $\sim 0.3$ dex in the
metallicity estimate. The latter is inferred through various Ca ii T EW -
[Fe/H] relations, as described in Section 4.

Spectra have been reduced with the standard routines for GMOS data in the
Gemini IRAF package111http://www.gemini.edu/observing/phase-iii/understanding-
and-processing-data/data-processing-software/gemini-iraf-general (Tody, 1986,
1993). Spectra have been corrected for bias and flat-fields, and the
wavelength has been calibrated using the CuAr emission lamp. Heliocentric
motion has been corrected in the observations. The 4 exposures have been
stacked using GemCombine. The flux of the combined spectra have been
normalised to one. The log of the observations is reported in Table 3.

Figure 2 displays the the Ca ii region of the combined and normalised GMOS
spectra of Target 1 and 2.

Figure 2: Normalised GMOS spectra in the Ca T region. Table 3: Total exposure
time, number of exposures, signal-to-noise ratio (SNR) measured at the Ca ii T
region, and the observation dates are reported for each target. The SNR is
calculated as the ratio between the median flux and its standard deviation in
the region [8550,8650] Å.

Target | $t_{\rm exp}$ | Nexp | SNR | Obs. date
---|---|---|---|---
| (s) |  | @Ca ii T | YY/MM/DD
Target 1 | 3600 | 4 | 66 | 22/11/15
Target 2 | 3600 | 4 | 78 | 22/10/30

##  3 Stellar parameters

To estimate the effective temperatures, the Mucciarelli, Bellazzini & Massari
(2021) colour-temperature relationship for giants is used in this work. This
is based on the Infrared flux method from González Hernández & Bonifacio
(2009) and adapted to Gaia EDR3 photometry. The input parameters are the Gaia
EDR3 (BP-RP) de-reddened colour and a metallicity estimate. The 2D Schlafly &
Finkbeiner (2011) map222https://irsa.ipac.caltech.edu/applications/DUST/ has
been used to correct the photometry for extinction333To convert from the
E(B-V) map to Gaia extinction coefficients, the $\rm A_{V}/E(B-V)=3.1$
(Schultz & Wiemer, 1975) and the $\rm A_{G}/A_{V}=0.85926$, $\rm
A_{BP}/A_{V}=1.06794$, $\rm A_{RP}/A_{V}=0.65199$ relations (Marigo et al.,
2008; Evans et al., 2018) are used.. A value of ${\rm[Fe/H]}=-2.0\pm 0.5$ is
adopted as input metallicity, compatible with the distribution of [Fe/H] in
Sculptor.

Surface gravities were inferred from the Stefan-Boltzmann
equation444$L_{\star}=4\pi R_{\star}^{2}\sigma T_{\star}^{4}$; the radius of
the star can be calculated from this equation, then the surface gravity is
inferred assuming the mass.. This calculation requires the effective
temperature, the distance of the object (see Table 1), the Gaia EDR3 G de-
reddened photometry, and the bolometric corrections on the flux (Andrae et
al., 2018). A Monte Carlo algorithm has been applied to the input parameters
with their uncertainties to estimate the uncertainties on the stellar
parameters. The input quantities were drawn from a Gaussian distribution,
except for the stellar mass. The latter is treated with a flat prior from 0.5
to 0.8 ${\rm\,M_{\odot}}$, which is consistent with the mass of very metal-
poor stars. The mean uncertainty on the effective temperature is $\sim 72$ K,
while on the surface gravity is $\sim 0.08$ dex, both mainly driven by the
distance’s uncertainty.

Stellar parameters from these methods have been shown to be in agreement with
spectroscopic methods applied to very metal-poor stars (e.g., Kielty et al.,
2021; Lardo et al., 2021; Sestito et al., 2023a; Waller et al., 2023). The
stellar parameters are reported in Table 4.

Table 4: Effective temperatures and surface gravities of the targets. Target | Teff | log g
---|---|---
| (K) |
Target 1 | $4535\pm 75$ | $1.11\pm 0.07$
Target 2 | $4583\pm 70$ | $1.24\pm 0.08$

##  4 Metallicity and radial velocity estimation

###  4.1 Radial velocity measurement

The radial velocity has been calculated by cross-correlating the combined
normalised spectra with a synthetic spectrum using the fxcor routine in IRAF.
The synthetic spectrum has been created from synth MOOG (Sneden, 1973) using a
spectral line list from linemake555https://github.com/vmplacco/linemake
(Placco et al., 2021). MARCS model atmospheres (Gustafsson et al., 2008; Plez,
2012) are adopted with the same stellar parameters as the observed stars. The
synthetic spectrum has been created at the same resolution of the GMOS R831
grating.

To test for possible RV offsets due to problems in the wavelength solution,
telluric absorption lines have been extracted from the spectra of our targets
and cross-correlated. There is no RV offset between the spectra of the two
targets within the dispersion of GMOS spectrograph ($\sim 13{\rm\,km\
s^{-1}}$).

###  4.2 Metallicities from Ca ii Triplet

Multiple calibrations to measure [Fe/H] from the near-infrared Ca ii Triplet
lines ($\lambda\lambda 8498,8542,8662$ Å) are available in the literature. In
this work, calibrations from Starkenburg et al. (2010, hereafter S10), Carrera
et al. (2013, hereafter C13), and Longeard et al. (2023, hereafter L23) are
selected. All of them are in the form:

${\rm[Fe/H]}=a+b\cdot\rm{M_{V}}+c\cdot\rm{\Sigma Ca}+d\cdot\rm{\Sigma
Ca}^{-1.5}+e\cdot\rm{\Sigma Ca}\cdot\rm{M_{V}},$ (1)

where $a,b,c,d$ are coefficients, $\rm{M_{V}}$ is the absolute V magnitude of
the giant star, $\rm{\Sigma Ca}$ is a linear function of the EWs of the
components of the Ca ii T lines. In the case of the calibration from S10,
$\rm{\Sigma Ca}$ is the sum of the EWs of the second and third components of
the Triplet, $\rm{\Sigma Ca}=\rm{EW}_{2}+\rm{EW}_{3}$. The calibration from
C13 uses $\rm{\Sigma Ca}$ as the weighted sum of the three components,
$\rm{\Sigma
Ca}=0.19\cdot\rm{EW}_{1}+0.47\cdot\rm{EW}_{2}+0.34\cdot\rm{EW}_{3}$. L23 use
the first and the second lines of the triplet, $\rm{\Sigma
Ca}=\rm{EW}_{1}+\rm{EW}_{2}$.

$\rm{M_{V}}$ is derived converting the Gaia EDR3 magnitudes to the Johnson-
Cousin filter following Riello et al. (2021) and adopting a distance modulus
of $19.67\pm 0.14$ mag (see Table 1). The EW is measured using the splot
routine in IRAF, fitting the spectral lines with multiple profiles. The median
and the standard deviation of the multiple measurements have been adopted as
values for the EW and its uncertainty. A Monte Carlo with $10^{6}$
randomisations is performed on the input variables of Equation 1 and on the
coefficients drawing them from a Gaussian distribution. In case of the
coefficients from L23, the randomisations have been done taking into account
their strong correlations (see L23). While the relative uncertainty on the
coefficients from S10 and C13 is set to 10 percent. The [Fe/H] and its
uncertainty are the median and the standard deviation of the draws,
respectively. Metallicities and radial velocities are reported in Table 5. The
mean of the three metallicity estimates is reported as well. We note that the
metallicities from S10 and C13 differ by $\sim 1.35\sigma$, while the ones
from S10 and L23 by only $\sim 0.33\sigma$. One source of discrepancy might be
the different weight in $\rm{\Sigma Ca}$ given to the Ca ii T lines by the
different calibrations. While they have the same weights in S10 and L23, the
first component of the triplet has a less impact in C13 (first : second is
$1:2.47$ and first : third is $1:1.79$). As discussed below, the small
discrepancy in the metallicities does not affect the membership of the two
candidates.

Figure 3 displays the metallicities and radial velocities of our targets and
known Sculptor members (from DART and APOGEE DR17) as a function of their
elliptical distances (left panels); the [Fe/H] vs. RV space and their
histograms (central and right panels). The three calibrations provide [Fe/H]
for both Target 1 and 2 in agreement with the Sculptor’s range. As regards
Target 1, we note that our RV determination places this star at the edge of
the distribution for this system. Its RV is compatible with the value from a
few known members observed by the DART survey and located in the inner region
($\lesssim 6r_{h}$) of the system (see top left panel of Figure 3). The RV of
Target 2 is in agreement with the Sculptor distribution and also with the
previous determination from Westfall et al. (2006).

###  4.3 Comparison with the Besançon model

Figure 3 also shows the Besançon simulation666https://model.obs-besancon.fr
(Robin et al., 2003, 2017) of the Milky Way halo towards Sculptor, with the
same cut on sky position as in the left panel of Figure 1. Nearby foreground
stellar particles are removed selecting a heliocentric distance greater than 5
kpc. Of the initial 230 stellar particles, only 55 inhabit the same
RV$-$[Fe/H] region as Sculptor. Of the 55 stellar particles, only one has a
proper motion compatible with Sculptor members. However, its photometry in the
G band differs more than a dex from the Sculptor stars at the same (BP $-$ RP)
colour. None of the simulated stellar particles can reproduce all the Sculptor
properties simultaneously. We confirm Target 1 and 2 as new Sculptor members.

Figure 3: Distributions of the Sculptor stars. Left panels: Radial velocities
(top) and metallicities (bottom) as a function of the elliptical distance.
Central panel: distribution of Sculptor stars in the [Fe/H] vs. RV space.
Corner plots: histograms of the RV (top) and metallicity (right)
distributions. Targets 1 and 2 are marked with a diamond and a star marker,
respectively. The metallicities from three calibration are denoted by dark
blue (Starkenburg et al., 2010), grey (Carrera et al., 2013), and magenta
(Longeard et al., 2023) markers. RV of Target 2 from Westfall et al. (2006) is
denoted by a dark orange star, which is in agreement with our measurement.
Lime green dots are the cleaned sample from DART, while hot pink squares are
Sculptor members from APOGEE DR17. Stellar particles from the Besançon
simulation are marked with black in the RV vs. [Fe/H] and the histograms
panels. Table 5: Radial velocities, equivalent widths, and metallicities of
the targets. The metallicities from the three calibrations are reported, S10
(Starkenburg et al., 2010), C13 (Carrera et al., 2013) and L23 (Longeard et
al., 2023), and the mean from the three inferences.

Target | RV | EW1 | EW2 | EW3 | [Fe/H] | [Fe/H] | [Fe/H] | $\overline{{\rm[Fe/H]}}$
---|---|---|---|---|---|---|---|---
| (km s${}^{-1})$ | (mÅ) | (mÅ) | (mÅ) | S10 | C13 | L23 |
Target 1 | $148.6\pm 2.1\pm 13.3$ | $940\pm 41$ | $1900\pm 105$ | $1640\pm 31$ | $-1.95\pm 0.33$ | $-2.66\pm 0.37$ | $-2.08\pm 0.11$ | $-2.23\pm 0.38$
Target 2 | $112.4\pm 2.5\pm 13.3$ | $912\pm 46$ | $1999\pm 86$ | $1462\pm 35$ | $-1.93\pm 0.33$ | $-2.59\pm 0.37$ | $-2.03\pm 0.09$ | $-2.18\pm 0.36$

##  5 Discussion

###  5.1 New Sculptor members

The two candidate members have positions, photometry, proper motions, radial
velocities, and metallicities compatible with known stars in Sculptor (see
Figures 1 and 3). These properties cannot be reproduced at the same time by
the Besançon simulations. We conclude that the two targets analysed in this
work are new members of Sculptor. Therefore, Sculptor extends at least out to
a projected elliptical distance of 10 $r_{h}$, which is a physical distance of
$\sim 3$ kpc from its centre.

###  5.2 Metallicity gradients

Figure 3 shows that the metal-rich population (${\rm[Fe/H]}{}\gtrsim-1.5$) is
more centrally concentrated and kinematically colder than the metal-poor stars
(${\rm[Fe/H]}{}\lesssim-1.5$). As a matter of fact, 90 percent of the metal-
rich population lies within $2.0r_{h}$, while the same fraction of the more
metal-poor stars is contained within $3.7r_{h}$. This result is in agreement
with previous work (e.g., Tolstoy et al., 2004; Battaglia et al., 2008; Walker
& Peñarrubia, 2011; Tolstoy et al., 2023), in which the metal-rich and
centrally concentrated population has a smaller velocity dispersion ($\sim
6-7{\rm\,km\ s^{-1}}$) than the more metal-poor and diffuse population ($\sim
10-12{\rm\,km\ s^{-1}}$). This gradient in metallicity can be explained by an
outside-in star formation history or could be the result of an early merger
which dispersed the metal poor component before the subsequent formation of
the metal-rich stars. (e.g., Zhang et al., 2012; Benítez-Llambay et al., 2016;
Revaz & Jablonka, 2018). While the presence of a gradient in Sculptor can be
an effect of mechanisms in play at early times, the presence of stars in the
extreme outskirts might be a signature of recent tidal perturbations, as
discussed in the following section.

###  5.3 Sculptor is affected by tides

At face value, the identification of member stars as far out as 10 half-light
radii from the centre of Sculptor is rather surprising (see also Qi et al.,
2022). For an exponential profile, the surface number density of stars at that
radius is expected to drop by 6 orders of magnitude from its central value
and, therefore, the chance of detecting member stars that far from the centre
should be negligible. Indeed, as shown in the top-left panel of Figure 4, the
stars presented in our study are, at $\sim 100$ arcmin, part of a clear excess
over an exponential profile (solid red line) fit to the inner regions of
Sculptor’s surface number density profile. Because of this excess, the surface
number density at 10 half light radii is roughly 100 times higher than
expected from an exponential profile.

What is the origin of this excess? Is it innate to Sculptor, or is it the
result of internal processes that have driven stars out, such as gravitational
fluctuations driven by the gas inflows and outflows (e.g., El-Badry et al.,
2018), or of external processes, such as early accretion or mergers (Benítez-
Llambay et al., 2016)? Or is it perhaps driven by the effects of Galactic
tides? Although it is difficult to make a definitive case for either scenario
without further data or detailed simulation, we believe that the data favours
a tidal interpretation, as detailed below.

Scl is a satellite of the Milky Way, currently at $\sim 86$ kpc from its
centre (e.g., McConnachie, 2012), and moving out on an orbit with pericentric
and apocentric distance of $\sim 55$ kpc and $\gtrsim 145$ kpc, respectively
(Li et al., 2022; Battaglia et al., 2022; Pace, Erkal & Li, 2022). Despite the
relatively large inferred pericentre, Sculptor shows sign of internal velocity
gradients, which have been interpreted as a consequence of Galactic tides
(see, e.g., Martínez-García, del Pino & Aparicio, 2023). Tides are expected
also to leave their imprint on the surface number density profile ($\Sigma$)
of a dwarf. This imprint is best identified by considering the logarithmic
slope of the profile, $\Gamma=\rm{d}\log\Sigma/\rm{d}\log r$, which is shown,
as a function of radius, in the bottom left panel of Figure 4. For an
exponential profile, $\Gamma$ decreases linearly with radius, as shown by the
solid red line, with a slope equal to the inverse of the exponential scale
length.

As a result of a tidal perturbation, some stars may gain energy and move to
orbits with larger apocentric radii, creating an outer "excess" of stars that
resemble that of Sculptor. A signature of this perturbation is a sudden
departure of $\Gamma$ from a linear dependence on radius; i.e., a "kink" in
the $\Gamma$ profile like the one clearly seen for Sculptor at
r${}_{\rm{ell}}\sim 25$ arcmin from the centre. Beyond the kink, the tidal
excess should approach a power-law of slope $\Gamma=-4$, as expected if stars
have been moved to very loosely bound orbits, in effect populating phase space
all the way to zero energy, E $\sim 0$ (White et al., 1987; Jaffe, 1987;
Peñarrubia et al., 2009).

This behaviour is illustrated in Figure 4 by the dashed and dot-dashed black
curves, which are taken from the tidally-perturbed models of Peñarrubia,
Navarro & McConnachie (2008, see the top-left panel of their Figure 4). After
scaling to the same half-light radius, the tidally perturbed (dot-dashed)
profile fits remarkably well the outer excess seen in Sculptor, including the
predicted R-4 behaviour in the outer regions, even though this simulation was
not intended in any way to reproduce Sculptor.

In the simulations, the $R^{-4}$ region ends at a “break” radius where the
$\Sigma$ profile flattens and $\Gamma$ increases sharply (this occurs at
r${}_{\rm{ell}}\sim 80$ arcmin in the simulated profile). The location of this
"break" corresponds roughly to where the average motion of stars transitions
from equilibrium ($\langle v_{r}\rangle\approx 0$) to becoming dominated by
outward motion ($\langle v_{r}\rangle>0$), and corresponds to the radius where
the local crossing time equals the time elapsed since pericenter (see Aguilar
& White, 1986; Navarro, 1990; Peñarrubia et al., 2009). Outside the break,
stars are still moving to their new apocentres and have not yet reached
equilibrium.

In the case of Sculptor, assuming that the system last passed pericentre $\sim
0.41$ Gyr ago (Battaglia et al., 2022), the break radius is expected777The
relation from Peñarrubia et al. (2009) is used,
r${}_{b}~{}=~{}C~{}\cdot~{}\sigma_{\rm{v}}~{}\cdot~{}t_{\rm{peri}}$, where C
is a coefficient ($=0.55$), $\sigma_{\rm{v}}$ is the velocity dispersion of
the system ($9.2{\rm\,km\ s^{-1}}$, see Table 1), and $t_{\rm{peri}}$ is the
time since the last pericentric passage. For the latter, the value derived by
Battaglia et al. (2022) using the MW potential perturbed by the Large
Magellanic Cloud is adopted. to be at $\sim 2.1$ kpc ($\sim 85$ arcmin or
$\sim 6.9$ rh), very close to the two outermost points of the $\Sigma$
profile888It is quite difficult to estimate robustly the logarithmic slope of
the profile at that radius given current data. A finer binning of the J23 data
would show a hint of a plateau in the $\Sigma$ profile around the inferred
break radius. However, a finer grid leads to a much higher noise in the
$\Gamma$ profile. We defer a more detailed analysis to a future contribution..
This is a robust prediction of the tidal interpretation that could, in
principle, be verified in future work. The "kink" in the $\Gamma$ profile and
the $\Gamma=-4$ are not expected in models where the excess stars in the
outskirts are a result of mergers or potential fluctuations (e.g., Peñarrubia
et al., 2009; Benítez-Llambay et al., 2016), and clearly favour a tidal
interpretation.

Could the outer excess be instead just an innate feature of Sculptor’s density
profile rather than a result of a tidal perturbation? If the outer excess is
actually tidal in origin, then it should be absent in dwarf spheroidals that
have not been affected by tides. This is the case of Fornax, whose large
Galactocentric distance, together with its large orbital pericentre, indicate
that tides have played at best a minor role in perturbing its structure (e.g.,
Borukhovetskaya et al., 2022; Yang et al., 2022, and references therein). As
shown in the right-hand panels of Figure 4, Fornax shows a profile with little
evidence for any excess over the whole radial range plotted here, which spans
more than 5 decades in surface number density.

We conclude that the outer excess in the density profile of Sculptor has been
driven by tides. In addition, the orbital motion of Target 1 and 2 can be
perturbed by the tidal effects. Further work designed to confirm that the
tidal signatures described above are consistent with Sculptor’s relatively
large pericentre, and search for corroborating evidence for this conclusion,
such as the presence of a "break radius" in its outer profile, or even the
presence of faint tidal tails even further out, is clearly needed.

Figure 4: Surface number density distribution $\Sigma$ and its spatial
logarithmic derivative $\Gamma$. Left panels: Sculptor. Right panels: Fornax.
Top and bottom panels: Sculptor candidate members from J23 and MV20 are marked
with blue and empty red circles, respectively. Fornax candidate members
(magenta circles) are from J23. P08 models are marked with a dashed line (no
tide) and with a dash-dotted line (after first apocentric passage model),
which are reported only for Sculptor. Red line is the exponential fit to J23
data. Grey arrow marks the position of the two outermost members analysed in
this work. Pink and orange ticks denotes the core and the half-light radius of
the systems. Central panels: difference between J23 data and the exponential
fit. For Sculptor, the exponential fit is able to reproduce the inner surface
number density distribution of the MV20 data as well.

##  6 Conclusions

The two outermost candidate members of Sculptor dwarf galaxy have been
selected from the Bayesian algorithm of Jensen et al. (2023, in prep.) and
observed with GMOS spectrograph at Gemini South. All the candidate members
(J23, MV20) and known Scl members from DART and APOGEE DR17 have been
gathered. We find that:

  1. 1.

the two candidates are actually new members given their radial velocity and
metallicity. Therefore, Scl is now known to extend out to a projected
elliptical distance of $\sim 10$ rh, or $\sim 3$ kpc, from its centre (see
Figures 1 and 3);

  2. 2.

as noted in earlier work, the radial velocity and metallicity distributions of
Sculptor indicate that the most metal-poor stars are more dispersed spatially
and kinematically than the more metal-rich population (see Figure 3). This
could be explained by an outside-in star formation scenario or an early
merger;

  3. 3.

the surface number density distribution and its logarithmic derivative (see
Figure 4) suggest that Sculptor has been influenced by tidal effects, starting
from a distance of $\sim 25$ arcmin, or $\sim 2$ rh, or $\sim 620$ pc from its
centre;

  4. 4.

tides are likely responsible for the presence of Target 1 and 2 in the extreme
outskirts of Sculptor.

Recent work on the outskirts of ultra faint and dwarf galaxies (Chiti et al.,
2021; Filion & Wyse, 2021; Longeard et al., 2022; Yang et al., 2022; Waller et
al., 2023; Sestito et al., 2023b; Longeard et al., 2023; Jensen et al., 2023,
in prep., this work) has uncovered a new population of stars that would have
been undiscoverable without the powerful synergy between the exquisite Gaia
data, efficient algorithms to remove foreground stars, and large aperture
telescopes.

## Acknowledgements

We acknowledge and respect the lkwŋn peoples on whose traditional territory
the University of Victoria stands and the Songhees, Esquimalt and
$\underaccent{\bar}{\rm W}$SÁNEĆ peoples whose historical relationships with
the land continue to this day.

We thank Nicolas Longeard for sharing the Monte Carlo chain drawings from
their metallicity calibration.

We want to thank the supporter astronomers Venu Kalari and Juan Font-Serra for
their support during Phase II.

FS thanks the Dr. Margaret "Marmie" Perkins Hess postdoctoral fellowship for
funding his work at the University of Victoria. KAV thanks the National
Sciences and Engineering Research Council of Canada for funding through the
Discovery Grants and CREATE programs.

This work is based on observations obtained with Gemini South/GMOS, as part of
the Gemini Fast Turnaround Program GS-2022B-FT-205. Based on observations
obtained at the international Gemini Observatory, a program of NSF’s NOIRLab,
which is managed by the Association of Universities for Research in Astronomy
(AURA) under a cooperative agreement with the National Science Foundation. On
behalf of the Gemini Observatory partnership: the National Science Foundation
(United States), National Research Council (Canada), Agencia Nacional de
Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e
Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e
Comunicações (Brazil), and Korea Astronomy and Space Science Institute
(Republic of Korea).

This work has made use of data from the European Space Agency (ESA) mission
Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P.
Sloan Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions. SDSS-IV acknowledges support and resources from
the Center for High-Performance Computing at the University of Utah. The SDSS
web site is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the
Participating Institutions of the SDSS Collaboration including the Brazilian
Participation Group, the Carnegie Institution for Science, Carnegie Mellon
University, the Chilean Participation Group, the French Participation Group,
Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de
Canarias, The Johns Hopkins University, Kavli Institute for the Physics and
Mathematics of the Universe (IPMU) / University of Tokyo, the Korean
Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut
für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA
Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-
Institut für Extraterrestrische Physik (MPE), National Astronomical
Observatories of China, New Mexico State University, New York University,
University of Notre Dame, Observatário Nacional / MCTI, The Ohio State
University, Pennsylvania State University, Shanghai Astronomical Observatory,
United Kingdom Participation Group, Universidad Nacional Autónoma de México,
University of Arizona, University of Colorado Boulder, University of Oxford,
University of Portsmouth, University of Utah, University of Virginia,
University of Washington, University of Wisconsin, Vanderbilt University, and
Yale University.

This research has made use of the SIMBAD database, operated at CDS,
Strasbourg, France (Wenger et al., 2000). This work made extensive use of
TOPCAT (Taylor, 2005).

## Data Availability

GMOS spectra will be available at the Gemini Archive web page
https://archive.gemini.edu/searchform after the proprietary time. All data are
incorporated into the article.

## References

  * Abdurro’uf et al. (2022) Abdurro’uf et al., 2022, ApJS, 259, 35
  * Aguilar & White (1986) Aguilar L. A., White S. D. M., 1986, ApJ, 307, 97
  * Andrae et al. (2018) Andrae R. et al., 2018, A&A, 616, A8
  * Battaglia et al. (2008) Battaglia G., Helmi A., Tolstoy E., Irwin M., Hill V., Jablonka P., 2008, ApJ, 681, L13
  * Battaglia et al. (2022) Battaglia G., Taibi S., Thomas G. F., Fritz T. K., 2022, A&A, 657, A54
  * Benítez-Llambay et al. (2016) Benítez-Llambay A., Navarro J. F., Abadi M. G., Gottlöber S., Yepes G., Hoffman Y., Steinmetz M., 2016, MNRAS, 456, 1185
  * Borukhovetskaya et al. (2022) Borukhovetskaya A., Errani R., Navarro J. F., Fattahi A., Santos-Santos I., 2022, MNRAS, 509, 5330
  * Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
  * Carrera et al. (2013) Carrera R., Pancino E., Gallart C., del Pino A., 2013, MNRAS, 434, 1681
  * Chambers et al. (2016) Chambers K. C. et al., 2016, The pan-starrs1 surveys
  * Chiti et al. (2021) Chiti A. et al., 2021, Nature Astronomy, 5, 392
  * Deason, Wetzel & Garrison-Kimmel (2014) Deason A., Wetzel A., Garrison-Kimmel S., 2014, ApJ, 794, 115
  * Deason et al. (2022) Deason A. J., Bose S., Fattahi A., Amorisco N. C., Hellwing W., Frenk C. S., 2022, MNRAS, 511, 4044
  * Dey et al. (2019) Dey A. et al., 2019, AJ, 157, 168
  * Drlica-Wagner et al. (2021) Drlica-Wagner A. et al., 2021, ApJS, 256, 2
  * El-Badry et al. (2018) El-Badry K. et al., 2018, MNRAS, 477, 1536
  * Evans et al. (2018) Evans D. W. et al., 2018, A&A, 616, A4
  * Fattahi et al. (2018) Fattahi A., Navarro J. F., Frenk C. S., Oman K. A., Sawala T., Schaller M., 2018, MNRAS, 476, 3816
  * Fenner et al. (2006) Fenner Y., Gibson B. K., Gallino R., Lugaro M., 2006, ApJ, 646, 184
  * Filion & Wyse (2021) Filion C., Wyse R. F. G., 2021, ApJ, 923, 218
  * Frebel, Kirby & Simon (2010) Frebel A., Kirby E. N., Simon J. D., 2010, Nature, 464, 72
  * Frenk et al. (1988) Frenk C. S., White S. D. M., Davis M., Efstathiou G., 1988, ApJ, 327, 507
  * Gaia Collaboration et al. (2021) Gaia Collaboration et al., 2021, A&A, 649, A1
  * Genel et al. (2010) Genel S., Bouché N., Naab T., Sternberg A., Genzel R., 2010, ApJ, 719, 229
  * Gimeno et al. (2016) Gimeno G. et al., 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, Evans C. J., Simard L., Takami H., eds., p. 99082S
  * González Hernández & Bonifacio (2009) González Hernández J. I., Bonifacio P., 2009, A&A, 497, 497
  * Grebel, Gallagher & Harbeck (2003) Grebel E. K., Gallagher, John S. I., Harbeck D., 2003, AJ, 125, 1926
  * Gustafsson et al. (2008) Gustafsson B., Edvardsson B., Eriksson K., Jørgensen U. G., Nordlund Å., Plez B., 2008, A&A, 486, 951
  * Hartwig et al. (2023) Hartwig T., Ishigaki M. N., Kobayashi C., Tominaga N., Nomoto K., 2023, ApJ, 946, 20
  * Higgs et al. (2021) Higgs C. R., McConnachie A. W., Annau N., Irwin M., Battaglia G., Côté P., Lewis G. F., Venn K., 2021, MNRAS, 503, 176
  * Hill et al. (2019) Hill V. et al., 2019, A&A, 626, A15
  * Hook et al. (2004) Hook I. M., Jørgensen I., Allington-Smith J. R., Davies R. L., Metcalfe N., Murowinski R. G., Crampton D., 2004, PASP, 116, 425
  * Ishigaki et al. (2018) Ishigaki M. N., Tominaga N., Kobayashi C., Nomoto K., 2018, ApJ, 857, 46
  * Jablonka et al. (2015) Jablonka P. et al., 2015, A&A, 583, A67
  * Jaffe (1987) Jaffe W., 1987, in Structure and Dynamics of Elliptical Galaxies, de Zeeuw P. T., ed., Vol. 127, p. 511
  * Jensen et al. (2023, in prep.) Jensen et al., 2023, in prep.
  * Kazantzidis et al. (2011) Kazantzidis S., Łokas E. L., Callegari S., Mayer L., Moustakas L. A., 2011, ApJ, 726, 98
  * Kielty et al. (2021) Kielty C. L. et al., 2021, MNRAS, 506, 1438
  * Lardo et al. (2021) Lardo C. et al., 2021, MNRAS, 508, 3068
  * Li et al. (2022) Li T. S. et al., 2022, ApJ, 928, 30
  * Longeard et al. (2022) Longeard N. et al., 2022, MNRAS, 516, 2348
  * Longeard et al. (2023) —, 2023, arXiv e-prints, arXiv:2304.13046
  * Marigo et al. (2008) Marigo P., Girardi L., Bressan A., Groenewegen M. A. T., Silva L., Granato G. L., 2008, A&A, 482, 883
  * Martin et al. (2009) Martin N. F. et al., 2009, ApJ, 705, 758
  * Martínez-García, del Pino & Aparicio (2023) Martínez-García A. M., del Pino A., Aparicio A., 2023, MNRAS, 518, 3083
  * McConnachie (2012) McConnachie A. W., 2012, AJ, 144, 4
  * McConnachie et al. (2009) McConnachie A. W. et al., 2009, Nature, 461, 66
  * McConnachie & Venn (2020a) McConnachie A. W., Venn K. A., 2020a, AJ, 160, 124
  * McConnachie & Venn (2020b) —, 2020b, Research Notes of the American Astronomical Society, 4, 229
  * Moster, Naab & White (2013) Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121
  * Mucciarelli, Bellazzini & Massari (2021) Mucciarelli A., Bellazzini M., Massari D., 2021, A&A, 653, A90
  * Navarro (1990) Navarro J. F., 1990, MNRAS, 242, 311
  * Navarro, Frenk & White (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Pace, Erkal & Li (2022) Pace A. B., Erkal D., Li T. S., 2022, ApJ, 940, 136
  * Peñarrubia, Navarro & McConnachie (2008) Peñarrubia J., Navarro J. F., McConnachie A. W., 2008, ApJ, 673, 226
  * Peñarrubia et al. (2009) Peñarrubia J., Navarro J. F., McConnachie A. W., Martin N. F., 2009, ApJ, 698, 222
  * Placco et al. (2021) Placco V. M., Sneden C., Roederer I. U., Lawler J. E., Den Hartog E. A., Hejazi N., Maas Z., Bernath P., 2021, Research Notes of the American Astronomical Society, 5, 92
  * Plez (2012) Plez B., 2012, Turbospectrum: Code for spectral synthesis
  * Qi et al. (2022) Qi Y., Zivick P., Pace A. B., Riley A. H., Strigari L. E., 2022, MNRAS, 512, 5601
  * Revaz & Jablonka (2018) Revaz Y., Jablonka P., 2018, A&A, 616, A96
  * Riello et al. (2021) Riello M. et al., 2021, A&A, 649, A3
  * Robin et al. (2017) Robin A. C., Bienaymé O., Fernández-Trincado J. G., Reylé C., 2017, A&A, 605, A1
  * Robin et al. (2003) Robin A. C., Reylé C., Derrière S., Picaud S., 2003, A&A, 409, 523
  * Schlafly & Finkbeiner (2011) Schlafly E. F., Finkbeiner D. P., 2011, ApJ, 737, 103
  * Schultz & Wiemer (1975) Schultz G. V., Wiemer W., 1975, A&A, 43, 133
  * Sestito et al. (2019) Sestito F. et al., 2019, MNRAS, 484, 2166
  * Sestito et al. (2023a) —, 2023a, MNRAS, 518, 4557
  * Sestito et al. (2023b) —, 2023b, arXiv e-prints, arXiv:2301.13214
  * Shapley (1938) Shapley H., 1938, Harvard College Observatory Bulletin, 908, 1
  * Simon (2019) Simon J. D., 2019, ARA&A, 57, 375
  * Simon et al. (2015) Simon J. D., Jacobson H. R., Frebel A., Thompson I. B., Adams J. J., Shectman S. A., 2015, ApJ, 802, 93
  * Skúladóttir et al. (2021) Skúladóttir Á. et al., 2021, ApJ, 915, L30
  * Sneden (1973) Sneden C. A., 1973, PhD thesis, THE UNIVERSITY OF TEXAS AT AUSTIN.
  * Starkenburg et al. (2013) Starkenburg E. et al., 2013, A&A, 549, A88
  * Starkenburg et al. (2010) —, 2010, A&A, 513, A34
  * Starkenburg et al. (2017) —, 2017, MNRAS, 471, 2587
  * Tafelmeyer et al. (2010) Tafelmeyer M. et al., 2010, A&A, 524, A58
  * Taylor (2005) Taylor M. B., 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, Shopbell P., Britton M., Ebert R., eds., p. 29
  * The Dark Energy Survey Collaboration (2005) The Dark Energy Survey Collaboration, 2005, arXiv e-prints, astro
  * Tody (1986) Tody D., 1986, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 627, Society of Photo-Optical Instrumentation Engineers (SPIE), Crawford D. L., ed., p. 733
  * Tody (1993) —, 1993, in Astronomical Society of the Pacific Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, Hanisch R. J., Brissenden R. J. V., Barnes J., eds., p. 173
  * Tolstoy et al. (2006) Tolstoy E. et al., 2006, The Messenger, 123, 33
  * Tolstoy, Hill & Tosi (2009) Tolstoy E., Hill V., Tosi M., 2009, ARA&A, 47, 371
  * Tolstoy et al. (2004) Tolstoy E. et al., 2004, ApJ, 617, L119
  * Tolstoy et al. (2023) —, 2023, arXiv e-prints, arXiv:2304.11980
  * Venn et al. (2004) Venn K. A., Irwin M., Shetrone M. D., Tout C. A., Hill V., Tolstoy E., 2004, AJ, 128, 1177
  * Walker & Peñarrubia (2011) Walker M. G., Peñarrubia J., 2011, ApJ, 742, 20
  * Waller et al. (2023) Waller F. et al., 2023, MNRAS, 519, 1349
  * Wenger et al. (2000) Wenger M. et al., 2000, A&AS, 143, 9
  * Westfall et al. (2006) Westfall K. B., Majewski S. R., Ostheimer J. C., Frinchaboy P. M., Kunkel W. E., Patterson R. J., Link R., 2006, AJ, 131, 375
  * Wheeler et al. (2019) Wheeler C. et al., 2019, MNRAS, 490, 4447
  * White et al. (1987) White S. D. M., Davis M., Efstathiou G., Frenk C. S., 1987, Nature, 330, 451
  * White & Rees (1978) White S. D. M., Rees M. J., 1978, MNRAS, 183, 341
  * Yang et al. (2022) Yang Y., Hammer F., Jiao Y., Pawlowski M. S., 2022, MNRAS, 512, 4171
  * York et al. (2000) York D. G. et al., 2000, AJ, 120, 1579
  * Zhang et al. (2012) Zhang H.-X., Hunter D. A., Elmegreen B. G., Gao Y., Schruba A., 2012, AJ, 143, 47


acronymindexonlyfirsttrue [acronym]long-short acronymnohyperfirsttrue

#  Autonomous Forest Inventory with Legged Robots:
System Design and Field Deployment

Matias Mattamala1, Nived Chebrolu1, Benoit Casseau1, Leonard Freißmuth1,3,
Jonas Frey2,4, Turcan Tuna2, Marco Hutter2, Maurice Fallon1 1The authors are
with the Oxford Robotics Institute at the University of Oxford, UK. {matias,
nived, benoit<EMAIL_ADDRESS>authors are with the Robotic
Systems Lab, ETH Zurich. {jonfrey, tutuna<EMAIL_ADDRESS>author is
also with the Technical University of Munich, Germany<EMAIL_ADDRESS>4The author is also with the Autonomous Learning Group, Max Planck Institute
for Intelligent Systems.

###### Abstract

We present a solution for autonomous forest inventory with a legged robotic
platform. Compared to their wheeled and aerial counterparts, legged platforms
offer an attractive balance of endurance and low soil impact for forest
applications. In this paper, we present the complete system architecture of
our forest inventory solution which includes state estimation, navigation,
mission planning, and real-time tree segmentation and trait estimation. We
present preliminary results for three campaigns in forests in Finland and the
UK and summarize the main outcomes, lessons, and challenges. Our UK experiment
at the Forest of Dean with the ANYmal D legged platform, achieved an
autonomous survey of a $0.96\text{\,}\mathrm{ha}\text{/}$ plot in
$20\text{\,}\mathrm{min}\text{/}$, identifying over 100 trees with typical DBH
accuracy of $2\text{\,}\mathrm{cm}\text{/}$.

##  I Introduction

Systematically measuring tree attributes such as the tree height, dbh, and
canopy volume, is important for modern forestry [1]. This process, known as
forest inventory, has traditionally involved manual tree measurement using
measuring tapes and calipers and also handheld laser rangefinders [2].
Nowadays, remote sensing techniques such as tls [3] and aerial photograph [4]
have massively improved the accuracy and scale at which forest inventories are
built [5].

Robotics has the potential to scale this process even further. Technologies
such as slam have been used for handheld forest mapping [6, 7]. Fixed-wing
aircrafts [8] and hexacopters [9] have enabled autonomous data collection
above the tree canopy where there is no risk of collision. Meanwhile small-
scale drones [10, 11] and wheeled platforms [12] have been recently used for
under-canopy forest inventory missions. Legged robots potentially provide
significant advantages to aerial and wheeled platforms in terms of mobility,
endurance, and lower soil impact [13]. However, robots have not been used for
forestry tasks apart from tree harvesting [14].

In this work, we study the application of legged platforms for autonomous
forest inventory tasks. Our solution combines the advanced mobility skills
developed by the robotics research community [15], a robust LiDAR-based
mapping system with a multi-level autonomous mission planning architecture.
Our solution leverages these components to provide a consistent representation
of the environment, while also enabling online processing and concurrent
estimation of the forest inventory.

This paper presents ongoing work, reporting autonomous deployments executed in
May 2023, October 2023, and February 2024 using our solution. Our experiments
have been carried out in various forests in Finland and the UK, as shown in
Fig. 1, and reflect the evolving design decisions we have made after each
deployment. We discuss the findings and lessons learned, as well as the main
challenges towards future autonomous forest inventory missions in the context
of the DigiForest EU Horizon project111https://www.digiforest.eu.

Figure 1: We present an autonomous forest inventory solution with legged
platforms. Our system aims to autonomously drive a legged platform in an
unknown forest plot, while recording environmental data used for tree
segmentation and trait estimation.

##  II Related Work

###  II-A Robotics-enabled forestry

Robotics has the potential to scale up forest mensuration through technologies
and platforms that can automate intensive manual labor. slam can accelerate
ground-based forest mensuration by integrating LiDAR scans from a mobile
sensing payload. The produced maps can be post-processed to extract the
desired attributes, such as the tree positions and dbh [6, 16, 7]. The
incremental nature of slam and the data acquisition process has motivated the
concurrent construction of maps and inventories. Chen et al. [10] and Proudman
et al. [17], demonstrated real-time tree detection and dbh estimation using
small-scale quadcopters and handheld devices, respectively. We further explore
this idea, developing an autonomy system that builds upon similar online
forest inventory approaches [18].

Robotic platforms can support surveying tasks by freeing humans from
collecting the data themselves [19]. In the context of aerial surveying,
autonomous fixed-wing aircrafts [8], small-scale helicopters [20] and more
recently hexacopter platforms [9] have been used for over-canopy data
collection. As these platforms fly over the trees, most of the robotic
challenges concern achieving maximum area coverage given energy [21] and time
constraints [22].

Under-canopy data collection has shown more challenging, as it requires robots
to operate in cluttered, unstructured spaces. For aerial platforms, this
requires the use of smaller drones compared to the ones that can fly over
canopy—with smaller payloads and battery size [23]. Teleoperated small-scale
drones have been used for more than ten years for surveying [24, 25] though
autonomous operation has been recently demonstrated for forest inventory [10,
11] and DNA sampling [26].

On the contrary, ground platforms are able to provide under-canopy
measurements while carrying heavier sensing payloads, enabling longer surveys.
Wheeled and tracked robot platforms have been used for forest data collection
[6, 12, 27] as well as the harvesting process [28, 14, 29]. However, their
larger footprint and mass contribute to soil compaction damage due to
trampling [30, 31]. In this work we aim to evaluate the performance of smaller
platforms—up to $50\text{\,}\mathrm{kg}\text{/}$.

###  II-B Deployment of legged platforms in natural environments

The recent advances in hardware and control of legged robots makes them an
compelling alternative for ground-level surveying [15].

Legged robots have been used for teleoperated soil sampling [32] as well as
wildlife recording [33]. Autonomous operation in natural environments has also
been demonstrated for search and rescue operations [34], load carrying tasks
[35, 36], measurement of aeolian processes [37, 38], and planetary analog
exploration [39].

Forestry applications have been less well unexplored, though their potential
advantages (heavier payload carrying, lower soil impact) have been discussed
in the past [13]. In this work, we aim to provide experimental demonstration
of their capabilities for forest inventory applications.

##  III Autonomous Forest Inventory with Legged Robots

Figure 2: System overview. The Autonomy system executes the mission from a
survey reference given by a human operator. State estimation provides a
consistent scene representation, as well as dense clouds which are used for
tree traits estimation. Forest analysis segments and estimates tree traits
from point cloud data. The main output of the system is a forest inventory
database with the main attributes of the surveyed forest.

Fig. 2 shows the system architecture of our solution, involving state
estimation, legged autonomy, and forest analysis. The specific tasks and
constituent modules are explained in the following sections.

###  III-A State Estimation

Our state estimation approach involves four modules: (1) LiDAR-Inertial
Odometry, (2) LiDAR SLAM, (3) dense mapping, and (4) local terrain mapping,
which provide representations for forest analysis as well as navigation.

LiDAR-Inertial Odometry: For odometry we use a LiDAR-inertial factor graph-
based odometry solution, VILENS [40], which provides high-frequency pose and
velocity estimates ($\sim$$200\text{\,}\mathrm{Hz}\text{/}$) in a fixed
inertial frame. The system performs sensor fusion of relative pose
displacements from LiDAR using the icp algorithm [41] with preintegrated IMU
measurements [42].

LiDAR SLAM: We use a pose graph LiDAR SLAM system [43], which incrementally
builds a lightweight but consistent representation of the environment. Each
graph node stores a LiDAR scan with the edge factors corresponding to
incremental odometry estimates from LiDAR-Inertial odometry. We detect loop
closures by checking the most recent scan against those linked to pose graph
nodes within a radius of $1015\text{\,}\mathrm{m}\text{/}$ from the robot. The
candidates are geometrically verified using icp, and the successful candidates
are added as additional factors in the pose graph. The graph is incrementally
optimized using the isam2 algorithm [44].

Dense Mapping: As the SLAM system only stores sparse LiDAR scans within the
pose graph, this map is not dense enough for further forestry analysis. We
implemented a local sub-mapping module to accumulate all the incoming LiDAR
scans and generate dense scans in the sensor frame, which we call _data
payloads_. The data payloads are gravity-aligned and published at a slower
rate ($\sim$$0.1\text{\,}\mathrm{Hz}\text{/}$), accumulating clouds for about
every $20\text{\,}\mathrm{m}\text{/}$ traveled. We process these data payload
to carry out online tree traits analysis (see Sec. III-C).

Local Terrain Mapping: In addition to the large-scale mapping modules, we also
carry out local terrain mapping [45] and use its output in the local planner.
We used a 2.5D grid-based representation with a resolution of
$4\text{\,}\mathrm{cm}\text{/}$.

###  III-B Autonomy System

Our autonomy system is multiple levels: spanning from the Survey Interface
(Level 3) down to the low-level locomotion control (Level 0).

Survey Interface (Level 3): The user interface represents the highest level of
planning of our system. The human operator directly interacts to deploy the
robot. We implemented it as a graphical user interface (GUI) on RViz using
interactive markers. Given the initial pose of the robot, the GUI allows the
operator to define a survey area and to initialize a survey plan, given by a
boustrophedon decomposition, also known as lawn mower pattern [46].

The survey plan is concretely defined as a sequence of equally-spaced 6 dof
waypoints, which enforce some minimum constraints to ensure that the slam
system is able to find loop closures during operation and can keep the
trajectory and the map consistent. However, this autonomy level makes no
assumptions about the traversability of the surveyed area and the reachability
of the waypoints, which are delegated to the lower levels of the autonomy
system.

The survey interface also provides mechanisms to interrupt the mission at any
time, which triggers a safe stopping behavior through the lower levels. It
also enables the operator to modify and resume the plan from any specified
goal.

Figure 3: Survey Interface: GUI for the operator, implemented on RViz using
interactive markers. The operator defines the area to be surveyed, while the
survey plan is automatically proposed by the Mission Planner.

Mission Planner (Level 2): The mission planner handles the execution of the
survey plan by interacting with the local planner. In particular, it schedules
the goals sent to the local planner and also processes any feedback obtained
by it. For example, if the local planner reports a goal as unreachable, the
mission planner chooses a new goal based on the survey plan’s state.

Local Planning (Level 1): The local planning level aims to reach the proposed
goals of the plan and checks its feasibility given the state of the local
environment encoded in the terrain map. A traversability estimation system
first processes the local terrain map and extracts geometric features to score
the terrain. The traversability score $s_{\text{trav}}$ is specifically
designed to enable navigation in the forest despite branches and twigs, though
visual traversability approaches could also be used [47].

The traversability representation is used by a reactive local planner [48],
which combines different vector fields for navigation. We extended the
approach to compute a gdf on a continuous cost-to-go $c_{\text{to-go}}$
obtained by the heuristic:

$c_{\text{to-
go}}=w_{\text{trav}}\,(1.0-s_{\text{trav}})+w_{\text{unkn}}\,s_{\text{unkn}}$
(1)

where $s_{\text{unkn}}$ is a fixed score assigned to empty cells in the
terrain map (unknown), and $w_{\text{trav}}$, $w_{\text{unkn}}$ are manually
defined weights to leverage the contribution of the traversability score of
known cells, as well as unknown ones.

The local planner measures the progress when trying to reach a waypoint, in
order to detect potential situations were it is unfeasible (e.g., if the
proposed waypoint is in the middle of bushes). If a goal is determined
unreachable, the local planner reports this situation to the mission planner,
triggering a replanning of the next goal.

Locomotion Controller (Level 0): The output of the local planner is a 3 dof
velocity command—linear $(x,y)$, angular (yaw)—, which is executed by a low-
level reinforcement learning-based locomotion controller [15]. The controller
specifies a nominal walking wait for flat ground, but it relies on the local
terrain map to negotiate uneven, rigid terrain. We used a pre-trained policy
and did not implement any specific modifications to the architecture or
training environment to operate in the forest.

###  III-C Online Forest Analysis

Our system integrates a state estimation-driven forest inventory pipeline,
presented in Freißmuth et al. [18]. The approach enables near real-time tree
segmentation and reconstruction during the mission execution. The main steps
are explained as follows, and please refer to the original work for more
technical details of the implementation.

Terrain Modeling: The input to the online forest analysis step are the data
payloads at $0.1\text{\,}\mathrm{Hz}\text{/}$. For each payload, we extract a
local terrain model using the _cloth simulation filtering_ method by Zhang et
al. [49]. This enables us to segment the points belonging to the terrain from
the potential tree points.

Tree Segmentation: To segment trees we use a method inspired by Cabo et al.
[50]. We first normalize the payload cloud z-values using the terrain model,
and only keep the points in a slice where we expect no tree foliage. We then
use a Voronoi-based clustering to obtain tree cluster candidates. We refine
the points associated to each tree by fitting a cylinder model and only
keeping the points within a distance from the cylinder’s surface. The filtered
tree points per tree candidates are finally de-normalized by height to obtain
tree candidate clouds in the original reference frame.

Spatio-temporal Aggregation: Our system is designed to incrementally build a
forest inventory, effectively improving the tree estimates and terrain model
as the robot autonomously collects more data. We achieve this by exploiting
the pose graph representation of the LiDAR-SLAM system, by attaching local
terrain models and tree instances clouds to it, and updating them with new
data after adding new poses or loop closure updates. Specifically, the global
terrain model is updated by locally averaging the points around the latest
terrain map using a weighted average. The tree instances clouds are managed
accordingly to decide if new instances correspond to an existing tree (due to
proximity) or a novel instance. We additionally check coverage conditions to
ensure that high points are captured as well as even coverage from different
viewpoints.

Tree Reconstruction and Trait Estimation: Lastly, we estimate the traits
required by the forest inventory. Particularly, we first aim to _reconstruct_
the tree by fitting circles along the stem, and then reconstructing the trunk
as a collection of oblique cone frustums. The reconstructions are then used to
recover geometric traits such as the dbh and tree height. The traits are
reported to the operator and also exported as a _marteloscope_ —a 2D
representation of the forest.

##  IV Experiments

Mission | Robot setup | Metrics |
---|---|---|---
Campaign |  Date |  ID |  Area [$m^{2}$] |  Goals |  Hardware |  Software |  Mission time [$s$] |  Distance traveled [$m$] |  Area covered [ha] |  Interv. [$\\#$] |  MDBI [$m$] |  MTBI [$m$]
Evo, Finland |  2023-05-03 |  M1 |  $40\times 25$ |  Autonomy system |  ANYmal C, Velodyne VLP-16 |  CompSLAM odometry |  575.6 |  270.3 |  0.33 |  2 |  84.8 |  176.1
|  2023-05-03 |  M2 |  $40\times 25$ |  |  |  |  432.4 |  233.6 |  0.32 |  0 |  233.6 |  432.4
|  2023-05-04 |  M3 |  $40\times 35$ |  |  |  |  816.8⋆ |  301.1 |  0.37 |  7 |  30.5 |  72.2
|  2023-05-04 |  M4 |  $35\times 35$ |  |  |  |  988.4 |  336.6 |  0.37 |  7 |  39.4 |  114.1
|  2023-05-05 |  M5 |  $70\times 25$ |  |  |  |  1275.5 |  609.7 |  0.64 |  10 |  51.5 |  99.5
Wytham Woods, UK |  2023-10-06 |  M6 |  $20\times 20$ |  Autonomy system, State estimation |  ANYmal C, Velodyne VLP-16 |  VILENS odometry, VILENS-SLAM |  436.9 |  215.0 |  0.29 |  2 |  69.9 |  136.4
Forest of Dean, UK |  2023-02-19 |  M7 |  $125\times 30$ |  State estimation, Forest inventory |  ANYmal D, Hesai QT64 |  VILENS odometry, VILENS-SLAM |  1283.5 |  665.5 |  0.96 |  8 |  65.9 |  127.7

  * $\star$

This mission was manually interrupted.

TABLE I: Missions summary. We report the specifications of each mission
executed across three campaigns in Finland and the UK, along with the robot
setup used, as well as the main metrics reported on each mission, such as the
Mean Distance Between Interventions (MDBI) and mtbi. Figure 4: Illustrative
examples of the robot deployments during the different campaigns. (a)
Coniferous forest in Evo, Finland. (b) Mixed forest in Wytham Woods, UK. (c)
Oak forest in the Forest of Dean, UK.

###  IV-A Campaign description

We executed three field campaigns in different forests in Finland and the UK,
which are summarized in Tab. I. Five missions (M1-M5) were executed in conifer
forests in Evo, Finland in May 2023. The other two were ran in the UK, in
mixed (M6, at Wytham Woods, in October 2023) and oak forests (M7, at the
Forest of Dean, in February 2024). While the two former had the objective of
testing the autonomy system, the latter also focused on the online forest
analysis integration. Fig. 4 illustrates some examples of the robot deployed
in the field.

###  IV-B Deployment methodology

All the missions followed the same procedure. A _robot operator_ was in charge
of defining the survey mission and supervising the overall execution. A
_safety operator_ , carrying a remote controller to overrule the autonomy
system, followed the robot at a safe distance during the mission to avoid
undesired situations (e.g. collisions and _cul-de-sacs_).

We logged all the signals from the robot and the autonomy system, as well as
the actions executed by the safety operator. The safety interventions were the
main signal we used for evaluating the performance of the autonomy system, as
it enables the use of metrics such as the mdbi and mtbi, as done in autonomous
driving literature [51]. However, instead of computing the metrics using the
total distance traveled over the number of interventions reported, we
determined the average by measuring the actual time and distance intervals
when the robot was effectively walking autonomously. We call these
measurements _autonomy segments_ , shown in blue in the summary plot (Fig. 5).

Figure 5: Summary of the seven survey missions executed across the three
campaigns. We illustrate the periods where the robots operated fully
autonomously, and where they required manual interventions from the safety
operator.

###  IV-C Evo Campaign

Our first campaign, in Finland, had the goal of evaluating the proposed
autonomy system in forest environments. This was executed with an ANYmal C
platform, using the Velodyne VLP-16
($30\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$ fov) as the main
environmental sensor. We used a state estimation setup derived from the
CERBERUS SubT stack [52], which used CompSLAM [53] as the main odometry
system. CompSLAM does not perform online pose graph optimization and loop-
closure detection to ensure consistency.

Fig. 5 shows the autonomy performance reported in these missions (M1-M5),
while Fig. 4 (a) shows an example from M2. Missions M1 and M2 were executed in
the same plots, obtaining consistent results. Mission M3 was executed in an
area nearby, but with more challenging terrain, which required to interrupt
the mission due to the robot getting trapped in a damp area (Fig. 6). Missions
M4 and M5 were executed in different areas of the forest, but achieving
similar area coverage per time unit ( $1.5\text{\,}\mathrm{ha}\text{/}$), as
shown in Tab. I. In M5, we additionally observed severe tracking problems with
the reference survey path due to drift in the state estimation system, which
were more evident in this longer sequence ($600\text{\,}\mathrm{m}\text{/}$
long, compared to $300\text{\,}\mathrm{m}\text{/}$ set for the other
missions).

Figure 6: Failure case in mission M3 in Evo, Finland. The robot was trapped
due to damp terrain, which was not detected by the onboard perception systems.

The missions lasted between $1020\text{\,}\mathrm{min}\text{/}$, with the
robot walking up to $600\text{\,}\mathrm{m}\text{/}$ on the longest missions.
About the autonomy performance, while we report the mdbi and mtbi in Tab. I,
we found the fine grained analysis of the distribution of autonomy segments a
more informative metric. This is shown in Fig. 7.

Figure 7: Distribution of distance and time between interventions for all
missions. _Top:_ Distance between interventions. _Bottom:_ Time between
interventions.

With exception of mission M2, which was fully autonomous, in most of the other
missions we reported a higher density of short autonomy segments
($\sim$$10\text{\,}\mathrm{m}\text{/}$ or
$\sim$$50\text{\,}\mathrm{s}\text{/}$). Many cases were due to short, frequent
‘pushes’ performed by the safety operator to move the local planner solution
out of local minima, or driving the robot around a dead end.

We additionally reported how long the interventions took on each mission (Fig.
8). We observed that most of the interventions were
$\sim$$15\text{\,}\mathrm{s}\text{/}$, which is within the time required for
the safety operator to get the robot out of dead ends that the mission and
local planner interaction was not able to solve automatically. Only in
specific cases—M3 and M5—, where the robot was deployed in considerably more
challenging environments (damp terrain, short trees, higher density of bushes)
we required longer interventions ($>$$1\text{\,}\mathrm{min}\text{/}$) to move
the robot to safer areas before continuing the mission.

Figure 8: Duration of the interventions reported for each mission. Most of the
interventions reported were short—$20\text{\,}\mathrm{s}\text{/}$ on average.

###  IV-D Wytham Woods Campaign

The second campaign was executed in the UK, in a mixed forest plot in Wytham
Woods, near Oxford. The same ANYmal C platform and LiDAR setup was used in
this campaign, as the goal was to improve the state estimation solution
problems reported in the Evo campaign. For this, we transitioned to the
proposed state estimation solution, which explicitly relied on a LiDAR-
inertial odometry system along with an online pose graph SLAM to ensure
consistency.

Fig. 4 (b) shows the testing environment, which was a clear sloped plot, with
a few loose branches and twigs. We set a smaller
$20\text{\,}\mathrm{m}\text{/}$ $\times$ $20\text{\,}\mathrm{m}\text{/}$
survey area, which corresponded to $0.29\text{\,}\mathrm{ha}\text{/}$ covered
for an effective $15\text{\,}\mathrm{m}\text{/}$ LiDAR range. The reported
coverage rate was consistent with the results obtained in the Evo campaign,
above $1.5\text{\,}\mathrm{ha}\text{/}\mathrm{h}$. Only two short
interventions were required (see Fig. 8, M6) to avoid falls due to loose
branches, which trapped the robot’s legs.

###  IV-E Forest of Dean Campaign

The last campaign, executed in the Forest of Dean in the UK, was used to test
the integration of the forest analysis system towards online forest inventory.
Additionally, we changed the robot platform from ANYmal C to ANYmal D—which
has longer shanks, though the locomotion policy does not exploit this fact
explicitly. We also improved the sensor payload, where we used a sensing unit
with a dedicated computer and a Hesai QT64 LiDAR, which as a larger vertical
fov ($104\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$), enabling better
coverage of the full trees. The state estimation and autonomy were left
unchanged apart from the sensor configuration, as they were validated in the
previous campaigns.

For the specific mission we executed, denoted M7, we specified a
$125\text{\,}\mathrm{m}\text{/}$ $\times$ $30\text{\,}\mathrm{m}\text{/}$
survey area, which effectively corresponded to a
$0.96\text{\,}\mathrm{ha}\text{/}$ plot. The robot completed the mission in
approximately $20\text{\,}\mathrm{min}\text{/}$, segmenting up to $100$ trees,
as shown in Fig. 9. Fig. 10 shows the marteloscope, which summarizes the tree
positions and dbh reported for the trees successfully reconstructed by our
proposed online forest analysis pipeline.

Figure 9: Illustrative example of the output produced by our system after an
autonomous mission in the Forest of Dean, UK. The robot surveyed
$0.96\text{\,}\mathrm{ha}\text{/}$ in $20\text{\,}\mathrm{min}\text{/}$, and
around $100$ trees were segmented during online operation. Figure 10:
Marteloscope produced after the deployment in the forest of Dean, indicating
the tree positions and the dbh of the successfully reconstructed trees.

##  V Current Challenges and Future Work

Our deployments so far have allowed us to continuously refine and improve our
proposed forest inventory solution, from a software and hardware perspective.
Our experiments to date indicate that, given the robot capabilities and the
sensing setup, our solution should enable the autonomous inventory at
$2\text{\,}\mathrm{ha}\text{/}\mathrm{h}$. This is a feasible bound given the
battery specifications of the ANYmal platform, assuming fairly clean forest
areas without challenging ground vegetation. Nevertheless, we still observe
open challenges to address to improve the performance of our systems towards
autonomous inventory and monitoring.

Autonomy performance: Firstly, we aim to improve the performance reported
through the mdbi and mtbi metrics. Across the missions we did not observe a
clear trend of improvement, partly caused by the changes in the constituent
components but also due to measuring the metrics in completely different
environments. To improve performance, we identified potential improvements in
the autonomy system. Particularly by designing new re-planning strategies to
detect and navigate around untraversable areas effectively, identifying and
avoiding dead ends, as well as maximizing loop closure corrections.

Locomotion and navigation: In spite of using a robust locomotion controller,
we observed that further improvements are still required to negotiate branches
and twigs in a safer manner, to avoid potential falls and failures due to
trapped legs. This could be solved by low-level strategies to vary the step
size, or by identifying such areas with the onboard perception systems,
effectively guiding the robot through longer-though-safer paths.

Forestry tasks integration: Lastly, we identified further improvements can be
achieved for the real-world usability of our solution. Practical issues such
as improving long-range communication in the forest, could enable better
feedback of the operation, and provide live reporting of the estimated traits
for the foresters. Additionally, studying the impact of legged locomotion on
the forest soil, as well as practical comparisons to alternative inventory
approaches could further contribute to the adoption of these platforms for
forestry.

##  VI Conclusion

In this work, we presented ongoing efforts to develop an autonomous forest
inventory system with small-scale legged platforms. Through seven different
missions in forests in Finland and the UK, we have demonstrated the
feasibility of using legged platforms for these tasks, providing a systematic
way to survey and study natural environments. We have also presented open
challenges, which we identified during our deployments. While some of them
will be improved for future field deployments, they also present research
questions that can be of interest to the field and legged robotics
communities.

##  VII Acknowledgments

This work is supported by EU Horizon 2020 (101070405 DigiForest) as well as a
Royal Society University Research Fellowship. We acknowledge PreFor Oy for
organizing the Evo campaign, as well as Forest Research UK for arranging the
Forest of Dean campaign.

## References

  * [1] C. L. Brack, “Forest measurement and modeling - measuring trees, stands and forests for effective forest management,” Australian National University, Tech. Rep., 2001. [Online]. Available: http://fennerschool-associated.anu.edu.au/mensuration/index.htm
  * [2] G. J. Hamilton, _Forest mensuration handbook_. Forestry Commission, 1988.
  * [3] R. O. Dubayah and J. B. Drake, “Lidar Remote Sensing for Forestry,” _Journal of Forestry_ , vol. 98, no. 6, pp. 44–46, June 2000.
  * [4] R. J. Hall, “The Roles of Aerial Photographs in Forestry Remote Sensing Image Analysis,” in _Remote Sensing of Forest Environments: Concepts and Case Studies_ , M. A. Wulder and S. E. Franklin, Eds. Boston, MA: Springer US, 2003, pp. 47–75.
  * [5] X. Liang, V. Kankare, J. Hyyppä, Y. Wang, A. Kukko, H. Haggrén, X. Yu, H. Kaartinen, A. Jaakkola, F. Guan, M. Holopainen, and M. Vastaranta, “Terrestrial laser scanning in forest inventories,” _ISPRS J. Photogramm. Remote Sens._ , vol. 115, pp. 63–77, 2016.
  * [6] M. Miettinen, M. Ohman, A. Visala, and P. Forsman, “Simultaneous localization and mapping for forest harvesters,” in _IEEE Int. Conf. Robot. Autom._ , 2007, pp. 517–522.
  * [7] J.-F. Tremblay, M. Béland, R. Gagnon, F. Pomerleau, and P. Giguère, “Automatic three-dimensional mapping for tree diameter measurements in inventory operations,” _J. Field Robot._ , vol. 37, no. 8, pp. 1328–1346, 2020.
  * [8] L. P. Koh and S. A. Wich, “Dawn of Drone Ecology: Low-Cost Autonomous Aerial Vehicles for Conservation,” _Tropical Conservation Science_ , vol. 5, no. 2, pp. 121–132, June 2012.
  * [9] D. R. A. Almeida, A. M. Almeyda Zambrano, E. N. Broadbent, A. L. Wendt, P. Foster, B. E. Wilkinson, C. Salk, D. D. A. Papa, S. C. Stark, R. Valbuena, E. B. Gorgens, C. A. Silva, P. H. S. Brancalion, M. Fagan, P. Meli, and R. Chazdon, “Detecting successional changes in tropical forest structure using GatorEye drone‐borne lidar,” _Biotropica_ , vol. 52, no. 6, pp. 1155–1167, Nov. 2020.
  * [10] S. W. Chen, G. V. Nardari, E. S. Lee, C. Qu, X. Liu, R. A. F. Romero, and V. Kumar, “Sloam: Semantic lidar odometry and mapping for forest inventory,” _IEEE Robot. Autom. Lett._ , vol. 5, no. 2, pp. 612–619, 2020\.
  * [11] X. Liu, G. V. Nardari, F. C. Ojeda, Y. Tao, A. Zhou, T. Donnelly, C. Qu, S. W. Chen, R. A. F. Romero, C. J. Taylor, and V. Kumar, “Large-scale autonomous flight with real-time semantic slam under dense forest canopy,” _IEEE Robot. Autom. Lett._ , vol. 7, no. 2, pp. 5512–5519, 2022.
  * [12] M. Pierzchała, P. Giguère, and R. Astrup, “Mapping forests using an unmanned ground vehicle with 3D LiDAR and graph-SLAM,” _Comput. Electron. Agr._ , vol. 145, pp. 217–225, Feb. 2018.
  * [13] D. J. Todd, _Applications of walking machines_. Boston, MA: Springer US, 1985, pp. 169–177.
  * [14] E. Jelavic, D. Jud, P. Egli, and M. Hutter, “Robotic precision harvesting: Mapping, localization, planning and control for a legged tree harvester,” _Field Robotics_ , 2021.
  * [15] T. Miki, J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning robust perceptive locomotion for quadrupedal robots in the wild,” _Sci. Robot._ , vol. 7, no. 62, 2022.
  * [16] J. Tang, Y. Chen, A. Kukko, H. Kaartinen, A. Jaakkola, E. Khoramshahi, T. Hakala, J. Hyyppä, M. Holopainen, and H. Hyyppä, “SLAM-Aided Stem Mapping for Forest Inventory with Small-Footprint Mobile LiDAR,” _Forests_ , vol. 6, no. 12, pp. 4588–4606, Dec. 2015\.
  * [17] A. Proudman, M. Ramezani, S. T. Digumarti, N. Chebrolu, and M. Fallon, “Towards real-time forest inventory using handheld lidar,” _Robotics and Autonomous Systems_ , vol. 157, p. 104240, 2022.
  * [18] L. Freißmuth, M. Mattamala, N. Chebrolu, S. Schaefer, S. Leutenegger, and M. Fallon, “Online Tree Reconstruction and Forest Inventory on a Mobile Robotic System,” 2024, (Under review).
  * [19] K. Anderson and K. J. Gaston, “Lightweight unmanned aerial vehicles will revolutionize spatial ecology,” _Frontiers in Ecology and the Environment_ , vol. 11, no. 3, pp. 138–146, 2013.
  * [20] J. R. Kellner, J. Armston, M. Birrer, K. C. Cushman, L. Duncanson, C. Eck, C. Falleger, B. Imbach, K. Král, M. Krůček, J. Trochta, T. Vrška, and C. Zgraggen, “New Opportunities for Forest Remote Sensing Through Ultra-High-Density Drone Lidar,” _Surveys in Geophysics_ , vol. 40, no. 4, pp. 959–977, July 2019.
  * [21] K. Shah, G. Ballard, A. Schmidt, and M. Schwager, “Multidrone aerial surveys of penguin colonies in antarctica,” _Sci. Robot._ , vol. 5, no. 47, p. eabc3000, 2020.
  * [22] D. C. Schedl, I. Kurmi, and O. Bimber, “An autonomous drone for search and rescue in forests using airborne optical sectioning,” _Sci. Robot._ , vol. 6, no. 55, p. eabg1188, 2021.
  * [23] X. Zhou, X. Wen, Z. Wang, Y. Gao, H. Li, Q. Wang, T. Yang, H. Lu, Y. Cao, C. Xu, and F. Gao, “Swarm of micro flying robots in the wild,” _Sci. Robot._ , vol. 7, no. 66, p. eabm5954, 2022.
  * [24] Y. Lin, J. Hyyppa, and A. Jaakkola, “Mini-UAV-Borne LIDAR for Fine-Scale Mapping,” _IEEE Geosci. Remote Sens. Lett._ , vol. 8, no. 3, pp. 426–430, May 2011.
  * [25] E. Hyyppä, J. Hyyppä, T. Hakala, A. Kukko, M. A. Wulder, J. C. White, J. Pyörälä, X. Yu, Y. Wang, J.-P. Virtanen, O. Pohjavirta, X. Liang, M. Holopainen, and H. Kaartinen, “Under-canopy UAV laser scanning for accurate forest field measurements,” _ISPRS J. Photogramm. Remote Sens._ , vol. 164, pp. 41–60, June 2020.
  * [26] E. Aucone, S. Kirchgeorg, A. Valentini, L. Pellissier, K. Deiner, and S. Mintchev, “Drone-assisted collection of environmental DNA from tree branches for biodiversity monitoring,” _Sci. Robot._ , vol. 8, no. 74, p. eadd5762, 2023.
  * [27] A. Tominaga, H. Eiji, and A. Mowshowitz, “Development of Navigation System in Field Robot for Forest Management,” in _Int. Conf. Soft Comput. Intell. Syst. and Int. Sym. Adv. Intell. Syst._ , Toyama, Japan, Dec. 2018, pp. 1142–1147.
  * [28] J. Rossmann, M. Schluse, C. Schlette, A. Buecken, P. Krahwinkler, and M. Emde, “Realization of a highly accurate mobile robot system for multi purpose precision forestry applications,” in _2009 International Conference on Advanced Robotics_ , June 2009, pp. 1–6.
  * [29] E. Jelavic, T. Kapgen, S. Kerscher, D. Jud, and M. Hutter, “Harveri : A Small (Semi-)Autonomous Precision Tree Harvester,” in _Innovation in Forestry Robotics: Research and Industry Adoption (ICRA 2022)_ , May 2022.
  * [30] S. Fountas, T. A. Gemtos, and S. Blackmore, _Robotics and Sustainability in Soil Engineering_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2010, pp. 69–80.
  * [31] A. Calleja-Huerta, M. Lamandé, O. Green, and L. Munkholm, “Vertical and horizontal stresses from a lightweight autonomous field robot during repeated wheeling,” _Soil and Tillage Research_ , vol. 233, p. 105790, Sept. 2023.
  * [32] C. G. Wilson, F. Qian, D. J. Jerolmack, S. Roberts, J. Ham, D. Koditschek, and T. F. Shipley, “Spatially and temporally distributed data foraging decisions in disciplinary field science,” _Cognitive Research: Principles and Implications_ , vol. 6, no. 1, p. 29, 2021.
  * [33] K. Melo, T. Horvat, and A. J. Ijspeert, “Animal robots in the african wilderness: Lessons learned and outlook for field robotics,” _Sci. Robot._ , vol. 8, no. 85, p. eadd8662, 2023.
  * [34] M. Tranzatto, M. Dharmadhikari, L. Bernreiter, M. Camurri, S. Khattak, F. Mascarich, P. Pfreundschuh, D. Wisth, S. Zimmermann, M. Kulkarni, V. Reijgwart, B. Casseau, T. Homberger, P. D. Petris, L. Ott, W. Tubby, G. Waibel, H. Nguyen, C. Cadena, R. Buchanan, L. Wellhausen, N. Khedekar, O. Andersson, L. Zhang, T. Miki, T. Dang, M. Mattamala, M. Montenegro, K. Meyer, X. Wu, A. Briod, M. Mueller, M. Fallon, R. Siegwart, M. Hutter, and K. Alexis, “Team CERBERUS Wins the DARPA Subterranean Challenge: Technical Overview and Lessons Learned,” 2022.
  * [35] D. Wooden, M. Malchano, K. Blankespoor, A. Howard, A. A. Rizzi, and M. H. Raibert, “Autonomous navigation for bigdog,” in _IEEE Int. Conf. Robot. Autom._ , 2010, pp. 4736–4741.
  * [36] D. M. Bradley, J. K. Chang, D. Silver, M. Powers, H. Herman, P. Rander, and A. Stentz, “Scene understanding for a high-mobility walking robot,” in _IEEE/RSJ Int. Conf. Intell. Robots Syst._ , 2015, pp. 1144–1151.
  * [37] F. Qian, D. Jerolmack, N. Lancaster, G. Nikolich, P. Reverdy, S. Roberts, T. Shipley, R. S. Van Pelt, T. M. Zobeck, and D. E. Koditschek, “Ground robotic measurement of aeolian processes,” _Aeolian Research_ , vol. 27, pp. 1–11, 2017.
  * [38] F. Qian, D. Lee, G. Nikolich, D. Koditschek, and D. Jerolmack, “Rapid in situ characterization of soil erodibility with a field deployable robot,” _Journal of Geophysical Research: Earth Surface_ , vol. 124, no. 5, pp. 1261–1280, 2019.
  * [39] P. Arm, G. Waibel, J. Preisig, T. Tuna, R. Zhou, V. Bickel, G. Ligeza, T. Miki, F. Kehl, H. Kolvenbach, and M. Hutter, “Scientific exploration of challenging planetary analog environments with a team of legged robots,” _Sci. Robot._ , vol. 8, no. 80, p. eade9548, 2023.
  * [40] D. Wisth, M. Camurri, and M. Fallon, “Vilens: Visual, inertial, lidar, and leg odometry for all-terrain legged robots,” _IEEE Trans. Robot._ , vol. 39, no. 1, pp. 309–326, 2023.
  * [41] P. J. Besl and N. D. McKay, “A method for registration of 3-d shapes,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 14, no. 2, pp. 239–256, 1992.
  * [42] C. Forster, L. Carlone, F. Dellaert, and D. Scaramuzza, “On-manifold preintegration for real-time visual–inertial odometry,” _IEEE Trans. Robot._ , vol. 33, no. 1, pp. 1–21, 2017.
  * [43] M. Ramezani, G. Tinchev, E. Iuganov, and M. Fallon, “Online LiDAR-SLAM for Legged Robots with Robust Registration and Deep-Learned Loop Closure,” in _IEEE Int. Conf. Robot. Autom._ , 2020, pp. 4158–4164.
  * [44] M. Kaess, H. Johannsson, R. Roberts, V. Ila, J. J. Leonard, and F. Dellaert, “iSAM2: Incremental smoothing and mapping using the Bayes tree,” _Int. J. Robot. Res._ , vol. 31, no. 2, pp. 216–235, 2012.
  * [45] P. Fankhauser, M. Bloesch, and M. Hutter, “Probabilistic terrain mapping for mobile robots with uncertain localization,” _IEEE Robot. Autom. Lett._ , vol. 3, no. 4, pp. 3019–3026, 2018.
  * [46] H. Choset, “Coverage of known spaces: The boustrophedon cellular decomposition,” _Auton. Robots_ , vol. 9, no. 3, pp. 247–253, 2000.
  * [47] J. Frey, M. Mattamala, N. Chebrolu, C. Cadena, M. Fallon, and M. Hutter, “Fast Traversability Estimation for Wild Visual Navigation,” in _Robot.: Sci. Syst._ , 2023.
  * [48] M. Mattamala, N. Chebrolu, and M. F. Fallon, “An efficient locally reactive controller for safe navigation in visual teach and repeat missions,” _IEEE Robot. Autom. Lett._ , vol. 7, no. 2, pp. 2353–2360, 2022.
  * [49] W. Zhang, J. Qi, P. Wan, H. Wang, D. Xie, X. Wang, and G. Yan, “An easy-to-use airborne lidar data filtering method based on cloth simulation,” _Remote. Sens._ , vol. 8, no. 6, p. 501, 2016.
  * [50] C. Cabo, C. Ordóñez, C. A. López-Sánchez, and J. Armesto, “Automatic dendrometry: Tree detection, tree height and diameter estimation using terrestrial laser scanning,” _Int. J. Appl. Earth Obs. Geoinformation_ , vol. 69, pp. 164–174, 2018.
  * [51] D. Paz, P. Lai, N. Chan, Y. Jiang, and H. I. Christensen, “Autonomous vehicle benchmarking using unbiased metrics,” in _IEEE/RSJ Int. Conf. Intell. Robots Syst._ , 2020, pp. 6223–6228.
  * [52] M. Tranzatto, T. Miki, M. Dharmadhikari, L. Bernreiter, M. Kulkarni, F. Mascarich, O. Andersson, S. Khattak, M. Hutter, R. Siegwart, and K. Alexis, “CERBERUS in the DARPA subterranean challenge,” _Sci. Robot._ , vol. 7, no. 66, 2022.
  * [53] S. Khattak, H. Nguyen, F. Mascarich, T. Dang, and K. Alexis, “Complementary multi–modal sensor fusion for resilient robot pose estimation in subterranean environments,” in _Int. Conf. on Unmanned Aircraft Systems (ICUAS)_ , 2020, pp. 1024–1029.