# Search for gamma-ray emission from a galactic supernova with the
anticoincidence system of SPI

M. Caixach,1,2 P. Jean,3,4 J. Isern1,2,5 and E. Bravo6
1Institut de Ciències de l'Espai (ICE, CSIC), Campus UAB, C/ de Can Magrans
s/n, 08193 Cerdanyola del Vallès, Spain
2Institut d'Estudis Espacials de Catalunya (IEEC), C/ Gran Capità 2-4, 08034
Barcelona, Spain
3Universiteé de Toulouse, UPS-OMP, IRAP, 31028 Toulouse, France
4Institut de Recherche en Astrophysique et Planétologie, 9 Av colonel Roche,
BP44346, 31028 Toulouse Cedex 4, France
5Fabra Observatory, Royal Academy of Sciences and Arts of Barcelona (RACAB),
La Rambla 115, 08002 Barcelona, Spain
6E.T.S.A.V., Univ. Politècnica de Catalunya, c/Pere Serra 1-15, 08173 Sant
Cugat del Vallès, Spain E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The detection of the very early gamma-emission of a Type Ia supernova (SNIa)
could provide a deep insight on the explosion mechanism and nature of the
progenitor. However this has not been yet possible as a consequence of the
expected low luminosity and the distance at which all the events have occurred
up to now. A SNIa occurring in our Galaxy could provide a unique opportunity
to perform such measurement. The problem is that the optical flux would
probably be so attenuated by interstellar extinction that would prevent
triggering the observations with gamma-spectrometers at the due time. In this
paper we analyse the possibility of using the anticoincidence system (ACS) of
the spectrometer SPI on board of the INTEGRAL space observatory for detecting
the early gamma-ray emission of a SNIa as a function of the explosion model
and distance as well as of pointing direction. Our results suggest that such
detection is possible at about 6 - 12 days after the explosion and, at the
same time, we can discard missing any hidden explosion during the lifetime of
INTEGRAL.

###### keywords:

supernovae: general – transients: supernovae – gamma-rays: stars

††pubyear: 2021††pagerange: Search for gamma-ray emission from a galactic
supernova with the anticoincidence system of SPI–Search for gamma-ray emission
from a galactic supernova with the anticoincidence system of SPI

##  1 Introduction

Type Ia supernovae (SNIa) are one of the most energetic explosion events in
the universe, with luminosity magnitudes comparable to the ones from galaxies.
These types of supernova are identified by the lack of H lines and for having
strong Si lines in the early spectra (Filippenko, 1997). The optical light
curve of SNIa peaks at about 20 days after the explosion and decreases 3
magnitudes after one month. They are the consequence of the thermonuclear
explosion of a C/O white dwarf (WD) that accretes mass from a companion in a
close binary (Hoyle & Fowler, 1960). The discovery of several subtypes of SNIa
these last years may be explained by SNIa possibly having more than one kind
of progenitor scenarios and explosion mechanisms (Hillebrandt & Niemeyer,
2000; Maoz et al., 2014; Isern et al., 2021).

Even though the progenitor system and explosion mechanism of such events are
not yet proved, one thing is known to be sure about SNIa: their optical
luminosities are powered by the decay chain $~{}^{56}$Ni $\rightarrow$
$~{}^{56}$Co $\rightarrow$ 56Fe. $~{}^{56}$Ni is synthesised during the rapid
C/O burning and it fuels the explosion by the release of high energy gamma-ray
photons that heat the debris and powers the optical and infrared light curves
of the event (see Colgate & McKee, 1969). The detection of gamma-ray lines
from this decay can give clues about the kinematics and morphology of the
explosion (Burrows & The, 1990; Gomez-Gomar et al., 1998; Jean et al., 1999;
Milne et al., 2004; Isern et al., 2008). The most prominent lines from 56Ni
decay are 158, 480, 750, and 812 keV and from 56Co decay are 847 and 1238 keV.
They can be detected once the ejecta is thin enough in order for the gamma-ray
photons to escape without getting thermalized. The decay of 56Co has a mean
lifetime of $\sim$111 days and its lines can be detected right after the
maximum of the optical light curve and the following months. The intensity and
broadening of the lines give information about the inner layers composition.
On the other hand, the decay of 56Ni has a mean lifetime of $\sim$8.8 days
which makes the detection of its gamma-ray lines somehow trickier, as for some
explosion models the opacity of the debris is too high for the photons to be
able to escape without being absorbed. However, if there are radioactive
elements in the outer layers, this gamma-ray emission can give the needed
light to the study of the early stage of SNIa.

The closer SNIa ever detected in gamma-rays is SN2014J. It was discovered by
Fossey et al. (2014). It was found in M82 at a distance of 3.5 Mpc. SPI
instrument on board of INTEGRAL (INTErnational Gamma-RAy Laboratory) collected
gamma-rays from day 16.3 to 164 after the explosion. The detection of 56Ni and
56Co lines of this nearby event allowed to prove the usefulness of gamma
detection to diagnose the dynamics and composition of the ejecta (Churazov et
al., 2014; Diehl et al., 2014; Isern et al., 2014; Churazov et al., 2015;
Isern et al., 2016). The analysis of these observations not only confirmed the
hypothesis that the light curve was powered by the disintegration of the 56Ni
radioactive chain but also allowed to compute in a direct way the total amount
of 56Ni synthesized during the event. Furthermore they showed, unexpectedly,
that $\sim 0.05$ M⊙ of 56Ni was present in the outer layers forming a non-
spherical structure. The non-spherical structure was indicated by the display
of redshifted and high intensity lines of 56Ni in the early spectrum.

A Galactic supernova would be ideal for more detailed studies, specially if
detected during the rising epoch of the light curve. The expected rate is
estimated to be $1.4^{+1.4}_{-0.8}$ events per century (Adams et al., 2013) in
the case of SNIa, and it is important to realize that the Galactic plane is
transparent to gamma-rays for which reason events severely obscured by dust
can be detected in gamma. Therefore, taking into account that INTEGRAL has
been operating during $\sim 20$ years the probability that such an event has
occurred during this period is not negligible.

The anticoincidence system of SPI (ACS) is used not only to reduce the
background of SPI but also to monitor the sky for possible sources thanks to
its large field of view and its spectroscopic capabilities (Vedrenne et al.,
2003). Jean et al. (1999) proposed the use of the SPI-ACS to search for
Galactic classical novae, Von Kienlin et al. (2001) to locate Gamma Ray Bursts
(GRBs), and Rodríguez-Gasén et al. (2014) to study solar flares. The scope of
this paper is to study the ability of SPI-ACS for detecting Galactic Type Ia
supernovae as soon as possible (Wang et al., 2019) and, in addition, to
scrutinize the presence of its signature in the already recorded ACS data.

This paper is organized as follows. We present the supernova models used to
simulate the detection by the ACS in section 2. Section 3 describes the ACS
data used for this study. The simulation method and the results of the
analyses to estimate the sensitivity of the ACS are presented in section 4.
Section 5 is dedicated to the search for SNIa signature in the existing
SPI/ACS data. Finally, in section 6 we present the discussion and conclusions.

##  2 Supernova models

Two models, DDTe and W7, of SNIa have been used for testing the possibility of
detection by ACS/SPI. The spectra as a function of time of these models were
obtained with 1D simulations as in Gomez-Gomar et al. (1998). These models
were compared with those obtained from SN 2011fe and SN 2014J by the
instruments on board of INTEGRAL (Isern et al., 2013, 2014). Both are
compatible with the upper limits deduced for SN 2011fe but the SN 2014J case
demands further considerations since the early spectra presented some features
that were better interpreted introducing non-spherical structures.

The DDTe is a sub-luminous SNIa model. It describes a situation in which the
burning front starts as a deflagration that turns out into a detonation when
reaches a critical density. As a result, material is processed into Fe-group
and into intermediate-mass elements. This model ejects 0.51 M⊙ of 56Ni with a
kinetic energy of 1.09 $\times$ 1051 ergs. The W7 model is a normal luminous
SNIa model (Nomoto et al., 1984). In this model a flame ignited near the
center propagates subsonically, thus allowing the expansion of the star and
leaving some unburnt material. The W7 model, in particular, is characterized
by a small portion of unburnt C-O and a large mass of burnt intermediate-mass
elements. This model ejects 0.56 M⊙ of 56Ni with a kinetic energy of 1.24
$\times$ 1051 ergs. Figure 1 displays the light curve of both models for SNIa
at 8 kpc. They have been obtained by integrating the spectrum models from 6
keV to 3684 keV for a SNIa at 8 kpc.

Figure 2 displays the gamma ray spectra of both models in a logarithmic scale
during the 9th day after the explosion. Even at such early stages the 750 keV
and 812 keV 56Ni lines decay are visible including several features from the
56Co decay.

The slow rising signal of the DDTe model as compared to that of W7 is a
convenient property for exploring if the differences in luminosity have an
impact on the detection sensitivity. In order to analyse the detection
sensitivity as a function of SNIa distance, models were placed at distances in
a range of 4 kpc to 16 kpc. For each model the output flux is provided with an
energy resolution of 0.5 keV, in the interval 6 keV to 3684 keV. The spectra
were calculated with a frequency of one day from days 5 to 200 after the
explosion, although just the spectra derived for the first days were used and
interpolated, as the main goal is to study the detection of the rising signal
as soon as possible. For further details about the models data preparation see
section 4.

Figure 1: Gamma ray flux variation for model W7 (blue) and DDTe (red) at a
distance of 8 kpc from day 5th to 200th after the explosion. Figure 2: Gamma-
ray spectra for day 9th of the models, W7 (blue) and DDTe (red) at 8 kpc. A
label on top of each line is added to indicate if they are from 56Ni or 56Co
decay.

##  3 The anticoincidence system (ACS) of SPI

The ACS of SPI aims first to reduce the instrumental background in the
germanium detectors due to charged particles and gamma-rays coming from
outside the spectrometer’s field of view. It is composed of an active shield
made with 91 scintillator blocks111Only 89 blocks are active (see Savchenko et
al., 2012) in bismuth germanate (BGO) and a plastic scintillator, optically
coupled to photomultiplier tubes (see Von Kienlin et al., 2001, for a detailed
description of the SPI ACS). When a particle releases energy in a scintillator
block, a veto signal is generated and transmitted to the on-board digital
acquisition system. By this way, the active shield allows to remove events
produced in coincidence by particles that deposit energy in the shield and in
the germanium detectors. The ACS veto rate is recorded with a sampling period
of 50 ms and is sensitive to photons releasing energy $\gtrsim$ 80 keV in a
scintillator block.

The ACS of SPI is also used to monitor astrophysical sources. For instance,
its large detection area allows the detection of gamma-ray bursts and it is
part of the INTEGRAL burst alert system providing the gamma-ray burst location
in the sky by triangulations with other space-borne gamma-ray burst monitors
(Von Kienlin et al., 2001; Von Kienlin et al., 2003; Savchenko et al., 2017).
A giant outburst from the soft gamma-ray repeater SGR 1806-20 was discovered
in 2004 by the analysis of the ACS rate (Borkowski et al., 2004; Mereghetti et
al., 2005). Rodríguez-Gasén et al. (2014) and Gros et al. (2004) explored the
capabilities of the SPI ACS to study solar flares. Jean et al. (1999) proposed
to use SPI ACS to search for the hard X-ray and gamma-ray emission from
classical novae in our Galaxy (see also Siegert et al., 2018). For this
purpose, they computed the effective area of the SPI ACS with GEANT 3 for
energies ranging from 80 keV to 500 keV and a rough angular binning.

For the present study, we re-calculated the response of the SPI ACS with GEANT
4, using the same INTEGRAL mass model, for energies ranging from 80 keV to 3.5
MeV and with a binning of 10 degrees for the zenithal and azimuthal angles
(see Fig. 3). The mass model of INTEGRAL (see Fig 4) was developped by
Ferguson et al. (2003). It includes the SPI mass model developed by the SPI
collaboration (see Sturner et al., 2003). For a given photon energy, the
response was calculated assuming a point source at infinity (i.e. parallel
beam) at a given angular position, by counting the number of events with
energy deposit larger than the threshold of 80 keV. The effective area is
equal to the number of triggering events divided by the number of simulated
photons and multiplied by the exposed area, perpendicular to the parallel
beam, in the simulation. The values of the effective area obtained in this way
are close (within $\sim\pm$30%) to the previous ones; e.g. $\sim$397 cm2 for E
$>$ 100 keV at a zenith angle of 100 deg and an azimuth angle of 0 deg,
compared to the value of 340 cm2 close to the SGR 1806-20 position (see
Mereghetti et al., 2005). They are also in agreement with the effective areas
calculated by Rodríguez-Gasén et al. (2014) (e.g. differences are less than
5%) who used the same INTEGRAL mass model to estimate the gamma-ray flux from
some selected solar flares.

Figure 3: Effective area (cm2) of SPI for energies of 0.081, 0.8 and 3.5 MeV
for all pointing directions in the sky. The effective area is lower for
regions where $z=0$ deg and the energies are low, this is due to SPI being
masked by IBIS.  Figure 4: Mass model of INTEGRAL used to calculate the
response of the anticoincidence system of SPI (in green) with GEANT 4 (see
text).

The anticoincidence system of SPI provides also the saturating event rate,
defined as the rate of events released in the BGO blocks with energy $\gtrsim$
150 MeV. The ACS saturating event rate, unlike the ACS rate, is not affected
by low energy particles from the radiation belt, from solar flares or from
gamma-ray sources that emit mainly in the low energy gamma-ray range (i.e. E
$<$ 100 MeV). Consequently, that makes it a good tracer for monitoring the
cosmic-ray intensity at the spacecraft level. When the spacecraft is not
affected by these low energy events, the temporal behaviour of the ACS
saturating rate is concordant with the total ACS rate. Figure 5 shows the
average ACS rates and the average ACS saturating rates for some revolutions
since the launch of INTEGRAL. In order to avoid the effects of radiation
belts, the rates shown are the rates averaged over the first to the third
quarter of each revolution. Figure 5 shows that the total and saturating rates
follow a linear relation, although the linear coefficients may slightly change
due to slight changes of the high energy particle spectrum impinging the ACS.
The saturating rate can be therefore used to trace the variation of the ACS
rate without source contribution and be used as background model to our study
(see section 4.2).

Figure 5: Total and saturating event rates of the anticoincidence system of
SPI averaged for some revolutions since the launch of INTEGRAL (see text).

##  4 Simulations and analyses

###  4.1 Simulated count rate

The count rate (CR), counts per second, from a source located at a zenithal
and azimuthal angles $\theta$ and $\varphi$ in the SPI frame at the instant
$t$ can be defined as:

$R_{\rm sn}(t)=\sum_{E_{i}}F(E_{i},t)\cdot
S_{eff}[E_{i},\theta(t),\varphi(t)]\cdot\Delta E$ (1)

where $F(E_{i},t)$ (counts/cm2/s/keV) is the flux for a given energy and time,
$S_{eff}[E_{i},\theta(t),\varphi(t)]$ is the effective area (cm2 – i.e. the
response matrix) of the ACS and $\Delta$E is the size of the energy bin (keV)
of the spectrum. Overall it provides the rate in counts/s, $R_{\rm sn}$, that
would be produced by a SN model. Two SNIa models (see section 2) have been
simulated, the W7 model and a delayed detonation model (DDTe), with an energy
binning of 0.5 keV and a time binning of 80 s to provide enough resolution.
Furthermore, the models have been scaled to Galactic distances: 4, 6, 8, 10,
12, 14 and 16 kpc.

The response matrix has been computed by Monte Carlo simulations (see Section
3). It gives the effective area for energies in the range of 81 keV to 3.5
MeV, with an azimuthal resolution of 20∘ and a zenithal resolution of 10∘.

An interpolation of the response matrix is used to take into account the
influence of the pointing variation during each revolution on the effective
area for each energy event reaching the satellite. INTEGRAL provides for each
revolution the right ascension (RA) and declination (Dec) of the x axis of the
instrument and the RA and Dec of the z axis. For running the tests we fix the
RA and Dec of the simulated source and we take into account the angular
evolution of the source with respect to the x and z axis of the instrument to
obtain the local coordinates and finally the effective area.

INTEGRAL provides the ACS data along the lifetime of the satellite for each
orbital revolution. The ACS rate from a chosen revolution is added as a
background to the SN model to make the simulated count rate realistic. We
choose revolutions without strong background variations due to solar flares or
early entrance in radiation belts. Each revolution has a duration between 2 to
3 days and the time of each revolution (and the SN model) is binned in
intervals of 80 s. Before adding the ACS rate to the SN model rate, the first
one needs to be converted from measured rate to true rate by taking into
account the dead time $\delta t$ of 0.6 $\mu$s (Savchenko et al., 2012):

$R_{\rm bck}=\frac{R_{\rm m}}{R_{\rm m}\cdot\delta t-1}$ (2)

where $R_{\rm bck}$ is the true background rate and $R_{\rm m}$ the measured
rate extracted from the data. The rate $R_{\rm sn}$ produced by the SNIa,
computed with equation (1), is added to the true background rate to yield to
the total rate $R_{\rm tot}$. The final simulated rate $R_{\rm sim}$ is
obtained with:

$R_{\rm sim}=\frac{R_{\rm tot}}{1+R_{\rm tot}\cdot\delta t}$ (3)

to take into account the dead time of the ACS.

###  4.2 Detection of the supernova signal

Once the ACS rate produced by the supernova is simulated, we aim to study the
sensitivity of the SPI ACS to detect the earliest gamma ray signature of our
SN models within a galactic distance. Two methods have been tested.

The first one consist of an "ON/OFF". It analyses a previous adjacent window
(or interval) in the data, called the "off" window which represents a null
hypothesis, and the subsequent window with the same size, the "on" window that
represents an hypothetical source (Li & Ma, 1983). A rising of significance is
expected when the "on" window have a rapid increase of flux compared to the
previous "off" window. This method, however, has been discarded as the
significance of the detections are not high enough. Indeed, when that method
is applied, the rate in the "off" period contains the low signal from events
produced in the anticoincidence system by the rising gamma-ray emission from
the SNIa. The background rate is therefore overestimated. The contribution of
the rising signal in the background rate reduces the significance of the
detection.

The second method estimates the detection of the rising signal of the
supernova by fitting the intensity of the supernova rate and the parameters of
a background rate model to the SPI/ACS data. This method takes into account
two assumptions. The first one considers that the background rate model (rbck)
is a linear function of the ACS saturating (Rsat) event rate :

$r_{\rm bck}=\frac{R_{\rm sat}-c}{f}$ (4)

with $c$ and $f$ as free parameters. The saturating event rate (rate ACS
saturating) is the rate of events with energy released in the BGO blocks with
energies $\gtrsim$ 150 MeV (see section 3). The second assumption is that the
rate produced by the rising gamma-ray flux of the supernova can be described
as the following power law:

$r_{\rm sn}=I\cdot(t/5)^{\beta}$ (5)

Where $I$ is the intensity and $t$ is the time of the model. Both are free
parameters. The power law model is justified by the shape of the gamma-ray
light curves when the emission is rising (see Fig.1). The slope of the power
law $\beta$ is equal to 18.9 for the W7 model and 14 for the DDTe one. These
two different values of the slope does not impact very much on the results.

With these two assumptions we obtain a total model $r_{\rm mod}$ = $r_{\rm
bck}$ \+ $r_{\rm sn}$, where $r_{\rm bck}$ is the background model (see Eq. 4)
and $r_{\rm sn}$ is the supernova signal model (see Eq. 5). The best fitting
parameters ($I$, $c$ and $f$) are the ones that minimize the chi square
calculated as:

$\chi^{2}=\sum\left(\frac{R_{\rm sim}-r_{\rm mod}(I,c,f)}{\sigma_{\rm
acs}}\right)^{2}+10\left(f-f_{0}\right)^{2}+\left(\frac{c}{2000}\right)^{2}$
(6)

where $R_{\rm sim}$ is the simulated ACS rate (see Eq. 3) and $\sigma_{\rm
acs}$ the standard deviation of the model and the ACS data. The second and
third terms in the sum are penalty functions that avoid wrong best fitting
values of the parameters $f$ and $c$. Without those penalty functions, the
minimization process could converge toward local chi square minimum yielding
to abnormal values of the parameters of the background model, $f$ and $c$. It
sometimes occurs when the signal of the supernova is strong. The parameter
$f_{0}$ in Eq. 6 is an estimated mean value of the factor that scales the ACS
total rate with the ACS saturated rate (see Fig. 5 and section 3). The factors
10 and 2000 in Eq. 6 were chosen to limit the variations of the f and c on a
reasonable range of values.

From these results we can estimate the detection date. To do so, we estimate
first the quality of the background model ($r_{\rm bck}$) computing a residual
($\delta$):

$\delta=\frac{R_{\rm sim}-r_{\rm bck}}{\sigma_{\rm acs}}$ (7)

Then, this residual is smoothed to avoid a false detection due to a single
short time scale spike in the data.

We consider that the supernova signal is detected when the following
conditions are satisfied: 1) tacs $>$ tstart \+ 0.5 day; 2) $\delta>$ 3\.
Where tacs is the time of the data, the sigma level is set to 3, and tstart is
the date at which the ACS data starts to be recorded at the beginning of an
orbital revolution of INTEGRAL. The addition of half a day in the first
condition aims to remove significant count rate excess induced by high energy
particles when INTEGRAL goes out of the radiation belts, after its perigee.

This procedure is done to test the DDTe SN model and the W7 SN model over 275
revolutions. These revolutions have been taken along the lifetime of INTEGRAL
to find the earliest time of detection. Each of them have the following
available data: ACS rate, ACS saturating rate, rate time, right ascension and
declination of the INTEGRAL X-axis, right ascension and declination of the
INTEGRAL Z-axis. The ACS rate is used as simulated background to the models,
as mentioned in section 4.1. The ACS saturating event rate is used on this
analysis to estimate the background ACS rate and detect the supernova signal,
as mentioned previously in this section. Each revolution has a length of 2 to
3 days that we fixed to $\approx$2 days as we cut the initial hours of data
because most of them still contain some noise from the radiation belt
phase222SPI is switched off between revolutions. A same length of the
supernova model is cut in order to compute the simulated count rate. The
initial time of the computed count rate equals to a time ranging from day 4th
to 6th of the SNIa model, for earlier days the flux’s model is not high enough
to make a contribution to the ACS rate. This initial time, from 4th to 6th, is
chosen randomly to give more realistic circumstances to the analysis.

When the count rate with a length of 2 days is analyzed two things may happen:
a supernova signal detection is found or not. If a detection is found in the
count rate then the revolution number, the time detection and its significance
are saved. In case it is not, we proceed to analyze the consecutive 2 days of
the computed count rate, which would account for the 2 following days of the
SNIa model and using the subsequent revolution as its background. This
procedure is repeated until a detection is found or until we reach an initial
time of the count rate higher than the 10th of the SNIa model. When this point
is reached, we start again the process with an initial count rate time ranging
from 4th to 6th but with a different revolution as background and so on. An
example output of this method is shown in Figure 6, where the count rate is
computed using the DDTe model with revolution number 1752. The computed
discovery date is at 6.73 day with $\approx$ 67$\sigma$, which is coherent if
we make a visual evaluation of the count rate data.

In some cases, due to high background or too many changes on pointing
direction in the same revolution, the rising of the SN model goes undetected
but on the following revolution, as the SN model flux is higher with time, the
flux is too high to find a good fit. In these cases an inspection of the data
makes clear that there is a high rising signal and the satellite would
recognize it as such. We want to count these cases as late detections,
therefore a condition is created where, if the maximum of the SNIa signal is
larger than 8.5$\%$ of the mean ACS rate, the detection time is saved even if
the fit is not good enough.

The analysis of the count rate has been performed for each of the models (DDTe
and W7) for distances of 4, 6, 8, 10, 12, 14 and 16 kpc. For each distance and
model several pointing directions have been tested to reproduce different
galactic areas backgrounds and to check the variability of the detection on
different positions in the Galaxy. The pointing directions, in galactic
coordinates of longitude (l) and latitude (b), are l = 0, 45, 90 and 180 deg
for b = 0 deg. All the position have a same value of latitude of b = 0 deg.
The reason for this fixed latitude is that the location for these events in
our galaxy would be in regions of old population stars, which are in the old
disk and in a spheroid centered to the galactic center. Therefore, they would
likely occur at $\pm$ 15 deg from the galactic plane. As seen in Figure 3, the
response of the ACS does not change significantly with the direction of the
source (except for low energies). Hence, simplifying the latitude to b = 0 deg
is good enough for our study.

Figure 6: Display of the SN model rate with revolution number 1752 added as
background in green. The background rate, without the SN signal, is shown in
black. The best fitting model to the SN with background rate is shown in red.
For this example, the rising of the SN rate is appreciated and allows the
detection time to be at 6.73 day.

###  4.3 Results

Following the methodology explained on Section 4.2, the time detection for the
W7 and DDTe models placed at distances of 4, 6, 8, 10, 12, 14 and 16 kpc and
at positions l = 0, 45, 90 and 180 deg with b = 0 deg has been computed. An
example of two different distances, 4 kpc and 16 kpc, and l, b = 0 deg, is
shown for models W7 and DDTe in Figure 7, top and bottom respectively. For
both models, the number of detections for higher distances is lower due to the
lower flux.

The mean value of the distribution is computed and shown for all distances and
pointing directions on Figure 8 for the W7 model and Figure 9 for the DDT
model. The time of the early signal detection for each distance is coherent
with the flux $\propto$ distance-2 relation, as for further distances the
detections are later due to lower flux. The W7 model has a slightly earlier
detection for each case because of its faster rise compared to the DDTe one.
The position of the source is not relevant for the early detection of the
model, as it does not change for each longitude tested, see bottom Figure 8
and 9. These results are due to the fact that INTEGRAL observatory has been
pointing in various directions during the $\sim$15 years of observation and
that the equivalent field of view of the ACS/SPI is very large, therefore the
position of a SNIa in the sky does not show any effect.

Figure 7: Top panel shows the distribution of the discovery dates for the W7
model at l = 0 deg and b = 0 deg for distances of 4 kpc (red) and 16 kpc
(black). Bottom panel displays the same for the DDTe model.

Figure 8: Mean discovery dates versus distance for model W7 (top panel) and
longitude (bottom panel).

Figure 9: Mean discovery dates versus distance for model DDTe (top panel) and
longitude (bottom panel).

##  5 Search of SNIa signature in the SPI/ACS data

We performed a systematic search for a SNIa signature in the ACS rate with the
data available at the date of the analysis: 1868 revolutions from revolution
26 to revolution 2063. The analyses were made for each revolution using ACS
data rebinned with a sampling period of 18 s. The detection of a SNIa
signature is triggered when the smoothed residuals333The residuals were
smoothed with a hanning window of 1.5 h to reduce the number of detection due
to gamma-ray bursts and false detection due to statistical fluctuations.
obtained after subtraction of the background model (see section 4.2), are
larger than 5$\sigma$ up to the end of the analyzed revolution. This
criterion, which was chosen for its simplicity, allows the detection of a flux
rising with a long timescale as to the one expected from a SNIa (e.g. see Fig.
6). A first analysis of all data yielded to a lot of detection due to
background variations, which spans up to the end of a revolution, that were
not well modeled with the ACS saturating event rate. Those background
variations are due to low energy particles of the radiation belt or from solar
events. The signature of such particles are clearly seen in the data of the
INTEGRAL Radiation Environment Monitor (IREM), which measures the flux of
particles with energy $>$ 0.5 MeV for electron and 10 MeV for proton (Hajdas
et al., 2003). The counter TC3 of the IREM reproduces well the variations of
ACS rate produced by low energy charged particles of radiation belts. In order
to reduce the number of detections due to low energy particles, we added the
rate of the TC3 counter in the background model for this specific analysis.

With such a detection method, we counted 58 revolutions with SNIa signature.
However, all of them are detected at dates larger than $\sim$ 2.2 days after
the beginning of the revolution. These excesses are explained by an imperfect
modelling of the enhancement of the rate produced by low energy particles
impinging on the scintillators of the ACS when INTEGRAL enters into the
radiation belt at the end of its orbital revolutions (e.g. see Figure 10).

Figure 10: Example of measured and modeled ACS rates (top) and the residuals
(bottom) for the revolution 1342. The imperfect modeling of the ACS rate at
the end of the revolution, when INTEGRAL enters into radiation belts, produces
a smoothed residual excess larger than 5$\sigma$.

##  6 Conclusions

This study demonstrates the capability of the ACS/SPI to detect a galactic
supernova during its early stages. Two theoretical SNIa models with different
luminosities rising times, W7 and DDTe, have been used in order to evaluate
the influence on the time of the detection. For such purpose we have developed
a method to detect the signature of the rising flux emitted by a supernova and
to derive a detection date as quickly as possible. Using simulations of the
count rate produced by the source models and the measured background ACS rates
of 275 orbital revolutions of INTEGRAL, we have found that our method allows
to discover a galactic SN as soon as 6 days up to 12 days after explosion,
depending of the distance and model.

The simulated sources have been placed at several distances (from 4 to 16
kpc), several galactic longitudes and zero latitude. The discovery date
increases with the distance, as the flux of the modelled sources is smaller.
This result follows the distance flux relationship. On the other hand, the
galactic longitude of the source has no influence on the discovery date, which
is logical taking into account the large field of view of the ACS/SPI.

As the discovery date depends on the flux, we have found there is a slightly
different range of discovery dates for the two models. The DDTe model is a
sub-luminous SN Ia model and the W7 a normal luminous one. The lower
luminosity of the DDTe resulted in a detection range of approximately half-day
later than the W7. Using our detection method we have analyzed all the
available ACS/SPI data from the beginning of the INTEGRAL mission to now and
we have not found any signature of a galactic SNIa, suggesting no event has
been missed.

This method is not able to provide a position of a source but it can provide
an early alert to activate other in orbit or ground-based observatories to
quickly browse the sky in search for it. In any case, the collected data would
provide very important information about the early development of the
supernova outburst no matter if the event is detectable in any other energy
range. INTEGRAL is on its last stages of life and the results obtained here
suggest that future missions should consider the possibility of using the ACS
as an all sky detector by improving as much as possible its sensitivity. The
new space telescope COSI 444see website: https://www.nasa.gov/press-
release/nasa-selects-gamma-ray-telescope-to-chart-milky-way-evolution, which
is expected to launch in 2026, will observe the gamma-ray sky in the 200 keV -
5 MeV range (Tomsick et al., 2021). It will be equipped of a compact Compton
telescope made with germanium detectors and of an anticoincidence system made
with BGO scintillators as for SPI. With its large field of view ($\sim$ 25 %
of the sky) and its pointing strategy, COSI will monitor the entire sky within
a day. Both the Compton telescope and the anticoincidence system will be
sensitive to the gamma-ray emission of a galactic SNIa.

## Acknowledgements

The authors thank the anonymous referee for helpful comments. We acknowledge
support from the Spanish Ministry of Science and Innovation and - FEDER UE
(MCI-AEI-FEDER,UE) through grants PID2019-108709GB-I00 (MC,JI) and
PGC2018-095317-B-C21 (EB), by grant 2014 SGR 1458 and CERCA Programe of the
Generalitat de Catalunya (MC,JI), and by the program Unidad de Excelencia
María de Maeztu CEX2020-001058-M (MC,JI).

## Data Availability

The SPI data used for this study are public data available in the ISDC
website: http://www.isdc.unige.ch/integral/archive#DataRelease.

## References

  * Adams et al. (2013) Adams S. M., Kochanek C. S., Beacom J. F., Vagins M. R., Stanek K. Z., 2013, The Astrophysical Journal, 778, 164
  * Borkowski et al. (2004) Borkowski J., Gotz D., Mereghetti S., Mowlavi N., Shaw S., Turler M., 2004, GRB Coordinates Network, 2920, 1
  * Burrows & The (1990) Burrows A., The L., 1990, ApJ, 360, 626
  * Churazov et al. (2014) Churazov E., et al., 2014, Nature, 512, 406
  * Churazov et al. (2015) Churazov E., et al., 2015, ApJ, 812, 62
  * Colgate & McKee (1969) Colgate S. A., McKee C., 1969, ApJ, 157, 623
  * Diehl et al. (2014) Diehl R., et al., 2014, Science, 345, 1162
  * Ferguson et al. (2003) Ferguson C., et al., 2003, A&A, 411, L19
  * Filippenko (1997) Filippenko A. V., 1997, ARA&A, 35, 309
  * Fossey et al. (2014) Fossey S. J., Cooke B., Pollack G., Wilde M., Wright T., 2014, Central Bureau Electronic Telegrams, 3792, 1
  * Gomez-Gomar et al. (1998) Gomez-Gomar J., Isern J., Jean P., 1998, MNRAS, 295, 1
  * Gros et al. (2004) Gros M., et al., 2004, in Schoenfelder V., Lichti G., Winkler C., eds, ESA Special Publication Vol. 552, 5th INTEGRAL Workshop on the INTEGRAL Universe. p. 669
  * Hajdas et al. (2003) Hajdas W., Bühler P., Eggel C., Favre P., Mchedlishvili A., Zehnder A., 2003, A&A, 411, L43
  * Hillebrandt & Niemeyer (2000) Hillebrandt W., Niemeyer J. C., 2000, ARA&A, 38, 191
  * Hoyle & Fowler (1960) Hoyle F., Fowler W. A., 1960, ApJ, 132, 565
  * Isern et al. (2008) Isern J., Bravo E., Hirschmann A., 2008, New Astronomy Reviews, 52, 377
  * Isern et al. (2013) Isern J., et al., 2013, A&A, 552, A97
  * Isern et al. (2014) Isern J., et al., 2014, The Astronomer’s Telegram, 6099, 1
  * Isern et al. (2016) Isern J., et al., 2016, A&A, 588, A67
  * Isern et al. (2021) Isern J., Hernanz M., Bravo E., Grebenev S., Jean P., Renaud M., Siegert T., Vink J., 2021, New Astron. Rev., 92, 101606
  * Jean et al. (1999) Jean P., et al., 1999, Astrophysical Letters and Communications, 38, 421
  * Li & Ma (1983) Li T. P., Ma Y. Q., 1983, ApJ, 272, 317
  * Maoz et al. (2014) Maoz D., Mannucci F., Nelemans G., 2014, ARA&A, 52, 107
  * Mereghetti et al. (2005) Mereghetti S., Götz D., von Kienlin A., Rau A., Lichti G., Weidenspointner G., Jean P., 2005, ApJ, 624, L105
  * Milne et al. (2004) Milne P. A., et al., 2004, ApJ, 613, 1101
  * Nomoto et al. (1984) Nomoto K., Thielemann F. K., Yokoi K., 1984, ApJ, 286, 644
  * Rodríguez-Gasén et al. (2014) Rodríguez-Gasén R., Kiener J., Tatischeff V., Vilmer N., Hamadache C., Klein K. L., 2014, Sol. Phys., 289, 1625
  * Savchenko et al. (2012) Savchenko V., Neronov A., Courvoisier T. J. L., 2012, A&A, 541, A122
  * Savchenko et al. (2017) Savchenko V., et al., 2017, A&A, 603, A46
  * Siegert et al. (2018) Siegert T., et al., 2018, A&A, 615, A107
  * Sturner et al. (2003) Sturner S. J., et al., 2003, A&A, 411, L81
  * Tomsick et al. (2021) Tomsick J. A., et al., 2021, arXiv e-prints, p. arXiv:2109.10403
  * Vedrenne et al. (2003) Vedrenne G., et al., 2003, A&A, 411, L63
  * Von Kienlin et al. (2001) Von Kienlin A., Arend N., Lichti G. G., 2001, in Costa E., Frontera F., Hjorth J., eds, Gamma-ray Bursts in the Afterglow Era. p. 427 (arXiv:astro-ph/0109119), doi:10.1007/10853853_118
  * Von Kienlin et al. (2003) Von Kienlin A., et al., 2003, A&A, 411, L299
  * Wang et al. (2019) Wang X., Fields B. D., Lien A. Y., 2019, MNRAS, 486, 2910


# Significance of the strain-dominated region around a vortex on induced
aerodynamic loads

Karthik Menon1<EMAIL_ADDRESS>Rajat Mittal1<EMAIL_ADDRESS>1Department of
Mechanical Engineering, Johns Hopkins University, Baltimore, MD 21218, USA

###### Abstract

The ability of vortices to induce aerodynamic loads on proximal surfaces plays
a significant role in a wide variety of flows. However, most studies of
vortex-induced effects primarily focus on analyzing the influence of the
rotation-dominated cores of vortices. In this work, we show that not only are
vortices in viscous flows surrounded by strain-dominated regions, but that
these regions are dynamically important and can sometimes even dictate the
induced aerodynamic loads. We demonstrate this for a pitching airfoil, which
exhibits dynamic stall and generates several force-inducing vortices. Using a
data-driven force partitioning method, we quantify the influence of vortices
as well as vortex-associated strain to show that our current understanding of
vortex-dominated phenomena, such as dynamic stall, is incomplete without
considering the substantial effect of strain-dominated regions that are
associated with vortices.

###### keywords:

Vortex dominated flows, Leading-edge vortex, Dynamic stall, Lift generation

##  1 Introduction

Vorticity play a singularly important role in fluid dynamics. This is because
(a) the generation/appearance of vorticity in a flow provides direct
information regarding the influence of viscosity; (b) flow instabilities and
transition to turbulence occur in vorticity dominant regions such as shear
layers, boundary layers and vortices; and finally, (c) vorticity organized in
the form of vortices is known to induce aerodynamic loads on proximal
surfaces,

The current paper concerns this third aspect of vorticity i.e. the ability of
vortices to induce aerodynamic loads on immersed surfaces – a fact that has
been exploited to great effect in nature (Ellington et al., 1996;
Triantafyllou et al., 2000) as well as in engineered wings and control
surfaces (Carr, 1988; Eldredge & Jones, 2019). Indeed, unsteady aerodynamics
is one arena where a number of studies have focused on determining the
aerodynamic loads induced by vortices in a flow (Dickinson & Götz, 1993; Sane,
2003; Shyy & Liu, 2007; Chen et al., 2010; Pitt Ford & Babinsky, 2013; Xia &
Mohseni, 2013; Mulleners & Raffel, 2013; Gharali & Johnson, 2013). This task
is made difficult not only by the fact that most flows consist of multiple
interacting vortices and shear layers, but also that there continues to be
significant ambiguity regarding the very definition of a vortex (Haller, 2015)
and its relationship to the induced pressure field.

One method for vortex identification that is commonly used is the so called
“$Q$-criterion” (Hunt et al., 1988; Jeong & Hussain, 1995). The quantity $Q$
is the second invariant of the velocity gradient tensor and can be defined as
$Q=\frac{1}{2}\left(||\bar{\Omega}||^{2}-||\bar{S}||^{2}\right)$. Here
$\bar{S}$ and $\bar{\Omega}$ are the strain-rate and rotation tensors given by
$\bar{S}=\frac{1}{2}\left[\boldsymbol{\nabla}\boldsymbol{u}+\left(\boldsymbol{\nabla}\boldsymbol{u}\right)^{T}\right]$
and
$\bar{\Omega}=\frac{1}{2}\left[\boldsymbol{\nabla}\boldsymbol{u}-\left(\boldsymbol{\nabla}\boldsymbol{u}\right)^{T}\right]$
respectively. A vortex, according to the $Q$-criterion, is then defined as a
connected region where $Q>0$. This corresponds to regions where the Frobenius
norm of the rotation tensor is greater than that of the strain-rate tensor.

Interestingly, while the quantity $Q$ has mostly been associated with the flow
kinematics, Jeong & Hussain (1995) noted that the pressure Poisson equation
for incompressible flow, derived by taking the divergence of the Navier-Stokes
equation, can be expressed as:

$\nabla^{2}p=-\rho\boldsymbol{\nabla}\cdot\left(\boldsymbol{u}\cdot\boldsymbol{\nabla}\boldsymbol{u}\right)\equiv\rho\boldsymbol{\nabla}\boldsymbol{u}:\left(\boldsymbol{\nabla}\boldsymbol{u}\right)^{T}\equiv
2\rho Q.$ (1)

Equation 1 shows that $Q$ appears in the source term of the pressure Poisson
equation and is thus, not merely a kinematic quantity, but one with direct
implications for the _dynamics_ of the flow as well. Furthermore, while
regions of positive $Q$, which are used to identify vortices, may draw much of
our attention, equation 1 shows that regions with negative $Q$, i.e. strain
dominated regions, are no less significant with regard to their effect on
pressure.

A question that immediately arises from equation 1 is: do vortices have
regions of negative $Q$ associated with them and, if so, how significant are
these regions? To answer this, we first note that the following Neumann
boundary condition completes the prescription of the above boundary value
problem:
$\boldsymbol{\hat{n}}\cdot\boldsymbol{\nabla}p=\boldsymbol{\hat{n}}\cdot\left[-\rho{\frac{d\boldsymbol{u}}{dt}}+\mu\nabla^{2}\boldsymbol{u}\right]$.
We now consider a vortex in an infinite quiescent flow for which, the boundary
condition at the far-field would be homogeneous. The divergence theorem
applied to the pressure Poisson equation would therefore require that
$\int_{V_{f}}Q\ dV\equiv 0$, where $V_{f}$ is the volume of the flow domain.
This implies that for such a vortex, the volume integral of the negative $Q$
region is equal to the volume integral of positive $Q$ in the vortex core.
Furthermore, immediately outside the region of positive $Q$, the velocity
gradients have to be non-zero and therefore, this region immediately
surrounding the vortex core should have $Q<0$. Thus, for any vortex in a
viscous flow, the positive-$Q$ core is necessarily associated with a
negative-$Q$ “corona” of equivalent significance. This can easily be confirmed
for exact vortex solutions of the Navier-Stokes equations such as the Rankine
vortex (Acheson, 1990), the Lamb-Oseen vortex (Lamb, 1993), and others. While
the presence of boundaries and other flow complexities might eliminate the
strict equipartitioning of $Q$ for a vortex, we nevertheless expect that every
vortex will have a significant corona of negative $Q$ around it. This
observation, combined with equation 1, suggests that the strain dominated
region around a vortex could have a significant contribution to the net
aerodynamic load induced by a vortex. However, although the overall effect of
these strain-dominated regions on aerodynamic loading has implicitly been
included in past studies, our physical intuition regarding aerodynamic loading
has generally been tied to the growth, motion, and overall evolution of the
rotational core of vortices. The role of the strain-induced region on force
production has, thus far, not garnered significant attention.

In this paper, we examine the aerodynamic loads induced by vortices on an
immersed surface, and in particular, we focus on the distinct roles played by
regions of positive as well as negative $Q$ associated with vortices. For this
analysis, we choose to focus on the problem of a sinusoidally pitching airfoil
that exhibits dynamic stall, and generates leading-edge (LEV) and trailing-
edge vortices (TEVs). We employ a force partitioning approach that allows us
to precisely quantify the force contributions of various vortices and vortical
regions in the flow. Using this analysis, we examine the phenomenon of dynamic
stall and show that our understanding of vortex-induced loads is incomplete
without also discussing the role played by the corona of negative $Q$ around a
vortex.

##  2 Methods and problem description

###  2.1 Computational model and numerical method

The two-dimensional computational model studied in this work consists of a
sinusoidally pitching NACA0015 airfoil with a rounded trailing edge, immersed
in a uniform free-stream flow with velocity $U_{\infty}$. The flow simulations
have been performed using the sharp-interface immersed boundary method of
Mittal et al. (2008) and Seo & Mittal (2011). This method is well-suited for
moving boundary problems as it allows the use of simple non-conformal
Cartesian grids to simulate a variety of different shapes and motions of the
immersed body. Further, the ability to preserve the sharp-interface around the
immersed boundary ensures very accurate computations of surface quantities.
The Navier-Stokes equations are solved using a fractional-step method. Spatial
derivatives are discretized using second-order central differences, and time-
stepping is achieved using the second-order Adams-Bashforth method. The
pressure Poisson equation is solved using a geometric multigrid method. This
code has been extensively validated in previous studies for a wide variety of
stationary as well as moving boundary problems (Mittal et al., 2008; Seo &
Mittal, 2011), where its ability to maintain local (near the immersed body) as
well as global second-order accuracy has been established.

The Reynolds number of the flow is defined as $Re=U_{\infty}C/\nu$, where $C$
is the chord-length of the airfoil and $\nu$ is the kinematic viscosity of the
fluid. The pitch oscillation is prescribed as:
$\theta=\theta_{0}+A_{\theta}\sin{(2\pi ft)}$ where $\theta$ is the
instantaneous pitch angle, $\theta_{0}$ is the mean pitch angle and
$A_{\theta}$ is the pitch amplitude. The dimensionless pitching frequency is
given by $f=FC/U_{\infty}$ where $F$ is the pitch frequency and $t$ is the
dimensionless time (also scaled by $U_{\infty}$ and $C$). This model is placed
symmetrically in a large $21C\times 23C$ computational domain with a grid of
size $384\times 320$ points. The grid around the airfoil is isotropic, and
corresponds to about $125$ points along the chord. A Dirichlet velocity
boundary condition is used at the inlet boundary of the domain, and zero-
gradient Neumann conditions are specified at all other boundaries. Grid
refinement studies and other details can be found in Menon & Mittal (2019,
2020a).

###  2.2 Force partitioning and vortex-induced aerodynamic loads

Figure 1: (a) Schematic of the setup for the force partitioning method (FPM),
along with relevant symbols. (b) Sample snapshot of the $\phi_{2}$ field
around an airfoil.

A key element of the current study is the ability to estimate the aerodynamic
loads induced by vortices. This is accomplished using the force partitioning
method (FPM) of Zhang et al. (2015) and Menon & Mittal (2020c), a concise
description of which is presented here. The first step in the FPM is the
calculation of an auxiliary potential field, $\phi_{i}$, which satisfies the
following equation:

${\nabla}^{2}\phi_{i}=0,\ \ \mathrm{in}\ \ V_{f}\ \ \mathrm{with}\ \
\boldsymbol{\hat{n}}\cdot\boldsymbol{\nabla}\phi_{i}=\begin{cases}n_{i}\;,\;\mathrm{on}\;B\\\
0\;\;,\;\mathrm{on}\;\Sigma\end{cases}\ \ \ \mathrm{for}\ i=1,2,3.$ (2)

Here, as shown in figure 1(a), $V_{f}$ is the fluid domain with outer boundary
$\Sigma$ and $B$ is the surface of the immersed body. The unit normal pointing
outward from the fluid domain and into the internal and external boundaries,
$B$ and $\Sigma$, is given by $\hat{n}$ and the index $i$ determines the
vector component of the aerodynamic force that is being partitioned. In the
current two-dimensional problem, $i=1$ and $2$ correspond to the drag and lift
forces respectively on the airfoil. We note that this auxiliary potential is a
function only of the instantaneous position and shape of the immersed body and
the outer domain boundary. Figure 1(b) shows an example $\phi_{2}$ field
around an airfoil.

The next step is to project the Navier-Stokes equation onto the gradient of
the auxiliary potential and integrate over the volume of the fluid domain.
This results in,

$-\int_{V_{f}}\boldsymbol{\nabla}p\cdot\boldsymbol{\nabla}\phi_{i}\
dV=\int_{V_{f}}\rho\left[\frac{\partial\boldsymbol{u}}{\partial
t}+\boldsymbol{u}\cdot\boldsymbol{\nabla}\boldsymbol{u}-\nu\nabla^{2}\boldsymbol{u}\right]\cdot\boldsymbol{\nabla}\phi_{i}\
dV\ \mathrm{for}\ i=1,2,3$ (3)

By making use of properties of $\phi_{i}$ as prescribed in equation 2, as well
as the incompressibility constraint and the divergence theorem, the above
expression can be rewritten for the pressure-induced force on the body as
follows:

$\displaystyle F_{i}=\int_{B}p\ n_{i}\ dS=$
$\displaystyle-\rho\int_{B}\boldsymbol{\hat{n}}\cdot\left(\frac{d\boldsymbol{U}_{B}}{dt}\phi_{i}\right)\
dS+\rho\int_{V_{f}}\boldsymbol{\nabla}\cdot\left(\boldsymbol{u}\cdot\boldsymbol{\nabla}\boldsymbol{u}\right)\phi_{i}\
dV$ (4)
$\displaystyle-\mu\int_{V_{f}}\left(\boldsymbol{\nabla}^{2}\boldsymbol{u}\right)\cdot\boldsymbol{\nabla}\phi_{i}\
dV-\rho\int_{\Sigma}\boldsymbol{\hat{n}}\cdot\left(\frac{d\boldsymbol{u}}{dt}\phi_{i}\right)\
dS\ \ \mathrm{for}\ i=1,2,3$

where $\boldsymbol{U}_{B}$ is the local velocity of the immersed surface. The
four terms on the right-hand-side of equation 4 represent distinct components
of the pressure-induced force on the surface. The first term contains the
well-known linear acceleration reaction (i.e. added mass) force (Menon &
Mittal, 2020b) as well as the centripetal acceleration reaction, which was
shown by Zhang et al. (2015) to be an important lift generation mechanism for
insect wings. The third term is the pressure force associated with viscous
diffusion in the flow and this term is usually significant only for low-
Reynolds number flows. The fourth term is associated with flow acceleration at
the outer boundary and this term is negligible for large domains and steady
freestream flows. The surface shear stress also induces a force on the surface
that is independent of the pressure-induced force, but this component of the
force is also very small for high Reynolds number flows. Further details of
this FPM applied to aerodynamic forces as well as moments can be found in
Zhang et al. (2015) and Menon & Mittal (2020b, c).

The second term in the above expression, which is the primary focus of the
current study, can be expressed as,

$F_{i}^{\omega}=\rho\int_{V_{f}}\boldsymbol{\nabla}\cdot\left(\boldsymbol{u}\cdot\boldsymbol{\nabla}\boldsymbol{u}\right)\phi_{i}\
dV\equiv-\rho\int_{V_{f}}2\ Q\ \phi_{i}\ dV\ \ \ \mathrm{for}\ i=1,2,3.$ (5)

Given the centrality of $Q$ in this force component, we identify this term as
the “vortex-induced force” and denote it as $\boldsymbol{F}^{\omega}$. This
vortex-induced force is found to be a dominant component in a variety of flows
including those over flapping wings (Zhang et al., 2015), bluff bodies (Menon
& Mittal, 2020b), pitching foils (Menon & Mittal, 2020c) and delta wings (Li
et al., 2020). This is also true for the particular pitching airfoil problem
relevant to this paper. As we will show in the next section, the order of
magnitude of vortex-induced force oscillations is $O(1)$ and this term
accounts for the majority of the total lift on the airfoil. In comparison, the
viscous and added-mass terms show oscillations with orders of magnitude
$O(10^{-1})$ and $O(10^{-2})$ respectively (not shown here for brevity). Plots
of these force contributions for sample pitching airfoil cases can be found in
Menon & Mittal (2020c).

It is also useful to note that the volume-integral form of the vortex-induced
force term in equation 5 allows us to compute the force contribution of any
region (or flow structure) of interest by constructing appropriate spatial
volumes enclosing the flow structure, and evaluating the integral over this
volume. Zhang et al. (2015) used this idea to determine the contribution to
lift of leading-edge vortices over flapping insect wings and Menon & Mittal
(2020b) employed this to show that during stationary-state, transverse flow-
induced vibration of a bluff body, the wake vortices serve to damp the
vibrations whereas the boundary layer on the oscillating cylinder surface
tends to promote oscillations.

The construction of such integral volumes in complex, time-varying flows is
itself a challenge and in this work, we perform this task using a novel
automated data-driven procedure described in Menon & Mittal (2020c). The input
to this framework is a time-resolved flow-field that is represented on a grid.
At each time-step, regions of strain and rotation-dominated flow are
identified as grid points which satisfy $Q<-\epsilon_{S}|Q_{max}|$ and
$Q>\epsilon_{\Omega}|Q_{max}|$ thresholds respectively, where $|Q_{max}|$ is
representative of the maximum magnitude of $Q$ in the flow. The thresholds
$\epsilon_{S}$ and $\epsilon_{\Omega}$ are chosen to be small numbers,
typically of $O(10^{-3})$, so as to provide an accurate identification of
strain and rotation dominant regions while being relatively insensitive to
numerical errors. In each case, this thresholding yields a subset of the total
grid which consists of several disconnected and dense clusters of grid points.
These grid points are then segmented into distinct flow structures using the
density-based clustering algorithm DBSCAN (Ester et al., 1996). This is
followed by a tracking procedure which allows us to consistently label each
segmented flow structure over multiple time-steps. This tracking is based on
comparing the location of clusters between time-steps using a simple model of
scalar convection. By isolating appropriate spatial volumes corresponding to
these structures as they evolve spatio-temporally, we can then rigorously
evaluate the forces induced on the airfoil by each of these structures as per
equation 5.

For computational feasibility, we perform this isolation, tracking and force
evaluation within a sub-domain of size $2.5C\times 2C$ around the airfoil
within which the computed vortex-induced lift is approximately equal to that
integrated over the entire flow domain. Furthermore, we partition the vortex-
induced lift into components associated with positive $Q$ regions and negative
$Q$ regions; these are referred to as “rotation-induced” and “strain-induced”
lift, respectively, and the corresponding lift coefficients denoted by
$C_{L}^{\Omega}$ and $C_{L}^{S}$, respectively.

Before the results are presented, it is instructive to place this method
within the context of past studies that have quantified the aerodynamic loads
induced by individual vortices in a flow. Two studies that are particularly
relevant here are those of Pitt Ford & Babinsky (2013) and Onoue & Breuer
(2016) who analyzed the loading induced by the leading and trailing-edge
vortices using experimental measurements in unsteady airfoil flows. The common
features of the approaches employed by Pitt Ford & Babinsky (2013) and Onoue &
Breuer (2016) are: (a) the assumption of two-dimensionality; (b) use of
inviscid models with point vortices in the field combined with conformal
mapping and image-vortices to satisfy no-penetration on the airfoil surface;
(c) parameterization of vortex circulation from experimental data; and (d)
calculation of induced aerodynamic load via the Blasius formula (Batchelor,
1967). The predictions from such models can be quite good and lead to
excellent insights into vortex-induced aerodynamic loads.

In comparison, the force partitioning method provides an approach to vortex-
induced load quantification that neither requires the assumption of two-
dimensionality (see Zhang et al. (2015) for application to a 3D flow) nor that
of inviscid flow. This eliminates the need for conformal mapping and the
difficulties associated flow regularization around sharp edges (by enforcing
the Kutta condition). It also removes errors associated with the inviscid flow
assumption and the use of point-vortex models. By alleviating these
constraints, the force partitioning method demonstrated here is therefore
rigorously applicable to general viscous flows. This method allows the
evaluation of aerodynamic loading due to individual flow structures in flows
that contain several such features interacting and deforming in complicated
ways. Furthermore, this method also provides the contribution of acceleration
reaction effects as well as viscous diffusion to aerodynamic loads, and can
therefore also be applied to flows at low Reynolds numbers. Finally, while the
force partitioning is applied here to data from viscous flow simulations, it
is, in principle, applicable to experimental data as well.

##  3 Results

The analysis is based on two cases of a pitching airfoil with different
dimensionless pitching frequencies: $f=0.35$ and $0.10$. The other parameters
which are kept the same for these cases are: $Re=1000$,
$\theta_{0}=15^{\circ}$, and $A_{\theta}=20^{\circ}$. Figure 2 shows the lift
versus pitch angle for the two cases. While both cases exhibit the classic
features of dynamic stall (instantaneous lift that significantly exceeds the
corresponding static lift followed by a rapid drop in lift), the differences
between the two due to the pitching frequency provide an effective substrate
to examine vortex-induced forces in a comprehensive manner.

Figure 2: Coefficient of lift ($C_{L}$) versus pitch-angle ($\theta$) plots
for the two cases: (a) $f=0.35$, $A_{\theta}=20^{\circ}$; (b) $f=0.10$,
$A_{\theta}=20^{\circ}$. Labels $t_{1}$–$t_{4}$ correspond to time-instances
examined in subsequent discussion and figures. The dashed line indicates the
lift coefficient for corresponding static airfoil (Menon & Mittal, 2020a).

###  3.1 Case 1: f=0.35

Figure 3: Flow snapshots and lift contributions for the case with $f=0.35$.
(a) Contours of spanwise vorticity at time $t_{1}$. (b)-(e) Filled contours of
$Q$, with levels $[-40,40]$, overlaid with line contours of $\phi_{2}$.
(f)-(i) Distribution of vortex-induced lift ($f_{2}^{\omega}=-2Q\phi_{2}$)
with contour levels $[-65,65]$. Dashed lines and labels show integration
volumes corresponding to each vortex which are used to compute the vortex-
induced lift associated with these regions as per equation 5. Snapshots
(b)-(e) and (f)-(i) are at time-instances $t_{1}$ – $t_{4}$ from top to
bottom. (j) Total lift coefficient ($C_{L}$;  ) compared with vortex-induced
lift ($F_{2}^{\omega}$;  ). Pitch angle is shown for reference ($\theta$;  ).
(k) Rotation-induced lift ($C^{\Omega}_{L}$;  ) and contributions from the
primary LEV ( ), TEV ( ) and secondary LEV ( ). (l) Rotation-induced lift due
to regions of $Q>0$ ($C^{\Omega}_{L}$;  ) and strain-induced lift due to
regions of $Q<0$ ($C^{S}_{L}$;  ). Time instances $t_{1}$ – $t_{4}$ are
indicated using dashed grey lines.

We first discuss the case with $f=0.35$, and specifically examine the time-
instances marked $t_{1}$–$t_{4}$ in figure 2. Figure 3(a) shows a snapshot of
vorticity at time $t_{1}$, and we see four key features: a large primary LEV
(pLEV) from the previous cycle which is in the near-wake of the foil, a large
TEV which is near the trailing edge, a secondary LEV (sLEV) which is located
near mid-chord and a shear layer on the suction side which is the initial
stage in the formation of the next pLEV.

Figures 3(b)-(e) show the $Q$ fields at time instances $t_{1}$ – $t_{4}$,
overlaid with contours of $\phi_{2}$. These $Q$ fields clearly identify the
key vortices that are generated during the pitching cycle as regions of
positive $Q$ and for time $t_{1}$, figure 3(b) can be correlated with the
vorticity plot in figure 3(a). Also noticeable in the plots of $Q$ are the
large regions of $Q<0$ that surround the vortices, which corroborates our
earlier discussion about the existence of such regions around vortices. These
regions of negative $Q$ tend to be more spread out than the vortex cores.
There is also a region of negative $Q$ near the leading edge of the airfoil
that is sustained throughout the pitching cycle.

The distribution of $\phi_{2}$ has some distinctive features that have
important implications for the vortex-induced lift. First, it has opposite
sign over the two surfaces of the airfoil. Thus a vortex adjacent to one
surface of the airfoil would generate a force in the opposite direction to
that induced by the same vortex adjacent to the other surface. Second, the
magnitude of $\phi_{2}$ is highest near mid-chord and decays towards the
leading and trailing edges. Therefore a vortex located near mid-chord would
generate a larger lift force than the same vortex located near the leading or
trailing-edges. Finally, $\phi_{2}$ decays with distance from the airfoil and
hence, as expected, the influence of a vortex on the lift would diminish with
distance from the airfoil.

Figures 3(f)-(i) show corresponding contour plots of the integrand of the
vortex-induced lift, $f_{2}^{\omega}=-2Q\phi_{2}$. This represents the vortex-
induced lift per unit volume of fluid and these plots can be interpreted based
on the distribution of $Q$ and $\phi_{2}$ shown in figures 3(b)-(e). We see
that the most lift-inducing flow structures are those with large $Q$
magnitudes and/or those near mid-chord due to the large magnitude of
$\phi_{2}$.

The dynamical influence of the various flow structures observed in the flow
snapshots discussed above are quantified in figures 3(j)-(l). Figure 3(j)
shows the time-variation of the total lift coefficient, shown as a solid black
line, and the vortex-induced lift coefficient (dashed black line) over two
pitching cycles. First, it is noted that the vortex-induced lift is nearly
identical to the total lift and this highlights the predominance of this
component of the lift force. As mentioned in section 2.2, the added mass and
viscous contributions have orders of magnitude $O(10^{-2})$ and $O(10^{-1})$
respectively. Second, the lift curve shows a rapid drop between $t_{1}$ and
$t_{3}$, which represents the typical signature of dynamic stall. The
generally accepted notion is that dynamic stall occurs when an attached
leading-edge vortex, which generates high lift, undergoes a “shedding” event
(Carr, 1988; Dickinson & Götz, 1993; Chen et al., 2010). We can evaluate this
notion rigorously by computing the distinct lift contributions of each vortex
highlighted above. This is done using an automated framework that segments and
tracks the volumes of the $Q>0$ regions occupied by each of the vortices.
These spatial volumes are shown as $Q=0$ contours using the dashed lines in
figures 3(f)-(i). The resulting “rotation-induced lift” for the various
vortices is presented in figure 3(k). We see that the pLEV is indeed the
dominant generator of lift on the airfoil, with the sLEV and the TEV
generating much smaller contributions (peak lift due to sLEV and TEV is about
15% of the pLEV). However, very surprisingly, we note that from time $t_{1}$
to $t_{3}$ while the total lift on the airfoil falls, the lift induced by the
pLEV actually increases and reaches a peak at $t_{3}$. This increase in pLEV-
induced lift is consistent with its growth while staying attached in a region
of relatively large $\phi_{2}$, as seen in figures 3(c)-(d) as well as
3(g)-(h). When the pLEV is released and starts to move away from the surface
of the foil, as seen in figures 3(e) and 3(i) at $t_{4}$, the lift induced by
it reduces, as expected. However, the total lift actually increases during
this time, from $t_{3}$ to $t_{4}$. Thus, the total lift on the airfoil is
nearly perfectly out of phase with the lift induced by the pLEV. This
observation runs contrary to our phenomenological understanding of the role of
LEVs in lift enhancement and dynamic stall. In fact, this phase-lag between
the growth of the pLEV and the lift maximum has been previously reported by
Xia & Mohseni (2013) and Gharali & Johnson (2013).

The cause of this seemingly anomalous behavior becomes clear when we examine
the components of the vortex-induced lift associated the negative $Q$ regions
(what we call “strain-induced lift”) as shown in figure 3(l). This plot shows
that even though the positive lift induced by the growing LEV does indeed
increase between phases $t_{2}$ and $t_{3}$, this increase is counteracted by
an increase in the negative lift induced by regions of $Q<0$ that surround the
LEV. These regions of negative $Q$ occur in areas where the $\phi_{2}$ has a
high magnitude, as seen in figures 3(c)-(e), and that helps accentuate their
effect on the total lift. Thus, the rapid drop in lift occurs not due to the
“shedding” of the LEV but due to the growth of the strain-dominated negative
$Q$ regions that develop around the growing LEV. Similarly, the rise in total
lift between $t_{3}$ and $t_{4}$, which occurs despite the reduction in pLEV-
induced lift, is driven by the slightly more rapid reduction in the negative
strain-induced lift. Some studies have in fact, reported a phase lag between
the onset of dynamic stall and attainment of peak LEV circulation (Xia &
Mohseni, 2013; Gharali & Johnson, 2013). For the current case, we show that
such a phase lag, which runs counter to our intuition regarding dynamic stall,
can be explained by the significant influence of strain-induced loading. Thus,
this case challenges our conventional understanding of phenomenological
connection between LEV evolution and dynamic stall, and points to the
significant role that the strain dominated regions can play in determining the
aerodynamic loads on the surface.

###  3.2 Case 2: f=0.10

Figure 4: Flow snapshots and lift contributions for the case with $f=0.10$.
(a) Contours of spanwise vorticity at time $t_{1}$. (b)-(e) Filled contours of
$Q$, with levels $[-40,40]$, overlaid with line contours of $\phi_{2}$.
(f)-(i) Distribution of vortex-induced lift ($f_{2}^{\omega}=-2Q\phi_{2}$)
with contour levels $[-65,65]$. Snapshots (b)-(e) and (f)-(i) are at time-
instances $t_{1}$ – $t_{4}$ from top to bottom. (j) Total lift coefficient
($C_{L}$;  ) compared with vortex-induced lift ($F_{2}^{\omega}$;  ). Pitch
angle is shown for reference ($\theta$;  ). (k) Rotation-induced lift
($C^{\Omega}_{L}$;  ) and contributions from the primary LEV ( ), TEV ( ) and
secondary LEV ( ). (l) Rotation-induced lift due to regions of $Q>0$
($C^{\Omega}_{L}$;  ) and strain-induced lift due to regions of $Q<0$
($C^{S}_{L}$;  ). Time instances $t_{1}$ – $t_{4}$ are indicated using dashed
grey lines.

We now contrast this behaviour with that at low-frequency oscillations of
$f=0.10$. The hysteresis plot in figure 2(b) shows that the pitch-angle
reversal in this case corresponds to a sharper initial drop in lift,
indicating that the vortices in this low frequency oscillation interacts
differently with the timescale of oscillation. The snapshots of $Q$, shown for
this case in figures 4(b)-(e), highlight the main vortex structures in the
flow. We see that it is dominated by the same three vortices as in the
previous case - the primary and secondary LEVs and the TEV. However, on
comparing figures 4(b) and 3(d), which correspond to time instances soon after
the pLEV-separation in the present and previous case respectively, we see that
the pLEV in the current low-frequency case is larger in size and farther away
from the surface of the airfoil. This is a result of the slower pitching
motion, which allows more time for the growth and downstream motion of the
pLEV.

This increased distance of the pLEV from the surface has important
consequences for the lift induced by the pLEV as well as its surrounding
region of strain. Figures 4(f)-(i) show snapshots of $f_{2}^{\omega}$ for this
case, and we see from figures 4(g) and 4(h) that the induced-force density of
the pLEV as well as the straining regions around it are significantly smaller
than in the $f=0.35$ case discussed earlier (note the difference in
$f_{2}^{\omega}$ contour levels between the two cases). In particular, a
comparison of the lift distribution in figures 4(g)-(i) with that of figures
3(g)-(i) shows a stark contrast in the distance of the pLEV from the wall soon
after pLEV-separation, as well as the resultant confinement of the strain
associated with the pLEV near the wall.

The implications of these observations are evident in figures 4(j)-(l). Figure
4(j) again shows that the vortex-induced lift dominates the overall lift
production. The contributions of total rotation-induced lift (i.e. $Q>0$
regions) as well as those of the individual vortices are quantified in figure
4(k). First, we note that the peak in the total rotation-induced lift aligns
(at time $t_{2}$) well with that for the total lift. This is in contrast to
the previous case where the the two were quite out-of-phase. Second, the
contribution of the pLEV peaks at time $t_{1}$, which is in-phase with $C_{L}$
but slightly ahead of the the peak rotation-induced lift. The sLEV appears
around this time and augments the rotation-induced lift, resulting in the peak
$C_{L}^{\Omega}$ at $t_{2}$. There is also a second peak in the total
rotation-induced lift at $t_{4}$ which, as evident from figure 4(k), is
generated by the trailing-edge vortex.

Finally, figure 4(l) shows the contribution of the rotation and strain-induced
lifts to the total vortex-induced lift. We note that the strain-induced lift
has a time-profile that is similar to the rotation-induced lift but remains
consistently lower in magnitude for most of the pitching cycle. This is
consistent with the observation that regions of large-strain around the
vortices are generated relatively further away from the airfoil in this case
as compared to the previous case. Furthermore, this also explains the fact
that $C_{L}$ is in-phase with the rotation-induced lift in this case, unlike
the behaviour seen for $f=0.35$.

Thus, this lower frequency case presents a more conventional phenomenology of
dynamic stall where the total lift correlates, as expected, with the growth
and movement of the rotational cores of the vortices. The strain-induced lift
force, despite being comparable in magnitude to the rotation-induced lift and
of significantly greater magnitude than the total lift, does not play as
significant a role for this case, as it does for the first one.

##  4 Conclusions

We have shown that while studies of vortices and their effects on aerodynamic
loads have mostly focused on the rotational cores of vortices, vortex cores
are in fact surrounded by dynamically-significant regions where strain
dominates over rotation. In particular, these strain-dominated regions exert a
considerable, and sometimes even dominant, influence on the induced
aerodynamic loads. By using simulations of flow past a pitching airfoil and
employing a powerful force partitioning method to separate the force
contributions of various regions of the flow, we have shown that the lift
induced by the strain dominated regions around vortices is of a magnitude
comparable to that induced by the rotational cores of of the vortices. For one
of the two cases presented here, we also show that dynamic stall, i.e. the
rapid drop in lift for an unsteady airfoil, is less connected with the
conventional notion of leading-edge vortex shedding and more with the
formation and evolution of the strain dominated “corona” of the leading edge-
vortices.

We close by pointing out that previous investigations (Pitt Ford & Babinsky,
2013; Xia & Mohseni, 2013; Wang & Eldredge, 2013; Onoue & Breuer, 2016) of
vortex-induced aerodynamic loads do implicitly include the effect of strain-
dominated regions of a vortex in estimating induced loads. However, the
accompanying interpretations of the mechanisms and phenomenology has been
based primarily on the rotational vortex cores and has largely ignored the
strain-dominated regions associated with vortices. The current study showcases
the significant, and sometimes dominant influence of the strain-dominated
regions around a vortex, and highlights the utility of explicitly accounting
for these effects when examining the mechanisms that govern aerodynamic
loading in vortex-dominated flows.

## Acknowledgments

We acknowledge support from AFOSR grant number FA9550-16-1-0404, monitored by
Dr. Gregg Abate, and XSEDE computational resources supported by the NSF
through allocation TG-CTS100002.

## Declaration of interests

The authors report no conflict of interest.

## References

  * Acheson (1990) Acheson, D.J. 1990 Elementary fluid dynamics. Oxford University Press.
  * Batchelor (1967) Batchelor, G K 1967 An introduction to fluid dynamics. Cambridge: Cambridge University Press.
  * Carr (1988) Carr, L. W. 1988 Progress in analysis and prediction of dynamic stall. J. Aircraft 25 (1), 6–17.
  * Chen et al. (2010) Chen, K. K., Colonius, T. & Taira, K. 2010 The leading-edge vortex and quasisteady vortex shedding on an accelerating plate. Phys. Fluids 22 (3), 1–11.
  * Dickinson & Götz (1993) Dickinson, M. H. & Götz, K. G. 1993 Unsteady Aerodynamic Performance of Model Wings At Low Reynolds Numbers. J. Exp. Biol. 174 (1), 45–64.
  * Eldredge & Jones (2019) Eldredge, J. D. & Jones, A. R. 2019 Leading-edge vortices: Mechanics and modeling. Annu. Rev. Fluid Mech. 51, 75–104.
  * Ellington et al. (1996) Ellington, C. P., van den Berg, C., Willmott, A. P. & Thomas, A. L. R. 1996 Leading-edge vortices in insect flight. Nature 384 (6610), 626–630.
  * Ester et al. (1996) Ester, M., Kriegel, H., Sander, J. & Xu, X. 1996 A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. In Proceedings of the 2nd International Conference on Knowledge Discovery and Data Mining, pp. 226–231.
  * Gharali & Johnson (2013) Gharali, K. & Johnson, D. A. 2013 Dynamic stall simulation of a pitching airfoil under unsteady freestream velocity. J. Fluids Struct. 42, 228–244.
  * Haller (2015) Haller, G. 2015 Lagrangian Coherent Structures. Annu. Rev. Fluid Mech. 47 (1), 137–162.
  * Hunt et al. (1988) Hunt, J. C. R., Wray, A. A. & Moin, P. 1988 Eddies, stream, and convergence zones in turbulent flows. Center for turbulence research report CTR-S88 pp. 193–208.
  * Jeong & Hussain (1995) Jeong, J. & Hussain, F. 1995 On the identification of a vortex. J. Fluid Mech. 285, 69–94.
  * Lamb (1993) Lamb, Horace 1993 Hydrodynamics. Cambridge University Press.
  * Li et al. (2020) Li, J., Zhao, X. & Graham, M. 2020 Vortex force maps for three-dimensional unsteady flows with application to a delta wing. J. of Fluid Mech. 900.
  * Menon & Mittal (2019) Menon, K. & Mittal, R. 2019 Flow physics and dynamics of flow-induced pitch oscillations of an airfoil. J. Fluid Mech. 877, 582–613.
  * Menon & Mittal (2020a) Menon, K. & Mittal, R. 2020a Aerodynamic characteristics of canonical airfoils at low Reynolds numbers. AIAA J. 58 (2), 977–980.
  * Menon & Mittal (2020b) Menon, K. & Mittal, R. 2020b On the initiation and sustenance of flow-induced vibration of cylinders: insights from force partitioning. J. Fluid Mech. 907 (A37).
  * Menon & Mittal (2020c) Menon, K. & Mittal, R. 2020c Quantitative analysis of the kinematics and induced aerodynamic loading of individual vortices in vortex-dominated flows: A computation and data-driven approach. arXiv preprint arXiv:2011.04632 .
  * Mittal et al. (2008) Mittal, R., Dong, H., Bozkurttas, M., Najjar, F. M., Vargas, A. & von Loebbecke, A. 2008 A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. J. Comput. Phys. 227 (10), 4825–4852.
  * Mulleners & Raffel (2013) Mulleners, K. & Raffel, M. 2013 Dynamic stall development. Exp. Fluids 54 (2).
  * Onoue & Breuer (2016) Onoue, K. & Breuer, K. S. 2016 Vortex formation and shedding from a cyber-physical pitching plate. J. Fluid Mech. 793, 229–247.
  * Pitt Ford & Babinsky (2013) Pitt Ford, C. W. & Babinsky, H. 2013 Lift and the leading-edge vortex. J. Fluid Mech. 720, 280–313.
  * Sane (2003) Sane, S. P. 2003 The aerodynamics of insect flight. J. Exp. Biol. 206 (23), 4191–4208.
  * Seo & Mittal (2011) Seo, J. H. & Mittal, R. 2011 A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations. J. Comput. Phys. 230 (19), 7347–7363.
  * Shyy & Liu (2007) Shyy, W. & Liu, H. 2007 Flapping wings and aerodynamic lift: The role of leading-edge vortices. AIAA J. 45 (12), 2817–2819.
  * Triantafyllou et al. (2000) Triantafyllou, M. S., Triantafyllou, G. S. & Yue, D. K. P. 2000 Hydrodynamics of Fishlike Swimming. Annu. Rev. Fluid Mech. 32 (1), 33–53.
  * Wang & Eldredge (2013) Wang, C. & Eldredge, J. D. 2013 Low-order phenomenological modeling of leading-edge vortex formation. Theor. Comput. Fluid Dyn. 27 (5), 577–598.
  * Xia & Mohseni (2013) Xia, X. & Mohseni, K. 2013 Lift evaluation of a two-dimensional pitching flat plate. Phys. Fluids 25 (9), 91901\.
  * Zhang et al. (2015) Zhang, C., Hedrick, T. L. & Mittal, R. 2015 Centripetal acceleration reaction: An effective and robust mechanism for flapping flight in insects. PLoS ONE 10 (8), 1–16.