# Broadening and redward asymmetry of H$\alpha$ line profiles observed by
LAMOST during a stellar flare on an M-type star

Yuchuan Wu  School of Earth and Space Sciences, Peking University, Beijing,
100871, China.  Hechao Chen  School of Earth and Space Sciences, Peking
University, Beijing, 100871, China.  Hui Tian  School of Earth and Space
Sciences, Peking University, Beijing, 100871, China.  National Astronomical
Observatories, Chinese Academy of Sciences, Beijing, 100101, China.  Liyun
Zhang  College of Physics, Guizhou University, Guiyang 550025, China.  Key
Laboratory for the Structure and Evolution of Celestial Objects, Chinese
Academy of Sciences, Kunming 650011, China.  Jianrong Shi  National
Astronomical Observatories, Chinese Academy of Sciences, Beijing, 100101,
China.  Han He  CAS Key Laboratory of Solar Activity, National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China.  University
of Chinese Academy of Sciences, Beijing 100049, China.  Hongpeng Lu  School of
Earth and Space Sciences, Peking University, Beijing, 100871, China.  Yu Xu
School of Earth and Space Sciences, Peking University, Beijing, 100871, China.
Haifeng Wang  GEPI, Observatoire de Paris, Université PSL, CNRS, Place Jules
Janssen 92195, Meudon, France Hechao Chen; Hui Tian<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Stellar flares are characterized by sudden enhancement of electromagnetic
radiation in stellar atmospheres. So far much of our understanding of stellar
flares comes from photometric observations, from which plasma motions in flare
regions could not be detected. From the spectroscopic data of LAMOST DR7, we
have found one stellar flare that is characterized by an impulsive increase
followed by a gradual decrease in the H$\alpha$ line intensity on an M4-type
star, and the total energy radiated through H${\alpha}$ is estimated to be on
the order of $10^{33}$ erg. The H$\alpha$ line appears to have a Voigt profile
during the flare, which is likely caused by Stark pressure broadening due to
the dramatic increase of electron density and/or opacity broadening due to the
occurrence of strong non-thermal heating. Obvious enhancement has been
identified at the red wing of the H$\alpha$ line profile after the impulsive
increase of the H$\alpha$ line intensity. The red wing enhancement corresponds
to plasma moving away from the Earth at a velocity of 100$-$200 km s-1.
According to the current knowledge of solar flares, this red wing enhancement
may originate from: (1) flare-driven coronal rain, (2) chromospheric
condensation, or (3) a filament/prominence eruption that either with a non-
radial backward propagation or with strong magnetic suppression. The total
mass of the moving plasma is estimated to be on the order of $10^{15}$ kg.

Stellar flares (1603); Late-type dwarf stars (906); Stellar coronal mass
ejections (1881); Time domain astronomy (2109); Spectroscopy (1558)

##  1 Introduction

Stellar flares are dramatic explosions in stellar atmospheres, which are often
characterized by a rapid rise followed by a slow decrease in the
electromagnetic radiation across a wide range of wavelength from radio to
$\gamma$ rays. On late-type main-sequence stars, these flare events are
generally thought to be triggered by sudden release of magnetic energy in the
magnetized stellar atmospheres (Haisch et al., 1991). According to the widely
accepted solar flare scenario (e.g., Benz & Güdel, 2010; Shibata & Magara,
2011), magnetic reconnection and its associated rapid energy release in the
corona play a key role in triggering solar flares. During this process an
enormous amount of energy is transported to the lower solar atmosphere and
dissipates in the dense chromosphere via thermal conduction and thermalization
of (relativistic) energetic electrons. This leads to intense localized plasma
heating, which in turn results in chromospheric evaporation and condensation
mainly in the impulsive phase of flares (e.g., Fletcher et al., 2011). In the
meantime, continuous increase of the local temperature, electron density and
plasma velocity would cause strong enhancement and broadening of chromospheric
spectral lines such as the H$\alpha$ line. Hence, spectroscopic observations
of chromospheric lines have long been important diagnostics for the physical
processes behind solar and stellar flares.

After decades of multi-wavelength imaging observations of solar flares, solar
physicists have well realized that a solar flare is only part of a large-scale
eruptive energy releasing process, in which the flare occurs low in the
atmosphere and a coronal mass ejection (CME) may take place at a much larger
scale (e.g., Lin et al., 2003; Zhang & Low, 2005; Shibata & Magara, 2011;
Chen, 2011; Osten & Wolk, 2017). In fact, CMEs play the most important role in
driving hazardous space weather in the near-Earth space environment (Gosling,
1993). Statistically, flares and CMEs on the Sun are more closely associated
with each other with increasing flare class (e.g., Yashiro & Gopalswamy, 2009;
Aarnio et al., 2011). Recent photometric observations from the Kepler
spacecraft and the Transiting Exoplanet Survey Satellite (TESS) revealed a lot
of superflares on G-, K-, and M-type stars (e.g., Maehara et al., 2012; Notsu
et al., 2013; He et al., 2015, 2018; Namekata et al., 2017; Doyle et al.,
2020; Tu et al., 2020, 2021; Yan et al., 2021). The total released energy of
one such superflare ($10^{33}-10^{38}$ erg) often far exceeds that of the
strongest solar flare (with energy up to $10^{32}$ erg)(Schaefer et al., 2000;
Shibata et al., 2013), indicating that they are possibly associated with
monster CMEs. Studying stellar flares and associated CMEs of these late-type
stars, especially active M dwarfs, are of particular interest. This is because
most known potentially habitable exoplanets are orbiting around these stars
(Alvarado-Gómez et al., 2021). If frequent high-energy flares and associated
CMEs take place on these stars, they may impact or even threaten the
habitability of nearby exoplanets (Khodachenko et al., 2007; Yelle et al.,
2008). However, detection of CMEs on these stars is still challenging, since
direct imaging of stellar coronae is almost impossible even for the nearest
stars (see a brief review by Moschou et al., 2019).

On the Sun, flare processes can be spatially resolved in great detail. As a
result, characteristics of the line broadening and asymmetry during different
phases of solar flares have been well studied through many observations (e.g.,
Švestka, 1972; Schmieder et al., 1987; Zarro et al., 1988). For instance, a
symmetric chromospheric line broadening is often believed to be caused by the
Stark effect or turbulence (e.g., Kowalski et al., 2017; Zhu et al., 2019).
While an asymmetric line broadening is usually explained by bulk plasma motion
in the form of chromospheric evaporation (e.g., Canfield et al., 1990; Li &
Ding, 2011; Li et al., 2021), chromospheric condensation (e.g., Ichimoto &
Kurokawa, 1984; Fisher et al., 1985; Gan et al., 1993; Ding et al., 1995; Tian
et al., 2015; Zhang et al., 2016; Tei et al., 2018; Yu et al., 2020) or
coronal rain (e.g., Antolin, 2020; Li et al., 2020; Chen et al., 2021; Schad
et al., 2021).

As an extension of solar spectroscopic studies, spectroscopic observations of
chromospheric lines on other stars have not only been used to understand the
dynamics of stellar flares, but also been used to detect flare-associated
CMEs, since CMEs generally host a chromosphere-temperature dense core, i.e.,
the erupting cool prominence/filament (e.g., Illing & Hundhausen, 1985). For
instance, Houdebine et al. (1990) found an obvious enhancement in the blue
wings of Balmer lines during a flare occurring on the well-known flaring dMe
dwarf, AD Leo. The blueward asymmetry corresponds to a high speed of 5800 km
s-1, which has been explained as a possible CME. With spectroscopic
observations, Crespo-Chacón et al. (2006) conducted a systematic flare study
on AD Leo, in which three strongest ones are found to be associated with
Balmer line broadening and two show weak redward asymmetries. Vida et al.
(2016) reported a complex CME with falling-back and re-ejected material on the
M4 dwarf V374 Peg, in which its associated blue shift (675 km s-1) is above
the stellar surface escape velocity ($\sim$580 km s-1). From high-resolution
spectra of 28 active M dwarfs, Fuhrmeister et al. (2018) identified broad,
potentially asymmetric, wings in H$\alpha$ during stellar flares. They
explained the observed red asymmetries as a process similar to chromospheric
downward condensation or coronal rain, and suggested that the line broadening
may be caused by the Stark effect. Recently, several authors have tried to
search for CMEs on some flaring M dwarfs with time-resolved spectra (Muheki et
al., 2020; Maehara et al., 2021). Two candidates of erupting
filaments/prominences have been reported on EV Lac and YZ CMi, respectively.
With photometric and follow-up spectroscopic observations, Wang et al. (2021)
reported detection of possible CMEs during flares of two M-dwarfs. In
addition, focusing on more than 5500 spectra of M-dwarf stars, Vida et al.
(2019) found 478 spectra with Balmer-line asymmetries, most of which are
associated with flux enhancement and their maximum velocities are on the order
of 100$-$300 km s-1. Similarly, Muheki et al. (2020) and Koller et al. (2021)
respectively used a large number of spectra to search for flares and CMEs on M
dwarfs, but consistently reported a high level of flare activity and a low
detection rate of CME candidates. More recently, Namekata et al. (2022)
successfully detected a definite evidence of a flare-filament eruption event
on the G-type star EK Draconis using time-resolved H${\alpha}$ spectroscopic
observations; moreover, they compared this stellar eruption signature with
solar filament eruptions in the Sun-as-a-star spectral view, revealing that
the former is the same as the magnified version of the latter.

From the spectra provided by the Large Sky Area Multi-Object Fibre
Spectroscopic Telescope (LAMOST), we have found one interesting stellar flare
with obviously broadened and asymmetric H$\alpha$ line profiles. In Section 2,
we describe the observational data and show analysis results of the event,
including the Voigt-Gauss function fitting and red-blue asymmetry analysis. In
Section 3, we estimate the flare energy and mass of the moving plasma, and
provide explanations for the Lorentz profile and asymmetry of H$\alpha$.
Finally, conclusions are given in Section 4.

##  2 Observation and data analysis results

LAMOST, also called the Guoshoujing Telescope, is a medium-aperture ( 4.6 m)
optical telescope in China, with the capability of performing spectral surveys
of the northern-sky stars in the Milky Way and other objects (Su et al., 1998;
Cui et al., 2012). In recent years, LAMOST provides spectra in time domain,
making it possible to search for stellar flares through the time evolution of
H$\alpha$ line intensity and profile asymmetry. Here we used data from LAMOST
DR7 to search for H$\alpha$ line asymmetries. In LAMOST DR7 medium-resolution
(R$\sim$7500) spectral catalogue, there are 2 828 655 time-domain spectra,
covering a wavelength range of 6300$-$6800 Å. These spectra were taken between
2017 September 28 and 2019 June 8.

Figure 1: H$\alpha$ , He i and Mg i line profiles observed at different times
during the flare of interest. The profiles are shifted in the y-axis for a
better illustration.

The stellar flare we study here was observed on 2019 January 16, and the
starting observation time is 23:37 (local time) , or 15:37 (UT). The obsid is
716907097. The original time-domain spectra of the H$\alpha$ , He i 6680 Å and
Mg i 5171 Å lines are shown in Figure 1. Each exposure lasts for 20 minutes,
and the time cadence of exposures is around 23 minutes. From Figure 1, we can
see that obvious radiation enhancement and profile broadening of H$\alpha$
line commence at 70 mins since the starting time and last to the end of this
observation. Other chromospheric lines like He i and Mg i become visible
during the flare, and their intensities decrease in the later phase of the
flare. Interestingly, a red wing enhancement of the H$\alpha$ line was also
presented in the later phase of the flare, indicating the presence of
materials (plasma) moving away from the Earth. To get a more clear view of the
red wing enhancement, please see Figure 2 in Section 2.1. LAMOST provides
Signal-to-noise ratios (SNRs) that were calculated by the equation of
$flux*(inverse\ variance)^{0.5}$. For the eight spectra of the star, their
SNRs are 10.02, 9.96, 10.13, 27.04, 23.64, 13.27, 11.62 and 10.13,
respectively. Note that after taking the Earth’s motion and the star’s
redshift into consideration, the corrected theoretical line center of
H${\alpha}$ profile is 6564.61 Å.

According to LAMOST catalog
111http://dr7.lamost.org/v2.0/spectrum/view?obsid=432616084, this target is an
M4-type star, and other relevant data are shown in Table 1. The star was also
observed by the Kepler’s second mission (Howell et al., 2014), and several
flares were detected on it 222The K2 observations of this star can be found
from https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html or
https://drive.google.com/drive/folders/1K46MbuEarPNwNwRIr0n8dR7eCppkYzDy?usp=sharing.
We used the effective temperature to estimate the mass and radius of the star.
We can find the average radius, mass and effective temperature of an M2 or M5
type main sequence star from Cox & Pilachowski (2000). We then assumed a
linear relationship of $lg(\frac{M}{M_{\sun}})$ and $lg(\frac{R}{R_{\sun}})$
versus $lgT_{\rm eff}$ for stars with a spectral type between M2 and M5. Here
$M$, $M_{\sun}$, $R$ and $R_{\sun}$ represent the masses and radii of the star
and our Sun, respectively, and $T_{\rm eff}$ is the effective temperature of a
star. Since the effective temperature of the observed star is 3404 K, we
obtained the following,

$\frac{M}{M_{\sun}}=0.326,\frac{R}{R_{\sun}}=0.411.$ (1)

Then we can calculate the escape velocity of the star, which is 550 km s-1.

Table 1:  Basic information of the star (716907097) Right ascension | Declination | Subclass | Redshift | Effective temperature | Surface gravity | Magnitudes | Magnitudes | Magnitudes
---|---|---|---|---|---|---|---|---
(degree) | (degree) |  |  | (K) | (log $g/cm~{}s^{-2}$) | of g filter | of r filter | of i filter
127.52 | 16.57 | dM4 | 0.00012 | 3404 | 5.02 | 15.82 | 14.53 | 13.04

###  2.1 Normalized spectra

The continuum levels of the original H$\alpha$ spectra vary with time, e.g.,
in the fourth and fifth exposures the continuum levels are about twice as
large as in other exposures. This variation is too big to be explained by a
white-light flare. Instead, it is more likely caused by the variation of
seeing. Due to the lack of absolute flux calibration, we thus assumed that the
continuum level remains unchanged before, during, and after the flare.
However, we also point out that the H$\alpha$ flare might be associated with a
white-light flare, especially considering that the broadened H${\alpha}$
profiles observed in the flare (see Figure 2) indicate a non-thermal heating
process (see details in Section 3.3 and also Namekata et al., 2020), but we
cannot confirm whether a white-light flare exists from the available data.
Before further analysis, we need to normalize the original spectra to the
continuum intensities. To obtain the continuum intensities, we used a
summation of a Voigt function, a Gauss function (only for the last four
spectra) and a constant background to fit each H$\alpha$ profile. The constant
was taken as the continuum level. The Gauss function was introduced here to
account for the additional bulk flow component, which can be inferred from the
redward asymmetry. The Voigt function is a convolution of a Lorentz function
and a Gauss function. A Voigt function fitting has often been applied to
analyze the spectral broadening of hydrogen Balmer lines (i.e., H${\alpha}$,
Hβ) (e.g., Luque et al., 2003; Tremblay & Bergeron, 2009; Falcon et al., 2015)
and helium lines (e.g., Schad et al., 2021; Kanodia et al., 2021) in
astrophysical and laboratory plasmas. Here, we used it to fit the H$\alpha$
line profile because the chromospheric line broadening is likely dominated by
a combination of pressure broadening and Doppler broadening during flares,
which can be described as a Voigt function (Zhu et al., 2019). The total
fitting function can be written as follows (Ginsburg & Mirocha, 2011),

$I(\lambda)=\frac{Re[wofz(\frac{\lambda-\lambda_{cen}+i\gamma}{\sqrt{2}\sigma})]}{\sqrt{2\pi}\sigma}I_{0}+ae^{-(\frac{\lambda-m}{s})^{2}}+C,$
(2)

where $\lambda$ and $\lambda_{cen}$ are the wavelength and centroid wavelength
of the H$\alpha$ line, respectively. $Re$ represents the real part of a
complex number. And $wofz$ is the Faddeeva function, which satisfies
$wofz(z)=e^{-z^{2}}erfc(-iz)$.

There are 7 free parameters: $I_{0}$ characterizes the total flux of the Voigt
function; $\gamma$ is the half width at half maximum (HWHM) of the Lorentzian
component. A larger $\gamma$ means that the Voigt function is more similar to
a Lorentz profile, indicating a larger effect of pressure broadening. The
parameter $\sigma$ is the standard deviation of the the Guassian component in
the Voigt function. The height, centroid and width of the Gauss function are
represented by $a$, $m$ and $s$, respectively. The continuum level is
represented by $C$. By dividing the original spectrum taken during each
exposure by $C$, we obtained the normalized spectrum $I_{n}(\lambda)$. All
normalized spectra are shown in Figure 2.

Figure 2: H$\alpha$ line profiles normalized to the continuum intensities.

###  2.2 Voigt-Gauss function fitting

We used the same function to fit the normalized spectra. After normalization,
the fitting function becomes the following one,

$I_{n}(\lambda)=\frac{Re[wofz(\frac{\lambda-\lambda_{cen}+i\gamma}{\sqrt{2}\sigma})]}{\sqrt{2\pi}\sigma}I_{n,0}+a_{n}e^{-(\frac{\lambda-m}{s})^{2}}+1.$
(3)

The subscript $n$ means normalized intensity. The normalized spectrum, fitting
function, and the residual between fitting and observation for each exposure
are shown in Figure 3. Note that we removed the continuum just for the sake of
simplicity. The variations of $I_{0}$, widths of the Gaussian and Lorentzian
components of the Voigt-Gaussian fitting function, as well as the fitting
parameters of the Gaussian function are shown in Figure 4.

Figure 3: Normalized H$\alpha$ line profiles and their fitting results. Blue
lines are the original data, orange lines are fitting results for Voigt
function, red lines are fitting results for Gauss function, and green lines
indicate residuals between fitting results for Voigt function and observation.
Figure 4: Time evolution of the fitting parameters. The top rows show the
total flux of the Voigt function ($I_{0}$), the Half width at half maximum
(HWHM) of the Gaussian component ($\sigma\sqrt{ln4}$), and the HWHM of the
Lorentzian component ($\gamma$) in the Voigt-Gauss function, respectively. The
bottom row only shows the height ($a_{n}$), the centroid ($m$), and the HWHM
($s\sqrt{ln2}$) in the Gaussian function for the last four spectra, where
profile asymmetries are observed. In each panel, the error bars represent
uncertainties calculated from the fitting errors.

From Figures 3 and 4, we noticed the following characteristics:

(1) The total flux of the Voigt function, $I_{0}$, increases quickly between
70 min and 117 min, indicating the occurrence of a flare.

(2) The HWHM of Gaussian component in the Voigt-Gauss fitting function is
close to 0 for three spectra taken at 70 min, 93 min and 117 min. At these
occasions, the wings of the H$\alpha$ line profiles are clearly characterized
as a Lorentz function. Interestingly, this change of spectral profile starts
as soon as the rapid increase of $I_{0}$. When the rapid increase of $I_{0}$
stops, the HMHW of Gaussian component in the Voigt-Gauss fitting function
begins to increase.

(3) $a_{n}$ keeps on increasing after 93 min, meaning that the redward
asymmetry of the H$\alpha$ line profile becomes more and more significant.

(4) $m$ slowly decreases after 93 min, indicating that the moving plasma
responsible for the redward asymmetry of H$\alpha$ possibly decelerates after
93 min (see Figure 3 (e)). From 93 min to 163 min, the average deceleration
rate was estimated as 7.6 m s-2.

The time variation of the integrated H$\alpha$ flux can reflect the change of
radiation power through the H$\alpha$ line. We first integrated each
normalized spectrum from 6557 Å to 6572 Å, then remove the continuum. The
remaining integral is the equivalent width of the H$\alpha$ line. Using the
fitting parameters of the Gauss function, we can also obtain the equivalent
width of the redshifted enhancement responsible for the line asymmetry in the
last four exposures. The time variations of the equivalent widths of the
H$\alpha$ line and the redshifted enhancement are shown in Figure 5. We can
see that the the equivalent width of H$\alpha$ line increases quickly between
70 min and 117 min, and decreases slowly after 117 min, which validates the
conclusion from the analysis of the total flux of the Voigt function. At the
same time, the equivalent width of the redshifted enhancement simultaneously
grows by a factor of $\sim$6 from 93 min to 163 min.

Figure 5:  Time variations of the equivalent widths of the H$\alpha$ line and
the redshifted enhancement. The black error bars indicate uncertainties of the
equivalent width of the H$\alpha$ line by summing the measurement errors
between 6557 – 6572 Å. The blue error bars represent the uncertainties of the
equivalent width of the redshifted enhancement as calculated from the fitting
errors.

###  2.3  Red-blue asymmetry analysis

A technique called red-blue (RB) asymmetry analysis has been often used to
analyze the asymmetries of solar spectral lines both under the quiet-Sun
condition (e.g., De Pontieu et al., 2009; Tian et al., 2011) and during solar
eruptions (Tian et al., 2012). The technique can reveal where the asymmetry is
most obvious. Here we used this method to analyze the asymmetry of the
H$\alpha$ line profiles. The RB asymmetry for an offset velocity can be
expressed as following:

$RB(u_{c})=\int_{u_{c}-\delta u/2}^{u_{c}+\delta
u/2}I(u)du-\int_{-u_{c}-\delta u/2}^{-u_{c}+\delta u/2}I(u)du,$ (4)

where $u_{c}$, $\delta u$, and $I(u)$ respectively represent the velocity from
the line center, the velocity range over which the RB asymmetry is determined,
and the ratio between spectral intensity and peak intensity. Here the
H$\alpha$ line center and peak intensity are $\lambda_{cen}$ and $I_{0}$,
respectively, as calculated in Section 2.1. We found that a smaller $\delta u$
leads to too much fluctuation, and a bigger $\delta u$ ignores too much data
in the spectra. By trial and error, we found that 0.5 Å is just appropriate
for our calculation. So we let $\delta u$ be the velocity corresponding to 0.5
Å Doppler shift, which is around 22.8 km s-1. We calculated $RB$ as a function
of $u_{c}$, where $u_{c}=22.8\times(0.5+n)$ km s-1, $n=0,1,2,\cdots$.

The results of the RB asymmetry analysis for the spectra taken during the last
four exposures are shown in Figure 6. We can see that an obvious redward
asymmetry, which is indicated by the positive peak of an RB asymmetry profile,
is clearly present after 93 min, when the H$\alpha$ line intensity shows an
impulsive increase. The redward asymmetry becomes more evident towards the end
of the observation. The peaks of the RB asymmetry profiles indicate a
redshifted velocity of 100 km s-1 to 200 km s-1, far below the escape velocity
of the star (550 km s-1).

Figure 6: RB asymmetry analysis results for the H$\alpha$ spectral profiles
taken during the last four exposures. The error bars indicate the
uncertainties propagated from the measurement errors.

##  3 Discussion

###  3.1 Radiation power and total energy of the flare

We tried to estimate the radiation power of the flare. We first estimated the
luminosity of the continuum ($L_{c}$) under the H$\alpha$ line. Since the
white light continuum of a star generally can be approximated by the Planck
function, $L_{c}$ can be calculated through the following equation:

$L_{c}=4\pi^{2}R^{2}\int_{6557\rm\AA}^{6572\rm\AA}\frac{1}{e^{\frac{hc}{\lambda
kT}}-1}\frac{2hc^{2}}{\lambda^{5}}d\lambda={7.58\times 10^{28}}\ erg\ s^{-1},$
(5)

where $T$ represents the effective temperature of the star, which is 3404.72
K. $L_{c}$ represents the luminosity of the continuum between 6557 Å and 6572
Å. $R$ is the radius of the star, which is $0.411R_{\sun}$.

Considering that the equivalent width of the normalized continuum in the
wavelength range of 6557$-$6572 Å is 15 Å and the equivalent width of the
H$\alpha$ line has a peak of 28.36 Å, we can calculate the peak luminosity
($L_{f}$) of the flare through the following relation,

$\frac{L_{f}}{L_{c}}=\frac{{28.36}}{15}.$ (6)

And $L_{f}$ was estimated to be $1.43\times 10^{29}$ erg $s^{-1}$.

We can also estimate the total energy of the flare from the integrated flux
calculated in Section 2.2. This observation did not cover the whole process of
the flare. We need to remove the H$\alpha$ contribution from non-flare state
of the star. We used the average of first three spectra to estimate it, and
the result in the form of equivalent width is 5.43 Å. During the observation,
from 70 min to 163 min, the flare energy radiated through H$\alpha$ is
$4.45\times 10^{32}$ erg. If we assume that the power of the flare and the
equivalent width of H$\alpha$ in the last three spectra decay at the same
speed of the last three spectra, we estimated that the flare lasts for 418
min, and the total flare energy radiated through H$\alpha$ was estimated to be
$2.42\times 10^{33}$ erg.

###  3.2  Mass of the moving plasma

We have tried to estimate the mass of the moving plasma. We decided to use the
method which was first introduced by Houdebine et al. (1990) and recently used
by Koller et al. (2021) to estimate the mass. We used the following formula to
obtain a lower limit of the mass (Houdebine et al., 1990),

${M\geq\frac{4\pi
d^{2}F_{m}(N_{tot}/N_{j})m_{H}\eta_{OD}}{h\upsilon_{j-i}A_{j-i}}=\frac{L_{m}(N_{tot}/N_{j})m_{H}\eta_{OD}}{h\upsilon_{j-i}A_{j-i}},}$
(7)

where $d$ is the distance from the star, $F_{m}$ is the integrated excess flux
of the moving plasma. Since we did not have absolute flux calibration, we used
$L_{m}$, the luminosity of the moving plasma, to estimate $4\pi d^{2}F_{m}$.
$N_{tot}/N_{j}$ is the ratio of the number density of hydrogen atoms and the
number density of hydrogen atoms at excitation level $j$. $m_{\rm H}$ is the
mass of a hydrogen atom. $\eta_{OD}$ is the opacity damping factor. $h$ is the
Plank constant. $\upsilon_{j-i}$ and $A_{j-i}$ are the frequency and Einstein
coefficient of the spectral line, respectively.

As for other parameters in the formula above, we used the same settings as
Koller et al. (2021) once used. This flare event under study was detected from
the H$\alpha$ line (transition from j = 3 to i = 2). However, they could not
find an appropriate calculation to estimate $N_{tot}/N_{3}$, so they decided
to use the relevant data of the H$\gamma$ line (j = 5, i = 2) for their
calculation. Houdebine et al. (1990) performed non-local thermal equilibrium
(NLTE) modeling for AD Leo and estimated $N_{tot}/N_{5}=2\times 10^{9}$. They
set $\eta_{OD}$ to 2, which indicates that half of the radiation escapes the
moving feature (Leitzinger et al., 2014). For H$\gamma$, $A_{j-i}$ is
$2.53\times 10^{6}$ (Wiese & Fuhr, 2009).

We used the same method mentioned in Section 3.1 to estimate the radiation
power of the moving plasma through the red wing enhancement of the H$\alpha$
line $L_{m,\alpha}$. In Section 2.3, we have calculated the equivalent width
of the redshifted enhancement, and the equivalent width reaches the maximum in
the last exposure. The equivalent width is 3.29 Å, which means
$\frac{L_{m,\alpha}}{L_{c}}=\frac{{3.29}}{15}$. Through this way we found
$L_{m,\alpha}={1.66}\times 10^{28}$ erg s-1. To obtain the radiation power of
the H$\gamma$ line $L_{m,\gamma}$, we need to divide $L_{m,\alpha}$ by a
factor of $\sim$3 (Butler et al., 1988), just as Koller et al. (2021) did.
Finally, we obtained a lower limit for the mass of the moving plasma, which is
$3.2\times 10^{15}$ kg.

###  3.3 Lorentz profile of the H$\alpha$ line

From 70 min to 117 min, the H$\alpha$ line appears to reveal a Lorentz profile
with a symmetric line broadening. Similar line broadening phenomena have been
reported in solar flare observations (e.g., Kowalski et al., 2017; Zhu et al.,
2019) and also in stellar flare studies (Fuhrmeister et al., 2018; Namekata et
al., 2020), which have been interpreted as non-thermal broadening or Stark
(pressure) broadening. There are two classical theories for the explanation of
Stark broadening: the impact theory (Weisskopf, 1932) and the statistical
theory. The former results in a Lorentz profile and the latter leads to a
Holtsmark profile (Holtsmark, 1919). Since the observed H${\alpha}$ spectra
can be well fitted by a Voigt-Gauss function that includes a Lorentz profile,
here we take the former to explain our current observation. According to the
impact theory, each time the radiating atom collides with a charged particle,
the radiation is interrupted. After Fourier transform, these interrupted
radiation form a Lorentz profile in frequency domain. Moreover, the RADYN
numerical simulation of M-dwarf flares from Namekata et al. (2020) recently
suggested that due to the non-LTE formation properties and highly variable
opacity of the H${\alpha}$ line, the broadened H${\alpha}$ profile observed
during flares is determined not only by the Stark broadening but also
sometimes relates to “opacity broadening”. Indeed, through changing the energy
spectra from soft- into hard-energy input cases, their simulation results
revealed that once hard high-energy electron beams appear during flares and
deposit their energy in the deep chromosphere, strong self-absorptions and a
Stark effect can simultaneously cause a wider line broadening. Therefore, this
radiative-transfer-induced broadening process may also be important for the
broadened flaring H${\alpha}$ line, because stellar and solar flares often
associate with intense non-thermal radiation/heating.

Clear Lorentz profiles were observed at 70 min, 93 min and 117 min, which
correspond to the impulsive increase of the H$\alpha$ line intensity. On the
basis of the theory mentioned above, pressure broadening implies that the
collision between electrons and hydrogen atoms becomes more frequent, that is
to say, electron density increases significantly. It is well known that in the
impulsive phase of solar flares, energetic electrons are accelerated in the
corona, then propagate downward along flare loops and finally collide with
hydrogen atoms in the (deep) chromosphere. Meanwhile, as revealed in (Namekata
et al., 2020), if a strong non-thermal heating occurs in this stellar flare,
the forming height of H${\alpha}$ line can further extend into the lower
chromosphere and then a larger opacity (or strong self-absorptions) would also
contribute to a wider line broadening.

After 117 min, the H$\alpha$ line intensity stops increasing, indicating the
end of the impulsive phase. From Fig. 4 we can see that the HWHM of the
Gaussian component in the Voigt-Gauss fitting function becomes non-zero again,
meaning that the spectral line profiles start to deviate from a Lorentz
function and the electron density begins to decrease , but the spectral line
profiles are still wider than their pre-flare state. This likely indicates
that a non-thermal heating, i.e., energetic electron beams, is still
continuing but with a softer/weaker energy injection. Indeed, Namekata et al.
(2020) recently reported a good correspondence between the white-light
emission and the H${\alpha}$ line width variation in stellar flares of the
M-type star AD Leo, suggesting H${\alpha}$ line broadening has a fairly good
relationship with non-thermal heating. In their observed flares, a similar
decease in the H${\alpha}$ line width from the impulsive to decay phase was
also found. Such a scenario is in accord with the physical explanation of the
flare decaying phase in the standard solar flare scenario (e.g., Benz & Güdel,
2010; Shibata & Magara, 2011).

###  3.4 Redward asymmetry of the H$\alpha$ line profile

On the Sun, numerous spectroscopic observations of flares have revealed that
chromospheric lines often exhibit line asymmetries during flares. Redward
asymmetries indicative of velocities of several tens of km s-1 have been
frequently reported during the flare impulsive phase. This is commonly thought
to be indued by chromospheric condensation (e.g., Fisher et al., 1985), the
downward movement of the chromospheric plasma caused by the high pressure in
the chromosphere. Downflows of dense plasma draining back to the solar surface
along post-flare loops also lead to redward asymmetries of chromospheric
spectral lines (e.g., Tian et al., 2015). Blue asymmetries or blue shifts of
spectral lines have also been detected in the early phase of solar flares,
which indicate upflows with speeds of a few to a few hundreds km s-1 in
different spectral lines. This has often been attributed to the chromospheric
evaporation, namely warm/hot plasma upflows driven by intense chromospheric
heating. In the meantime, CMEs may also produce absorption line asymmetries if
erupting filament materials are contained in CMEs. For example, Den &
Kornienko (1993) reported a strong mass ejection during a solar flare through
obvious blue wing absorption in H${\alpha}$, indicating a radial velocity of
over 600 km s-1. Similarly, Ding et al. (2003) investigated a series of
filament-induced absorption in the H${\alpha}$ blue wing, revealing a good
temporal correspondence between the acceleration process of a filament and the
impulsive phase of a two-ribbon solar flare.

During flares on other stars, Doppler shifts or asymmetries of spectral lines
have also been frequently detected (e.g., Houdebine et al., 1990; Gunn et al.,
1994; Crespo-Chacón et al., 2006; Vida et al., 2016; Fuhrmeister et al., 2018;
Cao et al., 2019). Using X-ray observations, Argiroffi et al. (2019) recently
reported a long-lasting blueshift of $\sim$ 90 km s-1 of cool plasma (3 MK) at
the decay phase of a large stellar flare on an evolved giant star HR 9024.
They explained this cool upflow as a CME candidate. At optical wavelengths,
spectral line asymmetries are more frequently reported (e.g., Leitzinger et
al., 2014; Honda et al., 2018; Koller et al., 2021), although their
explantations sometimes are not straightforward. For instance, Honda et al.
(2018) found that the H${\alpha}$ line in a stellar flare on the dMe star EV
Lac exhibited a long-lived blueward asymmetry, with a velocity of $\sim$ 100
km s-1, even during the whole duration of the flare. They speculated that the
blue asymmetries might be caused by possible filament activation/eruption and
evaporation, but little is known about the definite origin of such long-
duration blue asymmetries yet. Koller et al. (2021) have searched asymmetries
of the Balmer lines in time-resolved spectroscopy provided by the Sloan
Digital Sky Survey data release 14 and reported six CME candidates with excess
flux in line wings, with the highest projected velocities of 300$-$700 km s-1.
In fact, most (5/6) of these asymmetries appear at the red wings. Koller et
al. (2021) suggested that these red-wing asymmetries may be caused by several
possible processes, i.e., backward-directed CMEs occurring close to the limb,
falling cool material during the eruption, or severe chromospheric
condensation.

In our case, apart from obvious line broadening, red-wing asymmetries also
appear during the flare process. From Figure 5, it is clear that the redward
asymmetry becomes more evident after the impulsive increase in the H${\alpha}$
intensity. It strongly suggests that such line asymmetries are a flare-induced
phenomenon. The line asymmetries we detected last towards the end of our
observations, when the H${\alpha}$ intensity starts to decrease. This is timed
to coincide with the line broadening after the flare peak, namely the
continuing non-thermal heating. Moreover, the line asymmetries demonstrate an
increasing amplitude until the end of our observations, with a redshifted
velocity at 100$-$200 km s-1. According to the current knowledge of solar
flares, the long-lived red wing enhancement may originate from three possible
scenarios: (1) a flare-driven coronal rain process. The typical downward
velocity of solar coronal rain is more than 100 km s-1; as revealed in many
solar flare observations, coronal rain are resulted from the radiative cool
process of hot flaring plasma that trapped in post-flare loops, thus they
often become more prominent in the decay phase of flares. These observed
features are basically consistent with the properties of red-wing asymmetries
we detected here. (2) a flare-induced chromospheric condensation process. In
solar flares, chromospheric condensation results from flare non-thermal
heating processes and it usually happens in the impulsive phase, with a
typical downflow velocity of tens of km s-1 (e.g., Ding et al., 1995; Tian et
al., 2015; Yu et al., 2020) but sometimes even up to more than 100 km s-1 as
well (e.g., Ichimoto & Kurokawa, 1984; Tei et al., 2018). Here, the red-wing
asymmetry is accompanied by obvious line broadening after the flare peak, thus
suggesting a continuing non-thermal heating might still be responsible for the
appearance of condensation downflows in the decay phase. (3) a flare-
associated filament/prominence eruption either with a non-radial backward
propagation or with strong magnetic suppression. In case of a backward ejected
filament, the flare must occur near or at the limb of the stellar disk where
its footpoint-produced H${\alpha}$ emission (especially line broadening) at
least can be visible from the Earth. This special spatial location requirement
suggests that a backward filament eruption at least is a very quite rare case.
Another situation is a failed filament eruption, which is first reported on
the Sun (Ji et al., 2003) and recently supposed to be common on M-type stars
due to their strong magnetic suppression (Alvarado-Gómez et al., 2018; Li et
al., 2020). In this case, the ejected filament material may be quickly
confined by overlying magnetic field after the flare peak and then demonstrate
more and more prominent downward mass falling to the stellar disk, which thus
might cause the increasing amplitude of red-wing asymmetries in our
observation. Meanwhile, the ejecting filament at the flare onset is expected
to cause a short-lived blueshift before we detected the red-wing asymmetries.
But this blueshift was absent in the LAMOST observation, possibly due to its
poor temporal resolution (around 23 mins). Caution that only using
chromospheric spectral observations, our current observation is very difficult
to distinguish between these three possible scenarios or their combinations,
since the observed H${\alpha}$ radiation from the star could be a
superposition of flare, post-flare loop, and prominence eruption. In future,
to better understand these physical processes in stellar flares, new time-
resolved spectroscopic observations using multi-temperature lines are needed.

The speed of the moving plasma is 100 km s-1 to 200 km s-1, far below the
star’s escape velocity. However, we should remember that this speed only
represents the lower limit of the true velocity of the moving material.
Because the Doppler shift can only be used to measure the velocity component
along the line of sight, and ejected filament material originating from the
limb very likely has a small velocity component in the line of sight. Thus, we
can not exclude the possibility that the backward filament eruption might
drive a CME if the red-wing asymmetries were to be an eruption. Moreover, as
shown in Figure 4 (e), the moving plasma responsible for the redward asymmetry
likely decelerates. The average deceleration rate is computed as 7.6 m s-2,
which is much smaller than the gravity of the star. In the filament/prominence
eruption scenario, this very small deceleration rate might supports that the
ejecting filament material is flying along a non-radial direction that is near
perpendicularly to the line of sight, suffering a severe projection effect.

On the other hand, the acceleration of plasma during solar coronal rain is
about 1/3 $g_{\rm eff}$ (Antolin, 2020). If we assume that the falling
materials move downward with a constance acceleration of 1/3 $g_{\rm eff}$ and
the propagation distance scales with the radius of the star, the maximum
velocity of the falling plasma on a star can be estimated as follows,

$v=v_{\sun}\sqrt{\frac{gR}{g_{\sun}R_{\sun}}},$ (8)

where $v_{\sun}$ is the maximum velocity of plasma during coronal rain on the
Sun, which is about 100 km s-1 according to Antolin (2020). $g$ and $g_{\sun}$
are the gravitational accelerations on the star and the Sun, respectively. And
the stellar and solar radii are represented by $R$ and $R_{sun}$,
respectively. Our calculation yielded a value of 126 km s-1 for $v$, which is
consistent with the velocities inferred from our RB asymmetry analysis,
supporting the possible coronal rain scenario mentioned above.

##  4 Conclusion

From the LAMOST DR7 data archive, we identified a strong flare event on an
M4-type star. The key findings are summarized as follows,

(1) The flare shows an impulsive increase followed by a gradual decrease in
the H$\alpha$ line intensity. The impulsive increase starts at 70 min and the
H$\alpha$ intensity peaks at 117 min since the beginning of the observation.
The maximum radiation power of the flare through the H$\alpha$ line is
$1.43\times 10^{29}$ erg s-1. During the observation, the flare radiated
energy in H$\alpha$ is $4.45\times 10^{32}$ erg, and the total flare energy
radiated through H$\alpha$ is estimated to be on the order of $10^{33}$ erg.

(2) During the impulsive phase, the H$\alpha$ line appears to have a Lorentz
profile, which we think is caused by downward propagating nonthermal electrons
colliding with hydrogen atoms in the chromosphere. After t=117 min, the
spectral profile gradually deviates from a Lorentz profile due to the decrease
of non-thermal electron acceleration and the resultant decrease of electron
density in the chromosphere.

(3) H$\alpha$ line asymmetries corresponding to plasma moving away from the
Earth at a velocity of 100$-$200 km s-1 are detected. The asymmetries might
result from three possible scenarios: (1) flare-driven coronal rain, (2)
chromospheric condensation, or (3) a filament/prominence eruption either with
a non-radial backward propagation or with strong magnetic suppression. The
total mass of the moving plasma is estimated to be on the order of $10^{15}$
kg.

##  5 Acknowledgements

H.H. is supported by the National Key R&D Program of China (2019YFA0405000).
This work is also supported by NSFC grants 11825301, 11790304, & 11973059, the
fellowship of China National Postdoctoral Program for Innovative Talents
(BX20200013), and Peking University Undergraduate Research Program. We thank
the referee for the very helpful comments and constructive suggestions, and
thank Y.-J. Zhu, Dr. Y. Li, and Dr. H. Li for helpful discussions. The Large
Sky Area Multi-object Fiber Spectroscopic Telescope, LAMOST (Guoshoujing
Telescope) is a National Major Scientific Project built by the Chinese Academy
of Sciences. Funding for the project has been provided by the National
Development and Reform Commission. LAMOST is operated and managed by the
National Astronomical Observatories, Chinese Academy of Sciences.

## References

  * Aarnio et al. (2011) Aarnio, A. N., Stassun, K. G., Hughes, W. J., et al. 2011, Sol. Phys., 268, 195. doi:10.1007/s11207-010-9672-7
  * Antolin (2020) Antolin, P. 2020, Plasma Physics and Controlled Fusion, 62, 014016. doi:10.1088/1361-6587/ab5406
  * Alvarado-Gómez et al. (2018) Alvarado-Gómez, J. D., Drake, J. J., Cohen, O., et al. 2018, ApJ, 862, 93. doi:10.3847/1538-4357/aacb7f
  * Alvarado-Gómez et al. (2021) Alvarado-Gómez, J. D., Drake, J. J., Cohen, O., et al. 2021, arXiv:2111.09704
  * Argiroffi et al. (2019) Argiroffi, C., Reale, F., Drake, J. J., et al. 2019, Nature Astronomy, 3, 742. doi:10.1038/s41550-019-0781-4
  * Benz & Güdel (2010) Benz, A. O. & Güdel, M. 2010, ARA&A, 48, 241. doi:10.1146/annurev-astro-082708-101757
  * Bond et al. (2001) Bond, H. E., Mullan, D. J., O’Brien, M. S., et al. 2001, ApJ, 560, 919. doi:10.1086/322980
  * Butler et al. (1988) Butler, C. J., Rodono, M., & Foing, B. H. 1988, A&A, 206, L1
  * Canfield et al. (1990) Canfield, R. C., Penn, M. J., Wulser, J.-P., et al. 1990, ApJ, 363, 318. doi:10.1086/169345
  * Carmichael (1964) Carmichael, H. 1964, NASA Special Publication, 451
  * Cao et al. (2019) Cao, D., Gu, S., Ge, J., et al. 2019, MNRAS, 482, 988. doi:10.1093/mnras/sty2768
  * Chen (2011) Chen, P. F. 2011, Living Reviews in Solar Physics, 8, 1. doi:10.12942/lrsp-2011-1
  * Chen et al. (2021) Chen, H., Tian, H., Li, L., et al. 2021, arXiv:2112.01354
  * Cox & Pilachowski (2000) Cox, A. N. & Pilachowski, C. A. 2000, Physics Today, 53, 77. doi:10.1063/1.1325201
  * Crespo-Chacón et al. (2006) Crespo-Chacón, I., Montes, D., García-Alvarez, D., et al. 2006, A&A, 452, 987. doi:10.1051/0004-6361:20053615
  * Cui et al. (2012) Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, Research in Astronomy and Astrophysics, 12, 1197. doi:10.1088/1674-4527/12/9/003
  * De Pontieu et al. (2009) De Pontieu, B., McIntosh, S. W., Hansteen, V. H., et al. 2009, ApJ, 701, L1. doi:10.1088/0004-637X/701/1/L1
  * Den & Kornienko (1993) Den, O. E. & Kornienko, G. I. 1993, Astronomy Reports, 37, 76
  * Ding et al. (2003) Ding, M. D., Chen, Q. R., Li, J. P., et al. 2003, ApJ, 598, 683. doi:10.1086/378877
  * Ding et al. (1995) Ding, M. D., Fang, C., & Huang, Y. R. 1995, Sol. Phys., 158, 81. doi:10.1007/BF00680836
  * Doyle et al. (2020) Doyle, L., Ramsay, G., & Doyle, J. G. 2020, MNRAS, 494, 3596. doi:10.1093/mnras/staa923
  * Falcon et al. (2015) Falcon, R. E., Rochau, G. A., Bailey, J. E., et al. 2015, ApJ, 806, 214. doi:10.1088/0004-637X/806/2/214
  * Fisher et al. (1985) Fisher, G. H., Canfield, R. C., & McClymont, A. N. 1985, ApJ, 289, 425. doi:10.1086/162902
  * Fletcher et al. (2011) Fletcher, L., Dennis, B. R., Hudson, H. S., et al. 2011, Space Sci. Rev., 159, 19. doi:10.1007/s11214-010-9701-8
  * Fuhrmeister & Schmitt (2004) Fuhrmeister, B. & Schmitt, J. H. M. M. 2004, A&A, 420, 1079. doi:10.1051/0004-6361:20035644
  * Fuhrmeister et al. (2018) Fuhrmeister, B., Czesla, S., Schmitt, J. H. M. M., et al. 2018, A&A, 615, A14. doi:10.1051/0004-6361/201732204
  * Gan et al. (1993) Gan, W. Q., Rieger, E., & Fang, C. 1993, ApJ, 416, 886. doi:10.1086/173285
  * Ginsburg & Mirocha (2011) Ginsburg, A. & Mirocha, J. 2011, Astrophysics Source Code Library. ascl:1109.001
  * Gosling (1993) Gosling, J. T. 1993, J. Geophys. Res., 98, 18937. doi:10.1029/93JA01896
  * Guenther & Emerson (1997) Guenther, E. W. & Emerson, J. P. 1997, A&A, 321, 803
  * Gunn et al. (1994) Gunn, A. G., Doyle, J. G., Mathioudakis, M., et al. 1994, A&A, 285, 489
  * Haisch et al. (1991) Haisch, B., Strong, K. T., & Rodono, M. 1991, ARA&A, 29, 275. doi:10.1146/annurev.aa.29.090191.001423
  * He et al. (2015) He, H., Wang, H., & Yun, D. 2015, ApJS, 221, 18. doi:10.1088/0067-0049/221/1/18
  * He et al. (2018) He, H., Wang, H., Zhang, M., et al. 2018, ApJS, 236, 7. doi:10.3847/1538-4365/aab779
  * Honda et al. (2018) Honda, S., Notsu, Y., Namekata, K., et al. 2018, PASJ, 70, 62. doi:10.1093/pasj/psy055
  * Hirayama (1974) Hirayama, T. 1974, Sol. Phys., 34, 323. doi:10.1007/BF00153671
  * Holtsmark (1919) Holtsmark, J. 1919, Annalen der Physik, 363, 577. doi:10.1002/andp.19193630702
  * Houdebine et al. (1990) Houdebine, E. R., Foing, B. H., & Rodono, M. 1990, A&A, 238, 249
  * Howell et al. (2014) Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398. doi:10.1086/676406
  * Ichimoto & Kurokawa (1984) Ichimoto, K. & Kurokawa, H. 1984, Sol. Phys., 93, 105. doi:10.1007/BF00156656
  * Illing & Hundhausen (1985) Illing, R. M. E. & Hundhausen, A. J. 1985, J. Geophys. Res., 90, 275. doi:10.1029/JA090iA01p00275
  * Ji et al. (2003) Ji, H., Wang, H., Schmahl, E. J., et al. 2003, ApJ, 595, L135. doi:10.1086/378178
  * Kanodia et al. (2021) Kanodia, S., Ramsey, L. W., Maney, M., et al. 2021, arXiv:2111.14647
  * Khodachenko et al. (2007) Khodachenko, M. L., Ribas, I., Lammer, H., et al. 2007, Astrobiology, 7, 167. doi:10.1089/ast.2006.0127
  * Koller et al. (2021) Koller, F., Leitzinger, M., Temmer, M., et al. 2021, A&A, 646, A34. doi:10.1051/0004-6361/202039003
  * Kopp & Holzer (1976) Kopp, R. A. & Holzer, T. E. 1976, Sol. Phys., 49, 43. doi:10.1007/BF00221484
  * Kowalski et al. (2017) Kowalski, A. F., Allred, J. C., Uitenbroek, H., et al. 2017, ApJ, 837, 125. doi:10.3847/1538-4357/aa603e
  * Leitzinger et al. (2014) Leitzinger, M., Odert, P., Greimel, R., et al. 2014, MNRAS, 443, 898. doi:10.1093/mnras/stu1161
  * Luque et al. (2003) Luque, J. M., Calzada, M. D., & Sáez, M. 2003, Journal of Physics B Atomic Molecular Physics, 36, 1573. doi:10.1088/0953-4075/36/8/311
  * Li & Ding (2011) Li, Y. & Ding, M. D. 2011, ApJ, 727, 98. doi:10.1088/0004-637X/727/2/98
  * Li et al. (2020) Li, T., Hou, Y., Yang, S., et al. 2020, ApJ, 900, 128. doi:10.3847/1538-4357/aba6ef
  * Li et al. (2021) Li, D., Hong, Z., & Ning, Z. 2021, arXiv:2112.06118
  * Li et al. (2020) Li, L., Peter, H., Chitta, L. P., et al. 2020, ApJ, 905, 26. doi:10.3847/1538-4357/abc68c
  * Lin et al. (2003) Lin, J., Soon, W., & Baliunas, S. L. 2003, New A Rev., 47, 53. doi:10.1016/S1387-6473(02)00271-3
  * Maehara et al. (2021) Maehara, H., Notsu, Y., Namekata, K., et al. 2021, PASJ, 73, 44. doi:10.1093/pasj/psaa098
  * Maehara et al. (2021) Maehara, H., Notsu, Y., Namekata, K., et al. 2021, The 20.5th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun (CS20.5), 142. doi:10.5281/zenodo.4564272
  * Maehara et al. (2012) Maehara, H., Shibayama, T., Notsu, S., et al. 2012, Nature, 485, 478. doi:10.1038/nature11063
  * Martínez-Sykora et al. (2011) Martínez-Sykora, J., De Pontieu, B., Hansteen, V., et al. 2011, ApJ, 732, 84. doi:10.1088/0004-637X/732/2/84
  * Moschou et al. (2019) Moschou, S.-P., Drake, J. J., Cohen, O., et al. 2019, ApJ, 877, 105. doi:10.3847/1538-4357/ab1b37
  * Muheki et al. (2020) Muheki, P., Guenther, E. W., Mutabazi, T., et al. 2020, A&A, 637, A13. doi:10.1051/0004-6361/201936904
  * Muheki et al. (2020) Muheki, P., Guenther, E. W., Mutabazi, T., et al. 2020, MNRAS, 499, 5047. doi:10.1093/mnras/staa3152
  * Namekata et al. (2017) Namekata, K., Sakaue, T., Watanabe, K., et al. 2017, ApJ, 851, 91. doi:10.3847/1538-4357/aa9b34
  * Namekata et al. (2020) Namekata, K., Maehara, H., Sasaki, R., et al. 2020, PASJ, 72, 68. doi:10.1093/pasj/psaa051
  * Namekata et al. (2022) Namekata, K., Maehara, H., Honda, S., et al. 2022, Nature Astronomy. doi.org/10.1038/s41550-021-01532-8
  * Notsu et al. (2013) Notsu, Y., Shibayama, T., Maehara, H., et al. 2013, ApJ, 771, 127. doi:10.1088/0004-637X/771/2/127
  * Osten & Wolk (2017) Osten, R. A. & Wolk, S. J. 2017, Living Around Active Stars, 328, 243. doi:10.1017/S1743921317004252
  * Samanta et al. (2021) Samanta, T., Tian, H., Chen, B., et al. 2021, The Innovation, 2, 100083. doi:10.1016/j.xinn.2021.100083
  * Schad et al. (2021) Schad, T. A., Dima, G. I., & Anan, T. 2021, ApJ, 916, 5. doi:10.3847/1538-4357/ac01eb
  * Schaefer et al. (2000) Schaefer, B. E., King, J. R., & Deliyannis, C. P. 2000, ApJ, 529, 1026. doi:10.1086/308325
  * Schmieder et al. (1987) Schmieder, B., Forbes, T. G., Malherbe, J. M., et al. 1987, ApJ, 317, 956. doi:10.1086/165344
  * Shibata & Magara (2011) Shibata, K. & Magara, T. 2011, Living Reviews in Solar Physics, 8, 6. doi:10.12942/lrsp-2011-6
  * Shibata et al. (2013) Shibata, K., Isobe, H., Hillier, A., et al. 2013, PASJ, 65, 49. doi:10.1093/pasj/65.3.49
  * Sturrock (1966) Sturrock, P. A. 1966, Nature, 211, 695. doi:10.1038/211695a0
  * Su et al. (1998) Su, D. Q., Cui, X., Wang, Y., et al. 1998, Proc. SPIE, 3352, 76. doi:10.1117/12.319250
  * Tei et al. (2018) Tei, A., Sakaue, T., Okamoto, T. J., et al. 2018, PASJ, 70, 100. doi:10.1093/pasj/psy047
  * Tremblay & Bergeron (2009) Tremblay, P.-E. & Bergeron, P. 2009, ApJ, 696, 1755. doi:10.1088/0004-637X/696/2/1755
  * Tian et al. (2011) Tian, H., McIntosh, S. W., De Pontieu, B., et al. 2011, ApJ, 738, 18. doi:10.1088/0004-637X/738/1/18
  * Tian et al. (2012) Tian, H., McIntosh, S. W., Xia, L., et al. 2012, ApJ, 748, 106. doi:10.1088/0004-637X/748/2/106
  * Tian et al. (2015) Tian, H., Young, P. R., Reeves, K. K., et al. 2015, ApJ, 811, 139. doi:10.1088/0004-637X/811/2/139
  * Tu et al. (2020) Tu, Z.-L., Yang, M., Zhang, Z. J., et al. 2020, ApJ, 890, 46. doi:10.3847/1538-4357/ab6606
  * Tu et al. (2021) Tu, Z.-L., Yang, M., Wang, H.-F., et al. 2021, ApJS, 253, 35. doi:10.3847/1538-4365/abda3c
  * Vida et al. (2016) Vida, K., Kriskovics, L., Oláh, K., et al. 2016, A&A, 590, A11. doi:10.1051/0004-6361/201527925
  * Vida et al. (2019) Vida, K., Leitzinger, M., Kriskovics, L., et al. 2019, A&A, 623, A49. doi:10.1051/0004-6361/201834264
  * Wang et al. (2021) Wang, J., Xin, L. P., Li, H. L., et al. 2021, ApJ, 916, 92. doi:10.3847/1538-4357/ac096f
  * Weisskopf (1932) Weisskopf, V. 1932, Zeitschrift fur Physik, 75, 287. doi:10.1007/BF01341858
  * Wiese & Fuhr (2009) Wiese, W. L. & Fuhr, J. R. 2009, Journal of Physical and Chemical Reference Data, 38, 1129. doi:10.1063/1.3254213
  * Yan et al. (2021) Yan, Y., He, H., Li, C., et al. 2021, MNRAS, 505, L79. doi:10.1093/mnrasl/slab055
  * Yashiro & Gopalswamy (2009) Yashiro, S. & Gopalswamy, N. 2009, Universal Heliophysical Processes, 257, 233. doi:10.1017/S1743921309029342
  * Yelle et al. (2008) Yelle, R., Lammer, H., & Ip, W.-H. 2008, Space Sci. Rev., 139, 437. doi:10.1007/s11214-008-9420-6
  * Yu et al. (2020) Yu, K., Li, Y., Ding, M. D., et al. 2020, ApJ, 896, 154. doi:10.3847/1538-4357/ab9014
  * Zarro et al. (1988) Zarro, D. M., Canfield, R. C., Strong, K. T., et al. 1988, ApJ, 324, 582. doi:10.1086/165919
  * Zhang & Low (2005) Zhang, M. & Low, B. C. 2005, ARA&A, 43, 103. doi:10.1146/annurev.astro.43.072103.150602
  * Zhang et al. (2016) Zhang, Q. M., Li, D., & Ning, Z. J. 2016, ApJ, 832, 65. doi:10.3847/0004-637X/832/1/65
  * Zhu et al. (2019) Zhu, Y., Kowalski, A. F., Tian, H., et al. 2019, ApJ, 879, 19. doi:10.3847/1538-4357/ab2238
  * Švestka (1972) Švestka, Z. 1972, ARA&A, 10, 1. doi:10.1146/annurev.aa.10.090172.000245


# Extended Capture Point and Optimization-based Control for Quadrupedal Robot
Walking on Dynamic Rigid Surfaces

Amir Iqbal and Yan Gu Amir Iqbal and Yan Gu are with the Department of
Mechanical Engineering, University of Massachusetts Lowell, Lowell, MA 01854,
U.S.A. Emails<EMAIL_ADDRESS>and<EMAIL_ADDRESS>Research by A.
Iqbal and Y. Gu was supported in part by NSF under Grant no. CMMI-1934280.
Corresponding author: Y. Gu.

###### Abstract

Stabilizing legged robot locomotion on a dynamic rigid surface (DRS) (i.e.,
rigid surface that moves in the inertial frame) is a complex planning and
control problem. The complexity arises due to the hybrid nonlinear walking
dynamics subject to explicitly time-varying holonomic constraints caused by
the surface movement. The first main contribution of this study is the
extension of the capture point from walking on a static surface to locomotion
on a DRS as well as the use of the resulting capture point for online motion
planning. The second main contribution is a quadratic-programming (QP) based
feedback controller design that explicitly considers the DRS movement. The
stability and robustness of the proposed control approach are validated
through simulations of a quadrupedal robot walking on a DRS with a rocking
motion. The simulation results also demonstrate the improved walking
performance compared with our previous approach based on offline planning and
input-output linearizing control that does not explicitly guarantee the
feasibility of ground contact constraints.

##  I Introduction

Legged robots, mimicking nature’s design of terrestrial animals, have many
potential advantages compared to wheeled or tracked robots in traversing
unstructured terrains. Various control approaches have been created to enable
stable and robust locomotion on stationary surfaces [1, 2]. Yet, sustaining
legged locomotion on a nonstationary surface has not been fully addressed.

The objective of this study is to derive a control approach that stabilizes
locomotion on a dynamic rigid surface (DRS) by explicitly addressing the time-
varying surface movement and the feasibility of ground contact forces. A key
element of the proposed approach is an online footstep planner that considers
the surface motion through the extension of the concept of capture point [3,
4] from stationary to dynamic surfaces. Another key element is a quadratic
programming (QP) based control method that incorporates the surface motion.
Related work is reviewed next.

###  I-A Related work on Legged Locomotion Control for Rigid Stationary and
Deformable Surfaces

Previous control approaches for legged robot locomotion mainly focus on static
surfaces (e.g., pavement and stairs) [5, 6, 7, 8, 9, 10, 1, 2]. Yet, they may
not work well for locomotion on a deformable surface or a DRS because they do
not explicitly account for the perturbations caused by the surface deformation
or movement.

One of the pioneering efforts in addressing locomotion on deformable surfaces
is the modeling of the dynamic interaction between the robot’s legs and
granular terrains (e.g., sand and snow) [11]. Based on this model, Xiong et
al. [12] have developed a stability region criterion for granular terrain
walking to guide controller design. However, it is unclear how effective these
methods would be for sustaining locomotion on a DRS.

###  I-B Related Work on DRS Locomotion Control

Stabilizing legged locomotion on a DRS is fundamentally challenging due to the
high complexity of the associated robot dynamics. First, the robot dynamics
are inherently hybrid involving state-triggered discrete behaviors (e.g., foot
touchdowns) [13]. Second, if not properly addressed, the time-varying
perturbation at a foot-surface contact point caused by the surface motion can
destabilize the robot.

Recently, planning and control of DRS locomotion have been increasingly
extensively studied. For surfaces with a motion affected by the locomotor
(e.g., passive surfaces with relatively small inertia), planning and control
methods for bipedal walking on different-sized balls have been proposed based
on reduced-order [14] and full-order robot models [15]. Controller design for
legged locomotion on a floating island has also been studied based on an
eight-legged rimless wheel robot model [16]. Still, these methods cannot be
directly used to handle locomotion on a DRS whose motion is not affected by
the robot, due to the differences in the associated robot dynamics.

To address locomotion on a DRS whose motion is not affected by the robot
(e.g., vessels and aircraft), our previous work [17] has analyzed the effects
of surface movement on a robot’s locomotion stability based on the reduced-
order spring-loaded inverted pendulum (SLIP) robot model. We have also
recently derived and experimentally evaluated a provably stabilizing
controller for enabling quadrupedal locomotion on a DRS [18]. Yet, this
controller has two major limitations towards effective implementation in real-
world scenarios. First, the previous planning algorithm is not suitable for
online motion generation due to the heavy computational load caused by its
underlying full-order model. Second, the controller does not explicitly
account for the feasibility of ground-contact constraints during the actual
walking process.

###  I-C Related Work on Capture Point and QP-based Control

A capture point is defined as the point that an inverted pendulum model needs
to instantaneously extend its leg to in order to come to a complete stop
within the new step [3, 4]. Due to its computational efficiency, the concept
of capture point has been used to plan footsteps for legged locomotion on
static [19] and uneven terrains [20, 21] as well as for disturbance rejection
[22]. However, to the best of our knowledge, the capture point has not been
extended to solve the footstep planning problem for DRS locomotion. Due to its
capability of incorporating feasibility constraints and its low computational
load suitable for real-time implementation, QP-based control design has been
extensively studied for legged locomotion over high-slope terrains [23] and
unstructured terrains [2]. Recently, a QP-based controller utilizing a full-
order robot model has shown effectiveness in asymptotically tracking desired
trajectories [24]. Still, these controllers mainly focus on static surface
locomotion and thus do not explicitly consider the time-varying motion of a
DRS.

###  I-D Contributions

This paper aims to derive a planning and control approach that enables stable
and robust quadrupedal walking on a DRS by extending our previous work [18]
and solving its two major limitations, which are: a) lack of feasibility
guarantees in the controller design and b) heavy computational load of motion
planning that prevents online implementation.

The main contributions of this paper are: a) Extending the concept of capture
point to DRS locomotion by explicitly considering the surface motion, and
designing an online footstep planner based on the extended capture point. b)
Formulating a QP-based controller that ensures the feasibility of contact
forces and explicitly incorporates the surface motion in a time-varying
holonomic constraint. c) Validating the stability and robustness of the
proposed approach through simulations of quadrupedal walking on a DRS with a
rocking motion. d) Providing comparative simulation results to demonstrate the
improved robustness of the proposed method over our previous work.

##  II HYBRID FULL-ORDER ROBOT MODEL

This section presents the full-order model of a quadrupedal robot that walks
on a DRS, which captures the dynamic behaviors of the robot’s full degrees of
freedom (DOFs). The model serves as a basis for the proposed QP-based
controller design.

A walking robot’s dynamics are hybrid because walking naturally involves
continuous-time motion (e.g., foot swinging) and discrete events (e.g., foot
touchdowns). During locomotion on a DRS, the robot’s dynamics are also
affected the time-varying movement of the surface. The hybrid, time-varying
robot dynamics are explained next.

Continuous-phase dynamics. Let $\mathbf{q}$ be the robot’s generalized
coordinates, which comprise the base pose (i.e., position $\mathbf{p}_{b}$ and
orientation $\bm{\gamma}_{b}$) and joint angles
$\begin{bmatrix}q_{1}~{}q_{2}~{}...~{}q_{n}\end{bmatrix}^{T}$; that is,
$\mathbf{q}:=\begin{bmatrix}\mathbf{p}_{b}^{T},~{}\bm{\gamma}_{b}^{T},~{}q_{1},~{}q_{2},~{}...,~{}q_{n}\end{bmatrix}^{T}.$
Let $\mathcal{Q}\subset\mathbb{R}^{n+6}$ be the robot’s configuration space.
Then $\mathbf{q}\in\mathcal{Q}$. Let $\mathbf{u}\in U\subset\mathbb{R}^{m}$ be
the robot’s joint torques with $U$ the admissible set of joint torques.

Using Lagrange’s method, the continuous-phase robot dynamics can be obtained
as:

$\mathbf{M(q)\ddot{q}+C(q,\dot{q}):=Bu}+\mathbf{J}_{f}^{T}(\mathbf{q})\mathbf{F}_{f},$
(1)

where
$\mathbf{M}(\mathbf{q}):\mathcal{Q}\rightarrow\mathbb{R}^{(n+6)\times(n+6)}$
is the inertia matrix,
$\mathbf{\mathbf{C}(\mathbf{q},\dot{\mathbf{q}})}:\mathcal{TQ}\rightarrow\mathbb{R}^{(n+6)}$
represents the sum of centrifugal, Coriolis, and gravitational terms with
$\mathcal{TQ}$ the tangent space of $\mathcal{Q}$, and
$\mathbf{B}\in\mathbb{R}^{(n+6)\times m}$ is the actuator selection matrix.
The matrix $\mathbf{J}_{f}(\mathbf{q})$ is the contact Jacobian whose
expression varies during different phases because of different feet in
support. The vector $\mathbf{F}_{f}(\mathbf{q})\in\mathbb{R}^{n_{ct}}$
represents the external force acting at the support feet. During three-
dimensional ($3$-D) quadrupedal walking without slip, three of the four legs
are in secured contact with the surface, and therefore $n_{ct}=9$.

During robot walking on a DRS whose motion is not affected by the robot (e.g.,
elevators and vessels), the effects of the surface motion on the robot
dynamics can be modeled as a time-varying holonomic constraint.

Let $\mathbf{r}_{p}(t)\in\mathbb{R}^{n_{ct}}$ be the positions of the surface
at the support-foot locations. Note that for a DRS, the surface positions and
motions at the foot-contact points (i.e., $\mathbf{r}_{p}(t)$ and its
derivatives) are explicitly time-varying. Let
$\mathbf{r}_{f}(\mathbf{q})\in\mathbb{R}^{n_{ct}}$ be the positions of the
support feet. The contact Jacobian matrix, by definition, is expressed as
$\mathbf{J}_{f}(\mathbf{q}):=\frac{\partial\mathbf{r}_{f}}{\partial\mathbf{q}}(\mathbf{q})$.

When the support feet do not slip on the surface, the time-varying holonomic
constraint can be formed as:

$\mathbf{J}_{f}(\mathbf{q})\ddot{\mathbf{q}}+\mathbf{\dot{J}}_{f}((\mathbf{q}),\dot{\mathbf{q}})\dot{\mathbf{q}}=\ddot{\mathbf{r}}_{p}(t).$
(2)

Switching surface When a swing foot strikes the walking surface, it becomes a
new stance foot while another foot begins to swing in the air and acts as the
new swing foot. The foot landing event can be defined through the switching
surface,
$S:=\\{(\mathbf{q},\dot{\mathbf{q}},\mathbf{r}_{p},\dot{\mathbf{r}}_{p})\in\mathcal{TQ}\times\mathbb{R}^{n_{ct}}\times\mathbb{R}^{n_{ct}}:\phi(\mathbf{q},\mathbf{r}_{p})=0,~{}\dot{\phi}(\mathbf{q},\dot{\mathbf{q}},\mathbf{r}_{p},\dot{\mathbf{r}}_{p})<0\\},$
where the scalar function $\phi$ is the distance between the swing foot and
the surface.

##  III Capture-point based Online Planning

This section introduces the extension of the capture point from stationary to
DRS walking, as well as the proposed online planning method based on the
extended capture point. The desired trajectories produced by the planner will
be tracked by the proposed QP-based control as introduced in Sec. IV.

###  III-A Capture Point on a DRS

The concept of capture point for a static surface has been extensively studied
[3, 4]. In this work, we extend the capture point to a DRS based on a reduced-
order locomotor model, that is, a $2$-D linear inverted pendulum model (LIPM)
with a point mass atop a massless leg (see Fig. 1). The position of the point
mass coincides with the robot’s center of mass (CoM).

Let $\mathbf{r}_{wc}=(x_{wc},~{}x_{wc},~{}z_{wc})$ and
$\mathbf{r}_{ws}=(x_{ws},~{}x_{ws},~{}z_{ws})$ be the absolute positions of
the CoM and the leg’s far end (i.e., point $S$) in the world frame,
respectively. Then, the CoM position relative to point $S$ is given by:

$\mathbf{r}_{sc}=(x_{sc},~{}x_{sc},~{}z_{sc})=\mathbf{r}_{wc}-\mathbf{r}_{ws}.$
(3)

The equation of motion of the LIPM is expressed as:

$\displaystyle\ddot{x}_{wc}=\tfrac{f_{a}x_{sc}}{mr}$
$\displaystyle\sin\theta,~{}\ddot{y}_{wc}=\frac{f_{a}y_{sc}}{mr}\sin\theta,~{}\ddot{z}_{wc}=\frac{f_{a}}{m}\cos\theta-g,$
(4)

where $m$ is the mass of the LIPM, $\theta$ is the angle of the leg relative
to the vertical axis, $g$ is the magnitude of the gravitational acceleration,
and $r$ is the projected length of $\mathbf{r}_{sc}$ on the horizontal plane.
The scalar variable $f_{a}$ is the norm of the axial force that the massless
rod applies to the point mass.

Figure 1: Illustration of the 3-D LIPM model on a DRS. The model is
constrained to maintain constant height $z_{0}$ between the CoM and the point
S on the DRS.

It is assumed in this study that the vertical distance $z_{sc}$ between the
CoM and point $S$ is constant during walking (see Fig. 1); that is,
$z_{sc}=z_{0}$ with $z_{0}$ a positive number. This assumption is analogous to
the simplifying assumption of the LIPM model on a stationary surface that the
point-mass height over the surface is constant [3].

Under the assumption that $z_{sc}$ is constant, we have
$\dot{z}_{wc}=\dot{z}_{ws}$ and $\ddot{z}_{wc}=\ddot{z}_{ws}$. Then, from (4),
the expression of $f_{a}$ becomes
$f_{a}=m\frac{(\ddot{z}_{ws}+g)}{\cos\theta}.$ Accordingly, the equation of
motion of the LIPM on the horizontal plane is now given by:

$\displaystyle\ddot{x}_{wc}=(\ddot{z}_{ws}+g)\frac{x_{sc}}{z_{0}}~{}\text{and}~{}\ddot{y}_{wc}=(\ddot{z}_{ws}+g)\frac{y_{sc}}{z_{0}}.$
(5)

From (3), we obtain $\ddot{x}_{wc}=\ddot{x}_{ws}+\ddot{x}_{sc}$ and
$\ddot{y}_{wc}=\ddot{y}_{ws}+\ddot{y}_{sc}$. Then, using these relations in
(5) yields:

$\displaystyle\ddot{x}_{sc}-\frac{(\ddot{z}_{ws}+g)}{z_{0}}x_{sc}=\ddot{x}_{ws}~{}\text{and}~{}\ddot{y}_{sc}-\frac{(\ddot{z}_{ws}+g)}{z_{0}}y_{sc}=\ddot{y}_{ws}.$
(6)

The ordinary differential equations in (6) are non-homogeneous. If the
horizontal accelerations of point $S$ (i.e., $\ddot{x}_{ws}$ and
$\ddot{y}_{ws}$) are sufficiently small, the forcing terms can be approximated
as zero, and these equations become homogeneous:

$\displaystyle\ddot{x}_{sc}-\frac{(\ddot{z}_{ws}+g)}{z_{0}}x_{sc}=0~{}\text{and}~{}\ddot{y}_{sc}-\frac{(\ddot{z}_{ws}+g)}{z_{0}}y_{sc}=0.$
(7)

Analogous to the derivation of the capture point on a stationary surface [3],
we can define the apparent stiffness of the LIPM as
$k_{a}:=-\frac{(\ddot{z}_{ws}+g)}{z_{0}}$. Then, we reduce the LIPM’s equation
to a spring-mass system with unit mass; that is, $\ddot{x}_{sc}+k_{a}x_{sc}=0$
and $\ddot{y}_{sc}+k_{a}y_{sc}=0$.

The orbital energy of this system associated with the horizontal motion is:

$E_{x}=\frac{1}{2}(\dot{x}_{sc}^{2}+k_{a}x_{sc}^{2})~{}\text{and}~{}E_{y}=\frac{1}{2}(\dot{y}_{sc}^{2}+k_{a}y_{sc}^{2}).$

These energy terms predict whether the point mass will move over the support
point $S$ (i.e., $E_{x}+E_{y}>0$), stop at $S$ (i.e., $E_{x}+E_{y}=0$), or
reverse the direction of motion before reaching $S$ (i.e., $E_{x}+E_{y}<0$).

At the equilibrium (i.e., $E_{x}=0$ and $E_{y}=0$) of the system in (7), the
two eigenvectors are:

$\displaystyle\dot{x}_{sc}=\pm\sqrt{\frac{(\ddot{z}_{ws}+g)}{z_{0}}}x_{sc}~{}\text{and}~{}\dot{y}_{sc}=\pm\sqrt{\frac{(\ddot{z}_{ws}+g)}{z_{0}}}y_{sc}.$
(8)

The stable eigenvector corresponds to the case when $x_{sc}$ and
$\dot{x}_{sc}$ (and $y_{sc}$ and $\dot{y}_{sc}$) have opposite signs, which
lead to the following horizontal coordinates of the instantaneous capture
point of the $3$-D LIPM on a DRS:

$\displaystyle
x_{cap}:=\sqrt{\frac{z_{0}}{(\ddot{z}_{ws}+g)}}\dot{x}_{sc}~{}\text{and}~{}y_{cap}:=\sqrt{\frac{z_{0}}{(\ddot{z}_{ws}+g)}}\dot{y}_{sc}.$
(9)

The vertical coordinate of the instantaneous capture point lies on the DRS.
For a horizontal, flat DRS with a constant horizontal velocity but a varying
vertical velocity (e.g., elevators), the vertical coordinate of the capture
point will be zero, i.e., $z_{cap}=0$.

The capture point $(x_{cap},y_{cap},z_{cap})$ is the position of the new
support foot relative to the previous one. Note that the capture point no
longer exists when the term $(\ddot{z}_{wp}+g)$ is negative (i.e., the robot
loses contact with the DRS).

Analogous to the capture point on a stationary surface, if the LIPM
instantaneously places the far end of its leg (point $S$) at
$(x_{cap},y_{cap},z_{cap})$, it will come to a complete stop relative to the
new support point within the immediate successive step.

The expression of the horizontal coordinates of the capture point in (9)
explicitly considers the surface motion that was not incorporated in the
previous derivations of the capture point on a stationary surface [3, 21].
Note that when $\ddot{z}_{ws}=0$ (i.e., the surface is stationary or moving
with a constant speed in both the horizontal and vertical directions), the
proposed capture point on a DRS in (9) reduces exactly to the capture point on
a stationary surface.

###  III-B Capture-Point based Online Planning

The objective of the proposed footstep planner is to generate a robot’s
desired foot placement for maintaining a certain desired walking motion
instead of forcing the robot to stop. To this end, the capture point
$(x_{cap},y_{cap},z_{cap})$ is further extended to determine the desired
footsteps that help sustain the planned robot movement. For simplicity, the
desired motion is specified by a constant base velocity
$\mathbf{v}_{d}=(v_{dx},v_{dy},0)$ relative to point $S$.

We choose to design the desired horizontal foot placement
($\tilde{x}_{step},\tilde{y}_{step}$) as the difference between the capture
point associated with the current relative CoM velocity
$\mathbf{\dot{r}}_{sc}$ and a fictitious capture point corresponding to the
desired velocity $\mathbf{v}_{d}$ [2] multiplied by a positive gain
$K_{step}$:

$\displaystyle\tilde{x}_{step}$
$\displaystyle:=K_{step}\sqrt{\frac{z_{0}}{(\ddot{z}_{ws}+g)}}(\dot{x}_{sc}-v_{dx})~{}\text{and}~{}$
(10) $\displaystyle\tilde{y}_{step}$
$\displaystyle:=K_{step}\sqrt{\frac{z_{0}}{(\ddot{z}_{ws}+g)}}(\dot{y}_{sc}-v_{dy}).$

The gain $K_{step}$ can be used to tune how much the difference between the
current and desired CoM velocities affects the desired footstep location, thus
allowing us to consider practical hardware limitations (e.g., finite step
duration and limited leg length) in footstep planning.

To ensure that the robot moves with $v_{dx}$ and $v_{dy}$ in the horizontal
directions, Raibert’s heuristic [2] is used along with the instantaneous
capture point in (10) to plan the desired foot placement:

${x}_{step}={x}_{hip}+\frac{T}{2}v_{dx}+\tilde{x}_{step},~{}{y}_{step}={y}_{hip}+\frac{T}{2}v_{dy}+\tilde{y}_{step},$
(11)

where $(x_{hip},y_{hip})$ is the horizontal coordinates of the swing hip joint
and $T$ is the desired gait period.

The $z$-coordinate of the desired foot placement on a horizontal, flat DRS is
${{z}}_{step}=\tilde{z}_{step}=0$.

Given the desired CoM trajectory and foot placement, we utilize inverse
kinematics to generate the desired full-body trajectories. Since whole-body
trajectory generation is not the focus of this study, its details are omitted
for brevity.

##  IV Instantaneous QP-based Controller Design

This section introduces the proposed QP-based controller design that provably
tracks the desired walking motions $\mathbf{h}_{d}(t)$ on a DRS and guarantees
the feasibility for ground contact forces. The QP-based controller is
synthesized based on full-order dynamics of quadrupedal walking on a DRS, and
explicitly addresses the surface movement through the formulation of the
holonomic and the friction cone constraints.

Both QP and MPC can be used to formulate a controller with feasibility
guarantees. We choose to use QP because its computational load is low enough
for the typical real-time controller implementation rate (e.g., $500$ Hz or
higher) even when the QP is formulated based upon a full-order robot model. In
contrast, an MPC synthesized based on a full-order model typically computes
its output at a much lower rate, e.g., $30$Hz.

The key steps to the proposed QP formulation include the selection of
optimization variables and the design of the cost function and constraints,
which are explained next.

Optimization variables. As the control objective is to achieve reliable
trajectory tracking with the ground contact constraints respected, we define
the optimization variables $\mathbf{x}$ as the control input $\mathbf{u}$ and
ground contact force $\mathbf{F}_{f}$, i.e.,
$\mathbf{x}=[\mathbf{u}^{T},\mathbf{F}_{c}^{T}]^{T}$.

Cost function. The quadratic cost function to be minimized is set as
$\mathbf{x}^{T}\mathbf{Q}\mathbf{x}+(\tilde{\mathbf{x}}-\mathbf{x})^{T}\mathbf{W}(\tilde{\mathbf{x}}-\mathbf{x})$
with $\tilde{\mathbf{x}}$ the values of the optimization variables at the
previous optimization step. The weighting matrix $\mathbf{Q}$ is a symmetric,
positive definite matrix with an appropriate dimension. The minimization of
$\mathbf{x}^{T}\mathbf{Q}\mathbf{x}$ is to reduce the demanded joint torques.
The other weighting matrix, $\mathbf{W}$, is a symmetric, positive definite
matrix that can be used to penalize sharp variations in some of the
optimization variables.

Equality time-varying holonomic constraints: To ensure that the support feet
do not move relative to the surface, we consider the time-varying holonomic
contraint in (2) that explicitly contains the surface movement.

From the continuous-phase dynamics in (1), we have
$\mathbf{\ddot{q}}:=\mathbf{M}^{-1}(\mathbf{Bu}+\mathbf{J}_{f}^{T}\mathbf{F}_{f}-\mathbf{C}).$
Substituting it in the holonomic constraint in (2) gives:
$\mathbf{J}_{f}\mathbf{M}^{-1}(\mathbf{Bu}+\mathbf{J}_{f}^{T}\mathbf{F}_{f}-\mathbf{C})+\mathbf{\dot{J}}_{f}\dot{\mathbf{q}}=\ddot{\mathbf{r}}_{p}(t),$
which, upon rearrangement, gives:
$\mathbf{A}_{eq1}\mathbf{x}=\mathbf{b}_{eq1},$ where
$\mathbf{A}_{eq1}:=[\mathbf{J}_{f}\mathbf{M}^{-1}\mathbf{B},~{}\mathbf{J}_{f}\mathbf{M}^{-1}\mathbf{J}_{f}^{T}]$
and
$\mathbf{b}_{eq1}:=\mathbf{J}_{f}\mathbf{M}^{-1}\mathbf{C}-\mathbf{\dot{J}}_{f}\dot{\mathbf{q}}+\ddot{\mathbf{r}}_{p}(t)$.

Equality constraints based on input-output linearizing control: To realize
accurate tracking of the desired base and swing-foot trajectories
$\mathbf{h}_{d}(t)$, the controller takes the form of an input-output (I-O)
linearizing control law [18] with the output function defined as the tracking
error.

Let $\mathbf{h}$ denote the variables of interest (i.e., base pose and swing
foot position). Then, the output function representing the tracking error is
given by: $\mathbf{y}:=\mathbf{h}(\mathbf{q})-\mathbf{h}_{d}(t).$ Then,
$\mathbf{\dot{y}}=\frac{\partial{\mathbf{h}}}{\partial\mathbf{q}}\dot{\mathbf{q}}-\dot{\mathbf{h}_{d}}$
and
$\mathbf{\ddot{y}}=\frac{\partial}{\partial\mathbf{q}}(\frac{\partial{\textbf{h}}}{\partial\mathbf{q}}\dot{\mathbf{q}})\dot{\mathbf{q}}+\frac{\partial{\mathbf{h}}}{\partial\mathbf{q}}\ddot{\mathbf{q}}-\ddot{\mathbf{h}}_{d}$.
Define $\mathbf{J}_{h}:=\frac{\partial{\mathbf{h}}}{\partial\mathbf{q}}$, and
then
$\mathbf{\dot{J}}_{h}=\frac{\partial}{\partial\mathbf{q}}(\frac{\partial{\textbf{h}}}{\partial\mathbf{q}}\dot{\mathbf{q}})$.
We get the output function dynamics as
$\mathbf{\ddot{y}}=\mathbf{J}_{h}\ddot{\mathbf{q}}+\mathbf{\dot{J}}_{h}\dot{\mathbf{q}}-\ddot{\mathbf{h}}_{d}$.

To stabilize the output function dynamics, we utilize a proportional
derivative (PD) feedback term
$\mathbf{\ddot{y}}=-\mathbf{K}_{p}\mathbf{y}-\mathbf{K}_{d}\mathbf{\dot{y}}=:\mathbf{v},$
which, combined with the expression of $\ddot{\mathbf{q}}$ in (IV), yields:
$\mathbf{J}_{h}\mathbf{M}^{-1}(\mathbf{Bu}+\mathbf{J}_{c}^{T}\mathbf{F}_{c}-\mathbf{C})+\mathbf{\dot{J}}_{h}\dot{\mathbf{q}}+\ddot{\mathbf{h}}_{d}(t)=\mathbf{v}.$
It can be compactly written as $\mathbf{A}_{eq2}\mathbf{x}=\mathbf{b}_{eq2},$
where
$\mathbf{A}_{eq2}:=[\mathbf{J}_{hc}\mathbf{M}^{-1}\mathbf{B},~{}\mathbf{J}_{h}\mathbf{M}^{-1}\mathbf{J}_{f}^{T}]$
and
$\mathbf{b}_{eq2}:=\mathbf{J}_{hc}\mathbf{M}^{-1}\mathbf{C}-\mathbf{\dot{J}}_{h}\dot{\mathbf{q}}-\ddot{\mathbf{h}}_{d}(t)+\mathbf{v}$.

Inequality ground-contact constraints and joint torque limits. To secure
contacts between support feet and the surface, the QP incorporates the
friction cone (i.e., no foot slipping) and unilateral (i.e., no feet
penetrating the surface) constraints $\mathbf{F}_{f}\in\mathcal{F}_{gc},$
where $\mathcal{F}_{gc}$ is the set of ground-contact forces satisfying the
friction cone and unilateral constraints. Also, the solution to the QP must
respect the torque limits; i.e., the control input $\mathbf{u}$ must be within
the admissible torque set, i.e., $\mathbf{u}\in U$.

QP-based control law. With the cost function and constraints designed, the
proposed QP can be compactly expressed as:

$\displaystyle\min_{\mathbf{x}}$
$\displaystyle\mathbf{x}^{T}\mathbf{Q}\mathbf{x}+(\mathbf{x}_{p}-\mathbf{x})^{T}\mathbf{W}(\mathbf{x}_{p}-\mathbf{x})$
(12) subject to
$\displaystyle\mathbf{A}_{eq1}\mathbf{x}=\mathbf{b}_{eq1},~{}\mathbf{A}_{eq2}\mathbf{x}=\mathbf{b}_{eq2},~{}\mathbf{F}_{f}\in\mathcal{F}_{gc},~{}u\in\mathbf{U}.$

Stability property. As proven in our previous study [18], an input-output
linearizing controller with the proposed form can guarantee the (local)
asymptotic stability of the desired trajectory $\mathbf{h}_{d}(t)$ for the
hybrid, time-varying walking process, if the PD gains are chosen to render a
sufficiently fast error convergence rate during continuous phases.

##  V Simulations

This section presents MATLAB simulation results that validate the
effectiveness of the proposed control approach.

Robot. The simulated robot is a Laikago quadruped developed by Unitree
Robotics. Its total mass is 25 kg, and each leg weighs 2.9 kg. Each leg has
three independently actuated joints (i.e., hip-roll, hip-pitch, and knee-
pitch) with torque limits of 20 Nm, 55 Nm, and 55 Nm, respectively. The
length, width, and height of the robot are 0.56m, 0.35m, and 0.60m,
respectively.

Figure 2: System setup that the simulator emulates. The treadmill (①) has a
split belt (④) that moves at a constant speed of $8$cm/s while the treadmill
undergoes a sinusoidal rocking motion about the horizontal axis (③). The
quadruped is a Laikago robot (②) developed by Unitree Robotics.

Surface. The simulated DRS is an actuated platform (see Fig. 2) simultaneously
experiencing: a) a whole-body sinusoidal pitching motion, with an amplitude of
$\pm 5^{\circ}$ and frequency of $0.5$ Hz and b) a constant surface
translating motion. This DRS reasonably satisfies the assumption underlying
the proposed capture point extension, because the horizontal velocity of the
surface is approximately constant due to the small pitching amplitude. The
surface acceleration in the vertical direction is still relatively significant
for validating the proposed method.

QP-based controller. The PD gains are $\mathbf{K}_{p}=120\cdot\mathbf{I}_{9}$
and $\mathbf{K}_{d}=22\cdot\mathbf{I}_{9}$, where $\mathbf{I}_{n}$ is an
$n\times n$ identity matrix. The weighting matrices in the cost function are
set as:
$\mathbf{Q}=\text{blockdiag}(1000\cdot\mathbf{I}_{12},~{}\mathbf{I}_{9},)$ and
$\mathbf{W}=\text{blockdiag}(10^{-1}\cdot\mathbf{I}_{12},~{}10^{-4}\cdot\mathbf{I}_{9},)$.
The torque bounds are set based on the robot’s hardware limits. In the
unilateral constraint, the lower bound of the normal ground reaction force at
a support foot is set as $1$ N. In the friction cone constraint, the friction
coefficient is chosen as 0.5.

Planner. The inputs to the planner are: a) gait cycle $T=2$ s, b) maximum
swing foot height $15$ cm, c) walking speed $8$ cm/s, and d) the surface
movement as explained earlier. The gain used to tune the desired foot
placement is set as $K_{step}=1$.

Figure 3: Simulation results obtained from unperturbed walking (Case 1): (a)
joint torques, (b) ground reaction forces, and (c) base trajectories.  Figure
4: Simulation results obtained in the presence of lateral perturbations (Case
2): (a) joint torques, (b) ground reaction forces, and (c) base trajectories.
Figure 5:  Comparative simulation results on the satisfaction of the friction
cone constraint at the front right and rear right feet: (a) our previous
offline planning and I-O linearizing control [18] and (b) the proposed online
planning and QP-based control.

Case 1 (Unperturbed walking). Figure 3 shows the simulation results of DRS
walking without perturbations under the proposed method. As shown in Fig. 3
(c), the base orientation mildly varies, indicating that the proposed method
can sustain stable walking on the DRS in the absence of perturbations. Also,
Figs. 3 (a) and (b) show that the joint torques are within the limits and that
the ground contact constraints are satisfied.

Case 2 (Perturbed walking). Figure 4 displays the simulation results of robot
walking on the DRS under sinusoidal perturbation forces applied at the center
of the robot’s base. The magnitude of the external force is:

$F_{p}(t)=\begin{cases}30\sin(2\pi t)~{}\text{N}&\text{if $2\text{s}\leq t\leq
4\text{s}$}\\\ 0~{}\text{N}&\text{otherwise}\end{cases}$ (13)

The robot’s base trajectories are notably perturbed between 2 sec and 4 sec,
as shown in subplot (c). In particular, the lateral base position trajectory
is significantly displaced with a peak deviation of 0.12 m. Yet, the proposed
control approach is able to sustain stable walking during the perturbations,
as well as to drive the actual based trajectories to converge back to the
desired motion after the perturbation is over. Also, the joint torque and
ground contact constraints are respected throughout the walking process, as
demonstrated in subplots (a) and (b).

Case 3 (Comparative results). Figure 5 shows the simulation results of our
previous offline planning and I-O linearizing control [18] under unperturbed
walking. Although this controller explicitly addresses the surface motion, it
relies on offline planning and does not explicitly guarantee the feasibility
of ground contact forces. The unilateral ground contact constraint is indeed
respected by the controller during simulated walking, as shown in Fig. 5 (a).
Yet, the friction cone constraint is violated at the front right and rear
right feet. In contrast, the proposed online planning and QP-based control
approach is able to respect all ground contact constraints during the
unperturbed walking (subplot (a) of Fig. 5), as well as during perturbations
(subplots (b) in Fig. 4).

##  VI Conclusion

This paper has introduced an online planner and QP-based controller design
that sustains stable quadrupedal walking on a DRS (e.g., elevators and ships)
by explicitly considering the surface motion and guaranteeing the feasibility
of ground contact constraints. The key element of the online planner is a
real-time footstep generator synthesized based on the extension of the concept
of capture point from stationary surfaces to a DRS. The capture point was
derived based on a 3-D LIPM model walking on a DRS that undergoes a constant
horizontal velocity and a varying vertical motion. The controller was
formulated as a QP based on the full-order robot model that captures the
support foot movement caused by the surface motion. The effectiveness of the
proposed approach in walking stabilization and disturbance rejection was
demonstrated via simulations of 3-D quadrupedal walking over a DRS. The future
work is dedicated to experimental validation of the proposed approach on a
physical quadrupedal robot walking on a real-world nonstationary surface with
an uncertain motion.

## References

  * [1] M. Hutter, C. Gehring, D. Jud, A. Lauber, C. D. Bellicoso, V. Tsounis, J. Hwangbo, K. Bodie, P. Fankhauser, M. Bloesch _et al._ , “ANYmal-a highly mobile and dynamic quadrupedal robot,” in _Proc. of IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2016, pp. 38–44.
  * [2] G. Bledt, M. J. Powell, B. Katz, J. Di Carlo, P. M. Wensing, and S. Kim, “MIT Cheetah 3: design and control of a robust, dynamic quadruped robot,” in _Proc. of IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2018, pp. 2245–2252.
  * [3] J. Pratt, J. Carff, S. Drakunov, and A. Goswami, “Capture point: A step toward humanoid push recovery,” in _Proc. of IEEE-RAS International Conference on Humanoid Robots_ , 2006, pp. 200–207.
  * [4] T. Koolen, T. De Boer, J. Rebula, A. Goswami, and J. Pratt, “Capturability-based analysis and control of legged locomotion, Part 1: theory and application to three simple gait models,” _The International Journal of Robotics Research_ , vol. 31, no. 9, pp. 1094–1113, 2012.
  * [5] S. Kajita, F. Kanehiro, K. Kaneko, K. Fujiwara, K. Harada, K. Yokoi, and H. Hirukawa, “Biped walking pattern generation by using preview control of zero-moment point,” in _Proc. of IEEE International Conference on Robotics and Automation_ , vol. 2, 2003, pp. 1620–1626.
  * [6] Y. Gu, B. Yao, and C. S. G. Lee, “Exponential stabilization of fully actuated planar bipedal robotic walking with global position tracking capabilities,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 140, no. 5, 2018\.
  * [7] Y. Gu, B. Yao, and C. G. Lee, “Straight-line contouring control of fully actuated 3-D bipedal robotic walking,” in _Proc. Amer. Contr. Conf._ , 2018, pp. 2108–2113.
  * [8] Y. Gao and Y. Gu, “Global-position tracking control of multi-domain planar bipedal robotic walking,” in _Proc. of ASME Dynamic Systems and Control Conference_ , 2019.
  * [9] ——, “Global-position tracking control of multi-domain planar bipedal robotic walking,” in _Proc. ASME Dyn. Syst. Contr. Conf._ , 2019.
  * [10] V. Barasuol, J. Buchli, C. Semini, M. Frigerio, E. R. De Pieri, and D. G. Caldwell, “A reactive controller framework for quadrupedal locomotion on challenging terrain,” in _Proc. of IEEE International Conference on Robotics and Automation_ , 2013, pp. 2554–2561.
  * [11] C. Li, T. Zhang, and D. I. Goldman, “A terradynamics of legged locomotion on granular media,” _Science_ , vol. 339, no. 6126, pp. 1408–1412, 2013.
  * [12] X. Xiong, A. D. Ames, and D. I. Goldman, “A stability region criterion for flat-footed bipedal walking on deformable granular terrain,” in _Proc. of IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2017, pp. 4552–4559.
  * [13] S. Galeani, L. Menini, and A. Potini, “Robust trajectory tracking for a class of hybrid systems: An internal model principle approach,” _IEEE Transaction on Automatic Control_ , vol. 57, no. 2, pp. 344–359, 2011.
  * [14] Y. Zheng and K. Yamane, “Ball walker: A case study of humanoid robot locomotion in non-stationary environments,” in _Proc. of IEEE International Conference on Robotics and Automation_ , 2011, pp. 2021–2028.
  * [15] C. Yang, B. Zhang, J. Zeng, A. Agrawal, and K. Sreenath, “Dynamic legged manipulation of a ball through multi-contact optimization,” in _Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , Oct 2020, pp. 7513–7520.
  * [16] F. Asano, “Modeling and control of stable limit cycle walking on floating island,” in _Proc. of IEEE International Conference on Mechatronics (ICM)_ , 2021, pp. 1–6.
  * [17] A. Iqbal, Z. Mao, and Y. Gu, “Modeling, analysis, and control of slip running on dynamic platforms,” _ASME Letters in Dynamic Systems and Control_ , vol. 1, no. 2, 2021.
  * [18] A. Iqbal, Y. Gao, and Y. Gu, “Provably stabilizing controllers for quadrupedal robot locomotion on dynamic rigid platforms,” _IEEE/ASME Transaction on Mechatronics_ , vol. 25, no. 4, pp. 2035–2044, 2020.
  * [19] J. Englsberger, C. Ott, M. A. Roa, A. Albu-Schäffer, and G. Hirzinger, “Bipedal walking control based on capture point dynamics,” in _Proc. of IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2011, pp. 4420–4427.
  * [20] M. Morisawa, S. Kajita, F. Kanehiro, K. Kaneko, K. Miura, and K. Yokoi, “Balance control based on capture point error compensation for biped walking on uneven terrain,” in _Proc. of IEEE-RAS International Conference on Humanoid Robots (Humanoids 2012)_ , 2012, pp. 734–740.
  * [21] S. Caron, A. Escande, L. Lanari, and B. Mallein, “Capturability-based pattern generation for walking with variable height,” _IEEE Transactions on Robotics_ , vol. 36, no. 2, pp. 517–536, 2019.
  * [22] H.-M. Joe and J.-H. Oh, “Balance recovery through model predictive control based on capture point dynamics for biped walking robot,” _Robotics and Autonomous Systems_ , vol. 105, pp. 1–10, 2018.
  * [23] M. Focchi, A. Del Prete, I. Havoutis, R. Featherstone, D. G. Caldwell, and C. Semini, “High-slope terrain locomotion for torque-controlled quadruped robots,” _Autonomous Robots_ , vol. 41, no. 1, pp. 259–272, 2017.
  * [24] K. A. Hamed, J. Kim, and A. Pandala, “Quadrupedal locomotion via event-based predictive control and QP-based virtual constraints,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 3, pp. 4463–4470, 2020.