11institutetext: Indian Institute of Astrophysics, Koramangala II Block,
Bangalore 560034, India
22institutetext: Pondicherry University, R.V. Nagar, Kalapet, 605014,
Puducherry, India
33institutetext: Observatoire de Paris, LERMA, College de France, CNRS, PSL
University, Sorbonne University, F-75014 Paris, France

# A Triple AGN in the NGC 7733–7734 Merging Group

Jyoti Yadav jyoti [at] iiap.res.in 0000-0002-5641-8102 1122 Mousumi Das
0000-0001-8996-6474 11 Sudhanshu Barway  0000-0002-3927-5402 11 Francoise
Combes  0000-0003-2658-7893 33

###### Abstract

Context. Galaxy interactions and mergers can lead to supermassive black hole
(SMBH) binaries which become active galactic nuclei (AGN) pairs when the SMBHs
start accreting mass. If there is a third galaxy involved in the interaction,
then a triple AGN system can form.

Aims. Our goal is to investigate the nature of the nuclear emission from the
galaxies in the interacting pair NGC 7733–NGC 7734 using archival VLT/MUSE
Integral field spectrograph data and study its relation to the stellar mass
distribution traced by near-infrared (NIR) observations from the South African
Astronomical Observatory (SAAO).

Methods. We conducted near-infrared observations using the SAAO and identified
the morphological properties of bulges in each galaxy. We used MUSE data to
obtain a set of ionized emission lines from each galaxy and studied the
ionization mechanism. We also examined the relation of the galaxy pair with
any nearby companions with Far-UV observations using the UVIT.

Results.  The emission line analysis from the central regions of NGC 7733 and
NGC 7734 show Seyfert and LINER type AGN activity. The galaxy pair NGC 7733–34
also shows evidence of a third component, which has Seyfert-like emission.
Hence, the galaxy pair NGC 7733–34 forms a triple AGN system. We also detected
an Extended Narrow-line region (ENLR) associated with the nucleus of NGC 7733.

###### Key Words.:

Galaxies: individual: NGC7733, NGC7734 – Galaxies: interactions – Galaxies:
active – Galaxies: Seyfert – Galaxies: star formation – Techniques: imaging
spectroscopy

##  1 Introduction

Galaxy interactions and mergers are the major drivers of galaxy evolution in
our low redshift universe, leading to the growth of supermassive black holes
(SMBHs), bulges, and massive galaxies (Di Matteo et al., 2005). One of the
most favorable environments for such activity are galaxy groups where galaxies
are closely interacting, especially those that have significant reservoirs of
cold gas that can be used to fuel star formation and active galactic nuclear
(AGN) activity (Eastman et al., 2007; Georgakakis et al., 2008; Arnold et al.,
2009; Martini et al., 2009; Tzanavaris et al., 2014).

The tidal forces during galaxy interactions may trigger the formation of bars
and spiral arms that produce gravitational torques on the stars and gas in the
galaxy disks. These, in turn, drive gas into the nuclear regions, triggering
central star formation (Springel et al., 2005; Hopkins & Quataert, 2011). The
gas accretion onto the SMBHs can trigger AGN activity, but the details of how
the gas reaches the inner kpc region are still not clear. However, processes
such as nuclear spirals and kpc-scale bars likely play an important role
(Combes, 2001). Observations suggest that most galaxies are associated with
groups or have some companions (McGee et al., 2009). The galaxy interactions
within the group can strongly affect the galaxy morphologies, and trigger star
formation and AGN activity if the galaxies are gas-rich (Duc & Renaud, 2013).
Gas can also be pulled out due to tidal interaction and ram pressure
stripping, as seen in the tails of HI gas or star-forming knots in H$\alpha$
and ultraviolet (UV) (Kenney et al., 2004; von der Linden et al., 2010; Yadav
et al., 2021).

Figure 1: Upper panel shows the FUV, H$\alpha$ and J band images of galaxy
group NGC 7733–34. Bottom panel shows the zoomed in version of the red
rectangles. Blue box represents the possible merging galaxy 2MASS 2342. Yellow
ellipse shows the bridge between NGC 7733 and merger candidate 2MASS 2342.

As galaxies merge, the dynamical friction on the SMBHs in the individual
galaxies causes them to move closer together. If the SMBHs are accreting, they
can form AGN pairs (Rubinur et al., 2019). Studies show that there is a higher
fraction of dual AGN in interacting galaxies, which indicates that galaxy
interactions can trigger AGN activity (Satyapal et al., 2014; Kocevski et al.,
2015; Koss et al., 2018). Several dual AGN have been detected locally (Koss et
al., 2012; Rubinur et al., 2018; Imanishi et al., 2020), and at higher
redshifts (Myers et al., 2008; Hennawi et al., 2010; Silverman et al., 2020;
Silva et al., 2021). However, only a few binary AGN have been detected (Kharb
et al., 2017).

Galaxy interactions can also lead to triple merger systems, and if the SMBHs
of the individual galaxies are accreting, it will form a triple AGN system.
Although such systems are rare, a few have been discovered (Koss et al., 2012;
Pfeifle et al., 2019). It is important to detect more multiple AGN, in order
to understand how AGN activity can affect the merging processes on kpc scales.
In this letter, we present the detection of a kpc-scale triple AGN system in a
nearby galaxy pair NGC 7733–NGC 7734 (hereafter NGC 7733–34). The instruments,
methods and results are described below. The value of Hubble constant (H0)
used in this paper is 67.8 km s-1 Mpc-1.

Source | R.A | Dec. | z
---|---|---|---
| (hh mm ss) | (dd mm ss) |
NGC 7733 | 23:42:32.79 | -65:57:27.42 |  0.03392$\pm$0.00003
NGC 7734 | 23:42:43.20 | -65:56:44.19 |  0.03527$\pm$0.00006
NGC 7733N | 23:42:32.25 | -65:57:09.43 |  0.03619$\pm$ 0.00012
2MASS 2342 | 23:42:49 | -65:58:26 |
Table 1: Details of the galaxies. Column 2 and 3 are the coordinates of
galaxies in J2000.0. Column 4: redshift of NGC 7733 and NGC 7734 taken from
Mathewson & Ford (1996) and de Vaucouleurs et al. (1991) respectively.

##  2 NGC 7733–34 Group

NGC 7733–34 is a group of interacting galaxies at a distance of 154 Mpc (Arp &
Madore, 1987). NGC7734 is at a distance of 48.3 kpc and has a velocity of 406
km s-1 with respect to NGC7733 (de Vaucouleurs et al., 1991). NGC 7733 is a
barred spiral that shows knots in the arms and hosts a Seyfert 2 nucleus
(Tempel et al., 2016). NGC 7734 is a barred spiral with peculiar arms (de
Vaucouleurs et al., 1991; Paturel et al., 1989)). The atomic gas content of
the galaxy pair is $\leqslant 4.2\times 10^{9}$ M$\odot$ and the molecular
content of NGC 7733 and NGC 7734 is $\leqslant 15.4\times 10^{8}$ M$\odot$ and
$\simeq 1.073\times 10^{10}$ M$\odot$ respectively (Horellou & Booth, 1997).
Jahan-Miri & Khosroshahi (2001) identified the star-forming knots and their
ages in the individual galaxies and found that NGC 7733 has a much younger
stellar population than NGC 7734. The properties of the galaxies are given in
Table 1.

##  3 OBSERVATIONS AND DATA ANALYSIS

###  3.1 UV DATA

We observed NGC 7733–34 in the far-ultraviolet (FUV) band using the
Ultraviolet Imaging Telescope (UVIT) on board ASTROSAT (Kumar et al., 2012).
The UVIT has two telescopes that can simultaneously observe in three bands FUV
(1300–1800 Å), NUV (2000–3000 Å), and the visible band. Drift correction were
done with the visible channel. The UVIT has a spatial resolution of around 1″
and a field of view of 28$\arcmin$. We reduced the UVIT level 1 data
downloaded from the Indian Space Science Data Centre (ISSDC) using CCDLAB
(Postma & Leahy, 2017, 2020).

Figure 2:  [N ii] BPT Diagrams of NGC 7733 and NGC 7734. $\Delta\alpha$ and
$\Delta\delta$ are R.A. and declination with respect to the central position
of the source, which is indicated in Table 1. The solid curve in the BPT is
taken from Kauffmann et al. (2003), and the region below this line shows the
ionization due to star-forming regions. The dashed curve is taken from Kewley
et al. (2001), and the region above this line shows the ionization due to AGN.
The dashed-dotted line is taken from Schawinski et al. (2007), which divides
the LINER and Seyfert AGN. The bins with an amplitude-to-noise (A/N) ratio of
less than 4 in any of the lines are shown in gray.

###  3.2 Near Infrared DATA

We performed deep Near Infrared (NIR) imaging of the galaxy group NGC 7733–34
with the SIRIUS camera (Nagayama et al., 2003) on the Infrared Survey Facility
(IRSF) 1.4 m telescope in Sutherland, South Africa. The SIRIUS camera has
three 1024$\times$1024 HgCdTe detectors, which give simultaneous imaging in
the J, H, and K bands, with a field of view 7.7′ $\times$ 7.7′, and the pixel
scale is 0.45″ (Nagayama, 2012). The images were taken in automatic dithering
mode with $\sim$20″ steps with individual frame exposure times of 30 seconds
each, giving a total exposure time of typically $\sim$120 minutes for the
combined set of exposures. Data reduction was carried out with a pipeline for
the SIRIUS observations, including corrections for non-linearity, dark
subtraction, and flat fielding.

###  3.3 MUSE DATA

We have used Multi Unit Spectroscopic Explorer (MUSE) archival data to study
the optical emission lines from NGC 7733–34. MUSE is an integral field
spectrograph (IFU) on the Very Large Telescope (VLT) and gives 3D
spectroscopic data cubes with very high resolution. We used the data from the
wide-field mode, which has a field of view of 1′ $\times$ 1′. It has a spatial
sampling of 0.2″ $\times$ 0.2″ and a spectral resolution of 1750 at 4650Å to
3750 at 9300Å.

We ran the Galaxy IFU Spectroscopy Tool
(GIST111https://abittner.gitlab.io/thegistpipeline/ version 3.0.3; Bittner et
al. 2019) pipeline to analyze the MUSE data. GIST uses a python implemented
version of GANDALF (Sarzi et al., 2006; Falcón-Barroso et al., 2006; Bittner
et al., 2019) for fitting the gas emission lines and Penalized Pixel-Fitting
(pPXF) (Cappellari & Emsellem, 2004; Cappellari, 2017) method for stellar
continuum fitting. GIST gives stellar kinematics and gas emission line
properties based on binned data over a given wavelength range and signal-to-
noise (S/N). For stellar velocity calculation, the code convolves the linear
combination of spectral energy templates with line of sight velocity
distribution and fits that with the spectra from each bin.

We have Voronoi binned the data based on H$\alpha$ (6555–6575 Å) emission. For
binning, we used an S/N of 30. We have also corrected the spectra for the
Milky Way extinction. We adopted the method used in Comerón et al. (2021) and
fitted the continuum using a multiplicative eighth-order Legendre polynomial.
We used a minimum S/N of 3 to remove noisy data before binning the data. We
have used GIST emission line fluxes for [O iii], [N ii], H$\alpha$ and
H$\beta$ for further analysis.

Fig. 1 shows the FUV, H$\alpha$ and NIR J band images of the galaxy group. The
FUV and J band images reveal an extended source (2MASS 23424898-6558257
referred to as 2MASS 2342) southeast of the galaxy NGC 7733. The FUV image
reveals a bright bridge in the southeastern region connecting NGC 7733 to
2MASS 2342. The extended object in the southeastern region of NGC 7733 could
be interacting with NGC 7733 and forms part of the NGC 7733–34 group. The
projected spatial separation of 2MASS 2342 from NGC 7733 is $\sim$ 87 kpc
assuming that both galaxies are at the same redshift.

##  4 Nature of the emission from the nuclei in NGC 7733–34

To study the mechanism producing the ionized gas, we derived the BPT plots of
the two galaxies (Baldwin et al., 1981). The BPT diagram uses the [O
iii]$\lambda 5007/$H$\beta$ and the [N ii]$\lambda 6563/$H$\alpha$ line ratios
to classify the origin of the emission. Fig. 2 (a) and (b) show the BPT plots
of NGC 7733 and NGC 7734 respectively. The nuclear emission from NGC 7733 lies
in the Seyfert part of the BPT plot, while NGC 7734 shows LINER emission from
the central region. The BPT diagram of NGC 7733 also shows the ongoing star
formation in the arms and Seyfert-like emission along the minor axis of the
bar.

Fig. 3 shows the color composite image of H$\alpha$ (blue), [O iii] (green)
and [N ii] (red) emissions, respectively. The bright H$\alpha$ emission, in
blue along the arms of NGC 7733, confirms the ongoing recent star formation in
the arms of the galaxy. The [O iii] emission is extended towards the
northeast, which confirms that it is due to emission from the extended narrow-
line region (ENLR) along the minor axis of the galaxy. The orientation is due
to projection effects. The ENLR is associated with photoionized gas around the
AGN and can be traced out to $\sim$18 kpc in the northeastern region of NGC
7733. The velocity field and the trailing spirals suggest that the NE side is
the near side of the galaxy, so the [O iii] emission is extended towards us.
The southeastern ENLR is extended away from our line of sight and is not
visible due to obscuration by the dust in the galaxy. Ongoing recent star
formation activity in NGC 7733 (Jahan-Miri & Khosroshahi, 2001) suggests the
presence of gas and dust in this galaxy.

Figure 3: Color composite Image Of NGC 7733. Blue, green, and red represent
the H$\alpha$, [O iii] and [N ii] emissions respectively. The [O iii] emission
towards the northeastern side shows the ENLR.

##  5 Is there a third galaxy?

The FUV, J band, and MUSE H$\alpha$ images show a large bright knot in the
northern arm of NGC 7733. This region was identified as a star-forming knot in
Jahan-Miri & Khosroshahi (2001). The knot appears to be extended in the J band
as shown in the red box in Fig. 1. Fig. 4 shows the velocity field of the
galaxy NGC 7733. Most of the velocity field shows features typical of a
rotating disk, with the Doppler shifted velocities of the approaching (blue)
and receding sides (yellow and red). However, the velocity of the knot (shown
in the rectangle) is clearly different from the velocity field of the disc.

###  5.1 Decomposition using GALFIT

We decomposed the NIR images of NGC 7733, and NGC 7734 into bulge, disc, and
bar components using GALFIT (Peng et al., 2002) 2-D decomposition code. We
modeled both the bulges and the bars using Sersic profiles, and the discs with
exponential profiles. A $\chi^{2}$ minimization technique was used to carry
out this fit. We also used the Sersic profile to fit the knot on the northwest
arm of NGC 7733 using GALFIT. We found that the effective radius re of the
Sersic profile in the $J$ band is 3.204″ which corresponds to $\sim$2.4 kpc.
Such a large value for the radius of the knot suggests that it is likely to be
the bulge of a galaxy, instead of being a star-forming region within NGC 7733.

Figure 4: Stellar velocity map of NGC 7733 obtained from GIST. $\Delta\alpha$
and $\Delta\delta$ are R.A. and Declination with respect to the central
position of the source, which is indicated in Table 1. The black rectangle in
the northwest shows that the velocity of NGC 7733N is different from galactic
rotation.

###  5.2 Emission line Analysis of the Knot in NGC 7733

The knot shows a different velocity compared to the northern arm of NGC 7733
as indicated by the rectangle in Fig. 4. The MUSE spectra from this region
show that there are two sets of emission lines, as shown in Fig. 5. The knot
appears to be redshifted by $\sim$650 km s-1 with respect to the galaxy. This
is a strong indication that the knot is not a part of the galaxy and is
instead a separate system that is visually overlapping with the arm of the
galaxy. We will use ‘NGC 7733N’ to refer to the knot in the following text. To
check the nature of emission from NGC 7733N, we used the fluxes of H$\beta$,
[O iii], H$\alpha$ and [N ii] emission lines for all bins in the rectangle
shown in Fig. 4. We divided the spectrum into three wavelength ranges,
4851–4883, 4996–5029, and 6532–6608 Å. We fitted Gaussian profiles to
H$\beta$, [O iii], H$\alpha$ and [N ii] lines and a constant continuum as
shown in Fig. 5. The emission lines from some of the bins could not be fitted
due to poor S/N. The resultant fits confirmed that NGC 7733N has its own set
of emission lines. We calculated the mean velocity of H$\alpha$ in NGC 7733
and NGC 7733N after fitting the emission lines. We found that NGC 7733N is
moving with a mean velocity of +655$\pm$46 km s-1 with respect to NGC 7733 and
is at a redshift of 0.03619$\pm$ 0.00012.

The BPT diagram of two sets of emission lines is shown in Fig. 6. The red
markers and blue crosses show the emission from NGC 7733N, and NGC 7733,
respectively. The emission from the NGC 7733 component of this region is due
to the star formation, which is expected since it originates from the arm. The
emission from NGC 7733N lies in the Seyfert region, indicating an AGN in this
galaxy.

Figure 5:  The blue and red vertical dashed lines correspond to the spectra of
NGC 7733 and NGC 7733N, respectively. Observed and fitted spectra are shown in
red and black. Continuum is shown by dotted green horizontal lines. The
vertical dashed lines show the peaks of emission lines along with their
respective velocities. Figure 6: BPT diagram for NGC 7733 and its companion
NGC 7733N. Blue crosses and red dots color represents the emission from NGC
7733 and NGC 7733N respectively.

##  6 Is this a triple system?

###  6.1 Tidal features in UV and H$\alpha$

Fig. 1 shows an arm-like structure above the main arm in the northern part of
NGC 7733 in the UVIT FUV and MUSE H$\alpha$ images. The arm or tidal tail
starts from the western side and merges with the location of NGC 7733N. It is
likely to have formed via the tidal interaction between NGC 7733 and NGC
7733N. There are UV and H$\alpha$ bright knots in the tidal arm, which
suggests that there is ongoing massive star formation triggered by the
interaction. NGC 7733N then represents the bulge of the merging companion
galaxy (Fig. 1), which may have already lost most of its disk mass in the
tidal interaction.

Galaxy |  MK | Re |  log(M⋆) |  Black Hole Mass (108 M⊙)
---|---|---|---|---
| (mag) | (arcsec) |  |  M∙–LK |  M∙–$\sigma$ |  M∙-Mbulge
|  |  |  | MH2003 | KH2013 | Gu2009 | MM2013 | KH2013 | MH2003 | MM2013 | KH2013
NGC 7733 | -24.09 | 1.73 | 11.39 |  1.83${}^{+0.62}_{-0.43}$ |  4.68${}^{+1.14}_{-1.03}$ |  0.54${}^{+0.05}_{-0.05}$ |  0.64${}^{+0.03}_{-0.03}$ |  1.24${}^{+0.06}_{-0.07}$ |  1.31${}^{+0.89}_{-0.31}$ |  1.86${}^{+0.73}_{-0.56}$ |  3.04${}^{+0.91}_{-0.73}$
NGC 7734 | -23.74 | 1.58 | 11.81 |  1.28${}^{+0.36}_{-0.26}$ |  3.18${}^{+0.90}_{-0.77}$ |  1.13${}^{+0.20}_{-0.18}$ |  1.70${}^{+0.18}_{-0.17}$ |  2.63${}^{+0.28}_{-0.26}$ |  0.90${}^{+0.60}_{-0.21}$ |  1.25${}^{+0.57}_{-0.41}$ |  1.98${}^{+0.67}_{-0.52}$
NGC 7733N | -22.11 | 3.1 | 9.96 |  0.23${}^{+0.01}_{-0.01}$ |  0.51${}^{+0.24}_{-0.17}$ |  |  |  |  0.15${}^{+0.10}_{-0.04}$ |  0.19${}^{+0.14}_{-0.09}$ |  0.26${}^{+0.14}_{-0.09}$
Table 2: Properties of Galactic Bulges and central black holes. Column 2: K
band Absolute magnitude of bulges. Column 3: Effective radius of bulges.
Column 4: Stellar mass estimated in K band using $M/L=0.6$ (McGaugh &
Schombert, 2014). Column 5 and 6: Mass of central black holes using M∙ and K
band luminosity. Column 7, 8 and 9: Mass of central black holes using
M∙–$\sigma$ relation. Column 10, 11 and 12: Mass of central black holes using
M∙– Mass of bulge relation. MH2003: Marconi & Hunt (2003), KH2013: Kormendy &
Ho (2013), Gu2009: Gültekin et al. (2009), MM2013: McConnell & Ma (2013).

###  6.2 Masses of the central Black holes in NGC 7733, NGC 7733N and NGC 7734

We estimated the K band magnitude and effective radius of the bulges in NGC
7733, NGC 7733N, and NGC 7734 using GALFIT and corrected the magnitude for
Galactic extinction using Schlafly & Finkbeiner (2011). We estimated the
stellar mass (M⋆) of each galaxy using Mass–Luminosity relation (McGaugh &
Schombert, 2014). Based on the stellar mass ratios, NGC 7733-NGC 7734 is a
major while NGC 7733N-NGC 7733 is a minor merging system, respectively. The
nuclear black hole masses (M∙) were estimated using the following relations;
(i) M∙–luminosity in the K band (Marconi & Hunt, 2003; Kormendy & Ho, 2013),
(ii) the M–$\sigma$ relation (Gültekin et al., 2009; McConnell & Ma, 2013;
Kormendy & Ho, 2013) and the (iii) M∙–bulge mass relation (Marconi & Hunt,
2003; McConnell & Ma, 2013; Kormendy & Ho, 2013). The SMBH masses for the
three AGN and the stellar masses of the galaxies are listed in Table 2. The
bulge velocity dispersion of NGC 7733N is not reliable, so we could not
calculate the black hole mass in NGC 7733N using the M–$\sigma$ relation. The
nuclear black hole masses of NGC7733 (0.54–4.68 $\times 10^{8}$ M⊙) and
NGC7734 (0.90–3.18 $\times 10^{8}$ M⊙) have a similar range of values, and
their stellar masses are also comparable. NGC 7733N hosts the least massive
stellar disk and the least massive nuclear black hole (1.53–5.15 $\times
10^{7}$ M⊙) in the triple AGN system.

##  7 Implications for our triple AGN detection

Dual AGN have been detected in many galaxies, but triple AGN are rare. Our
study shows that multiple AGN systems can be present in galaxy groups that
show ongoing mergers. Since galaxy groups are fairly common (Tempel et al.,
2017), they can be good targets for detecting multiple AGN. Another important
implication of our study is the enhanced AGN feedback in such systems, which
can heat the intra-cluster medium in clusters and enrich the CGM around
galaxies (Fabian, 2012; Tumlinson et al., 2017). Although multiple AGN
feedback is relatively unexplored, it must play an important role in building
hot gas reservoirs around galaxies, especially at early epochs when galaxy
mergers were more frequent (Lackner et al., 2014; Man et al., 2016; Mundy et
al., 2017). Feedback processes such as jets or intense quasar outflows from
AGN can affect the evolution of the host galaxy (Salomé et al., 2016; Harrison
et al., 2018). Multiple AGN enhance feedback effects may lead to faster SMBH
growth (Volonteri et al., 2020). There may be additional star formation if the
AGN outflows interact with each other. Triple AGN systems such as the one we
have detected in NGC 7733–34 are thus ideal laboratories to study the growth
of hot gas, AGN feedback, and SMBH growth in mergers and should be relatively
common in gas-rich galaxy groups.

##  8 Conclusions

  1. 1.

We find that the paired galaxies NGC 7733 and NGC 7734 host Liner and Seyfert
type nuclei, respectively.

  2. 2.

We also detected an ENLR from the AGN in NGC 7733 towards the northeastern
region.

  3. 3.

The multiwavelength data shows a tidal arm above the northern arm of NGC 7733,
and the tidal arm connects with the location of a third galaxy which we call
NGC 7733N.

  4. 4.

We have confirmed the existence of the third galaxy NGC 7733N in the NGC
7733–34 group. It appears to overlap with the northern arm of NGC 7733. The
emission line analysis of NGC 7733N shows that it hosts a Seyfert nucleus.

  5. 5.

The UV and NIR images show a fourth galaxy (2MASS 2342) lying southeast of the
galaxy pair, which could be part of the merging system. We see a potential
bridge between NGC 7733 and 2MASS 2342 in the UV map. This bridge is composed
of hot gas that has been pulled out during the interaction phase. However,
redshift estimation is necessary to confirm this.

Triple-AGN are rare, however, detailed studies could enhance their observed
frequency in groups, and especially at high redshift. Although this study
focuses only on one system, the results suggest that small merging groups are
ideal laboratories to study multiple AGN systems.

###### Acknowledgements.

We thank the anonymous referee for the thoughtful review, which improved the
impact and clarity of this work. This publication uses the data from the UVIT,
which is part of the AstroSat mission of the Indian Space Research
Organisation (ISRO), archived at the Indian Space Science Data Centre (ISSDC).
This research has also used NIR data from IRSF at the South African
Astronomical Observatory (SAAO). Part of the results are based on observations
collected at the European Southern Observatory under ESO programme
0103.A-0637, run B. This research has made use of the NASA/IPAC Extragalactic
Database (NED) which is operated by the Jet Propulsion Laboratory, California
Institute of Technology, under contract with the National Aeronautics and
Space Administration. M.D. acknowledges the support of Science and Engineering
Research Board (SERB) MATRICS grant MTR/2020/000266 for this research.

## References

  * Arnold et al. (2009) Arnold, T. J., Martini, P., Mulchaey, J. S., Berti, A., & Jeltema, T. E. 2009, ApJ, 707, 1691
  * Arp & Madore (1987) Arp, H. C. & Madore, B. 1987, A catalogue of southern peculiar galaxies and associations
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Bittner et al. (2019) Bittner, A., Falcón-Barroso, J., Nedelchev, B., et al. 2019, A&A, 628, A117
  * Cappellari (2017) Cappellari, M. 2017, MNRAS, 466, 798
  * Cappellari & Emsellem (2004) Cappellari, M. & Emsellem, E. 2004, PASP, 116, 138
  * Combes (2001) Combes, F. 2001, in Advanced Lectures on the Starburst-AGN, ed. I. Aretxaga, D. Kunth, & R. Mújica, 223
  * Comerón et al. (2021) Comerón, S., Knapen, J. H., Ramos Almeida, C., & Watkins, A. E. 2021, A&A, 645, A130
  * de Vaucouleurs et al. (1991) de Vaucouleurs, G., de Vaucouleurs, A., Corwin, Herold G., J., et al. 1991, Third Reference Catalogue of Bright Galaxies
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604
  * Duc & Renaud (2013) Duc, P.-A. & Renaud, F. 2013, Tides in Colliding Galaxies, ed. J. Souchay, S. Mathis, & T. Tokieda, Vol. 861, 327
  * Eastman et al. (2007) Eastman, J., Martini, P., Sivakoff, G., et al. 2007, ApJ, 664, L9
  * Fabian (2012) Fabian, A. C. 2012, ARA&A, 50, 455
  * Falcón-Barroso et al. (2006) Falcón-Barroso, J., Bacon, R., Bureau, M., et al. 2006, MNRAS, 369, 529
  * Georgakakis et al. (2008) Georgakakis, A., Gerke, B. F., Nandra, K., et al. 2008, MNRAS, 391, 183
  * Gültekin et al. (2009) Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, ApJ, 698, 198
  * Harrison et al. (2018) Harrison, C. M., Costa, T., Tadhunter, C. N., et al. 2018, Nature Astronomy, 2, 198
  * Hennawi et al. (2010) Hennawi, J. F., Myers, A. D., Shen, Y., et al. 2010, ApJ, 719, 1672
  * Hopkins & Quataert (2011) Hopkins, P. F. & Quataert, E. 2011, MNRAS, 415, 1027
  * Horellou & Booth (1997) Horellou, C. & Booth, R. 1997, A&AS, 126, 3
  * Imanishi et al. (2020) Imanishi, M., Kawamuro, T., Kikuta, S., Nakano, S., & Saito, Y. 2020, ApJ, 891, 140
  * Jahan-Miri & Khosroshahi (2001) Jahan-Miri, M. & Khosroshahi, H. G. 2001, AJ, 122, 792
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055
  * Kenney et al. (2004) Kenney, J. D. P., van Gorkom, J. H., & Vollmer, B. 2004, AJ, 127, 3361
  * Kewley et al. (2001) Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121
  * Kharb et al. (2017) Kharb, P., Lal, D. V., & Merritt, D. 2017, Nature Astronomy, 1, 727
  * Kocevski et al. (2015) Kocevski, D. D., Brightman, M., Koekemoer, A. M., & Nandra, K. 2015, Are Compton-Thick AGN the Missing Link Between Mergers and Black Hole Growth?, HST Proposal
  * Kormendy & Ho (2013) Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 511
  * Koss et al. (2012) Koss, M., Mushotzky, R., Treister, E., et al. 2012, ApJ, 746, L22
  * Koss et al. (2018) Koss, M. J., Blecha, L., Bernhard, P., et al. 2018, Nature, 563, 214
  * Kumar et al. (2012) Kumar, A., Ghosh, S. K., Hutchings, J., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8443, Space Telescopes and Instrumentation 2012: Ultraviolet to Gamma Ray, ed. T. Takahashi, S. S. Murray, & J.-W. A. den Herder, 84431N
  * Lackner et al. (2014) Lackner, C. N., Silverman, J. D., Salvato, M., et al. 2014, AJ, 148, 137
  * Man et al. (2016) Man, A. W. S., Zirm, A. W., & Toft, S. 2016, ApJ, 830, 89
  * Marconi & Hunt (2003) Marconi, A. & Hunt, L. K. 2003, ApJ, 589, L21
  * Martini et al. (2009) Martini, P., Sivakoff, G. R., & Mulchaey, J. S. 2009, ApJ, 701, 66
  * Mathewson & Ford (1996) Mathewson, D. S. & Ford, V. L. 1996, ApJS, 107, 97
  * McConnell & Ma (2013) McConnell, N. J. & Ma, C.-P. 2013, ApJ, 764, 184
  * McGaugh & Schombert (2014) McGaugh, S. S. & Schombert, J. M. 2014, AJ, 148, 77
  * McGee et al. (2009) McGee, S. L., Balogh, M. L., Bower, R. G., Font, A. S., & McCarthy, I. G. 2009, MNRAS, 400, 937
  * Mundy et al. (2017) Mundy, C. J., Conselice, C. J., Duncan, K. J., et al. 2017, MNRAS, 470, 3507
  * Myers et al. (2008) Myers, A. D., Richards, G. T., Brunner, R. J., et al. 2008, ApJ, 678, 635
  * Nagayama (2012) Nagayama, T. 2012, African Skies, 16, 98
  * Nagayama et al. (2003) Nagayama, T., Nagashima, C., Nakajima, Y., et al. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, ed. M. Iye & A. F. M. Moorwood, 459–464
  * Paturel et al. (1989) Paturel, G., Fouque, P., Bottinelli, L., & Gouguenheim, L. 1989, A&AS, 80, 299
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266
  * Pfeifle et al. (2019) Pfeifle, R. W., Satyapal, S., Manzano-King, C., et al. 2019, ApJ, 883, 167
  * Postma & Leahy (2017) Postma, J. E. & Leahy, D. 2017, PASP, 129, 115002
  * Postma & Leahy (2020) Postma, J. E. & Leahy, D. 2020, PASP, 132, 054503
  * Rubinur et al. (2018) Rubinur, K., Das, M., & Kharb, P. 2018, Journal of Astrophysics and Astronomy, 39, 8
  * Rubinur et al. (2019) Rubinur, K., Das, M., & Kharb, P. 2019, MNRAS, 484, 4933
  * Salomé et al. (2016) Salomé, Q., Salomé, P., Combes, F., & Hamer, S. 2016, A&A, 595, A65
  * Sarzi et al. (2006) Sarzi, M., Falcón-Barroso, J., Davies, R. L., et al. 2006, MNRAS, 366, 1151
  * Satyapal et al. (2014) Satyapal, S., Ellison, S. L., McAlpine, W., et al. 2014, MNRAS, 441, 1297
  * Schawinski et al. (2007) Schawinski, K., Thomas, D., Sarzi, M., et al. 2007, MNRAS, 382, 1415
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Silva et al. (2021) Silva, A., Marchesini, D., Silverman, J. D., et al. 2021, ApJ, 909, 124
  * Silverman et al. (2020) Silverman, J. D., Tang, S., Lee, K.-G., et al. 2020, ApJ, 899, 154
  * Springel et al. (2005) Springel, V., Di Matteo, T., & Hernquist, L. 2005, ApJ, 620, L79
  * Tempel et al. (2016) Tempel, E., Kipper, R., Tamm, A., et al. 2016, A&A, 588, A14
  * Tempel et al. (2017) Tempel, E., Tuvikene, T., Kipper, R., & Libeskind, N. I. 2017, A&A, 602, A100
  * Tumlinson et al. (2017) Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389
  * Tzanavaris et al. (2014) Tzanavaris, P., Gallagher, S. C., Hornschemeier, A. E., et al. 2014, ApJS, 212, 9
  * Volonteri et al. (2020) Volonteri, M., Pfister, H., Beckmann, R. S., et al. 2020, MNRAS, 498, 2219
  * von der Linden et al. (2010) von der Linden, A., Wild, V., Kauffmann, G., White, S. D. M., & Weinmann, S. 2010, MNRAS, 404, 1231
  * Yadav et al. (2021) Yadav, J., Das, M., Patra, N. N., et al. 2021, ApJ, 914, 54


# Fully Implicit Spectral Boundary Integral Computation of Red Blood Cell Flow

Pei Chuan Chao
Department of Civil, Structural and Environmental Engineering
University at Buffalo
212 Ketter Hall, Buffalo, NY 14260
<EMAIL_ADDRESS>
& Ali Gürbüz
Department of Civil, Structural and Environmental Engineering
University at Buffalo
212 Ketter Hall, Buffalo, NY 14260
<EMAIL_ADDRESS>
& Frederick Sachs
Department of Physiology and Biophysics
University at Buffalo
304 Cary Hall, Buffalo, NY 14214
<EMAIL_ADDRESS>
& M. V. Sivaselvan
Department of Civil, Structural and Environmental Engineering
University at Buffalo
212 Ketter Hall, Buffalo, NY 14260
<EMAIL_ADDRESS>
Corresponding author

###### Abstract

An approach is presented for implicit time integration in computations of red
blood cell flow by a spectral boundary integral method. The flow of a red cell
in ambient fluid is represented as a boundary integral equation (BIE), whose
structure is that of an implicit ordinary differential equation (IODE). The
cell configuration and velocity field are discretized with spherical
harmonics. The IODE is integrated in time using a multi-step implicit method
based on backward difference formulas, with variable order and adaptive time
stepping controlled by local truncation error and convergence of Newton
iterations. Jacobians of the IODE, required for Newton’s method, are
implemented as Jacobian matrix-vector products that are nothing but
directional derivatives. Their computation is facilitated by the weakly
singular format of the BIE, and these matrix-vector products themselves amount
to computing a second BIE. Numerical examples show that larger time steps are
possible, and that the number of matrix-vector products is comparable to
explicit methods.

_Keywords_ Weakly singular boundary integral $\cdot$ Fully implicit time
integration $\cdot$ Adaptive time stepping $\cdot$ Spherical harmonics $\cdot$
Elastic thin-shell mechanics $\cdot$ Jacobian-vector multiplication

##  1 Introduction

Cell-level computational modeling of blood flow has been of great interest,
because such simulations offer insight into the rheological character of blood
in microcapillaries, aid in design of microfluidic devices for research and
diagnosis, and have potential to guide targeted drug delivery and treatment
[17, 30, 9]. Our interest in modeling red blood cell (RBC) flow stems from
seeking to explore and test hypotheses about the role of the PIEZO1
mechanosensitive ion channel (MSC) in normal functioning of RBC; evidence of
such a role is offered by clinical studies and by electrophysiological and
microfluidic experiments [3, 34, 4, 10, 12]. This short note describes a
byproduct of this exploration, an implicit time integration method for
simulation of RBC flow. Investigating the role of PIEZOs involves simulation
of RBC flow in microcapillaries of diameter less than the RBC size [19]. In
this note, however, the simpler computational problem of deformation of a
single RBC in ambient flow in an unbounded domain is considered, to present
the formulation and implementation of the implicit method.

RBC flow at the cellular level is typically modeled using Stokes equation to
represent the intra- and extracellular plasma, with the cell membrane
described by elastic thin shell mechanics. Within this modeling framework,
different mathematical formulations and corresponding numerical methods have
been adopted (see [17] for a comprehensive review). Among these, a commonly
employed approach is to formulate the Stokes flow of the intra- and
extracellular plasma as boundary integral equations (BIE) [26, 24, 25]. Since
RBC flow involves moving boundaries corresponding to the cell membranes, the
BIE formulation, where different fields are computed only on boundaries, and
not in the interiors, is beneficial. Further within the BIE setting, different
schemes have been employed to discretize fields on the cell membrane, such as
finite elements [14, 5] and Fourier-like orthogonal spherical harmonic basis
functions [35, 33]. In this short note, we use a BIE formulation and spherical
harmonic basis functions to discretize the RBC deformed geometry and the
velocity field on the cell membrane. The cell membrane is nearly inextensible,
and this is represented in elastic shell models either as a direct constraint
[33, 28], or by a penalty approach by means of a sufficiently large membrane
dilatation modulus [35, 16]. In this note, the latter approach is adopted.

The BIE for the cell membrane velocity, together with the membrane shell
mechanics, constitutes an implicit ordinary differential equation (IODE) [2]
in the membrane configuration (see equation (11))111When a vessel is present
as well in the model, it is an index-1 differential algebraic equation (DAE)
[2], with the traction and velocity fields on the vessel wall being algebraic
variables; if the vessel were also deformable, the BIE system would again be
an IODE. This is the case when cell membrane inextensibility is modeled using
a penalty approach; when inextensibility is imposed as a constraint as in [33,
28], then the BIE is a DAE.. Most commonly, this IODE is integrated in time
using explicit schemes such as forward Euler [29, 24, 25, 35, 36]. Explicit
schemes have been successfully used even in challenging scenarios, for example
simulating red cells passing through very small capillaries and slits, and
taking on complex configurations consisting of high buckled modes of
deformation [35, 16]. However, it has been demonstrated in [33] that the BIE
IODE/DAE system is stiff due to inextensibility of the membrane and higher
order shape derivatives that occur in computing bending deformation. Hence
when explicit schemes are used, time step is restricted by stability rather
than by accuracy. Furthermore, modeling the kinetics of the PIEZO1 MSC will
introduce additional time scales [3, 4]. Therefore, robustness with respect to
time step size is desirable, and we seek to implement an implicit integration
scheme.

Implicit methods for time integration of BIE arising in interface evolution
problems such as RBC flow have been explored less in the literature.
Dimitrakopoulos [13] presented an implicit method to compute the evolution of
a three-dimensional droplet in ambient fluid of a different viscosity. The
interface was modeled as having constant surface tension (unlike more general
elastic shell mechanics in [35, 16, 33, 28] and here). Special characteristics
of the droplet mechanics and a perturbation approach were used to reduce the
nonlinear algebraic equation in each time step to one in a scalar amplitude
field of a prescribed vector field search direction. A single-step implicit
Runge-Kutta method and a multi-step method based on backward difference
formulas (BDF) were both considered for time discretization. Veerapaneni et al
[33], building on their work on axisymmetric flows [32], proposed a semi-
implicit scheme for simulation of three-dimensional vesicle flows. The fluids
inside and the outside the vesicles had the same viscosity, and so the BIE
contained only the single layer potential. Inextensibility of the membrane was
modeled as a constraint, and the membrane mechanics was derived from a strain
energy function of mean curvature. Backward Euler discretization was used, and
the semi-implicit (or linearly implicit) scheme consisted of splitting the
computation of the bending component of the membrane force in such a way that
the new configuration can be computed by solving a linear equation. Rahimian
et al [28] extended this work to account for different viscosities inside and
outside the vesicle, so that the BIE also included a double layer potential.
The BIE was treated using a Galerkin approach (as also done in this note). The
implicit-explicit split was again applied to the bending force, to obtain a
linearly implicit scheme. A hierarchy of iterative methods were considered for
linear system solution, and their performances compared. Despite reduced
overall computational cost compared to the explicit scheme, large numbers of
iterations per time step were needed. Quaife and Biros [27] presented an
adaptive semi-implicit method for two-dimensional vesicle suspensions. By
adopting a spectral deferred correction method, the simulation was allowed to
automatically select time step size and method accuracy up to third order.

In this note, an implicit method is presented that is quite simple to
implement. The flow is formulated as a BIE and the membrane mechanics using
elastic thin shell theory, with in-plane inextensibility of the membrane
captured by a penalty approach. Spherical harmonic basis functions are used to
discretize different fields on the cell membrane, and as in [35], the package
SPHEREPACK [1] is used to perform computations related to spherical harmonics.
The BIE, including the stress kernel, is set up in a _weakly singular_ format
following [22], facilitating its linearization. The IODE structure of the BIE
is deliberately recognized, and an off-the-shelf IODE/DAE solver [21] is used.
The solver is passed two problem-specific routines — one that evaluates the
BIE producing a residual, and another that computes Jacobian matrix-vector
products. Computing a Jacobian matrix-vector product amounts to evaluating
another BIE (see equation (13)), and the matrix itself is not constructed. In
forming the Jacobian matrix-vector product, both the BIE kernels and the
membrane force are linearized. The IODE solver uses BDF for time
discretization, and adapts time step as well as BDF order based on truncation
error. The solver also balances truncation error and accuracy of solution of
the nonlinear algebraic equations in each time step, resulting in a comparable
number of Jacobian matrix-vector products per time step as needed in an
iterative solution of the linear system arising in explicit methods.

This short note is organized as follows. The BIE for RBC flow is presented in
Section 2 in weakly singular form. Section 3 summarizes key results from shell
mechanics. The IODE structure of the BIE is discussed in Section 4, which
concludes by stating the two functions — BIE evaluation and Jacobian matrix-
vector product computation, to be provided to the IODE solver. The BIE
directional derivatives needed for implicit solution are presented in Section
5. Numerical examples are presented in Section 6, highlighting various
characteristics of the implicit solution. Finally some key observations and
potential for further development of this implicit scheme implementation are
summarized in Section 7.

##  2 Boundary integral equation (BIE)

Figure 1: Reference and deformed geometries of the RBC, the shapes of which
are isomorphisms of the unit sphere (parent geometry); also shown is the
coordinate system on the unit sphere. The hatted versions of the functions are
the ones approximated by spherical harmonic basis functions

With the terminology and notation defined in Figure 1, the velocity field,
$\bar{u}$, on the surface, $\gamma$, of a red cell in a three-dimensional
ambient flow, $\bar{u}^{\infty}$, satisfies the following BIE

$\frac{\lambda+1}{2}\bar{u}(\chi)+\frac{\lambda-1}{8\pi}\\!\int_{\gamma}K(\chi,\xi)\bar{u}(\xi)\operatorname{d\\!}\gamma\\\
=\bar{u}^{\infty}(\chi)-\frac{1}{8\pi\mu}\\!\int_{\gamma}G(\chi,\xi)\bar{f}(\xi)\operatorname{d\\!}\gamma\;\;\forall\chi\\!\in\\!\gamma$
(1)

where $\mu$ is the viscosity of the extracellular plasma, $\lambda\mu$ that of
the intracellular plasma, $\bar{f}$ is the net force on the membrane applied
by the fluids on either side of it per unit area of the surface $\gamma$, and
the kernels,

$G(\chi,\xi)=\frac{1}{\mathinner{\\!\left\lVert
r\right\rVert}}(I+\tilde{r}\tilde{r}^{\top});\quad
K(\chi,\xi)=\frac{6}{\mathinner{\\!\left\lVert
r\right\rVert}^{2}}(\tilde{r}^{\top}\tilde{n})(\tilde{r}\tilde{r}^{\top})$ (2)

are fundamental solutions for Stokes flow, with $r=\chi-\xi$ and the unit
vector, $\tilde{r}=\frac{r}{\mathinner{\\!\left\lVert r\right\rVert}}$,
$\tilde{n}$ being the unit normal vector field on $\gamma$ pointing outward to
ambient fluid. Equation (1) is the same as equation (11) in [24] (see also
equation (2) on page 143 of [26]), equation (3a) in [28], and equation (6) in
[35], with appropriate mapping of notation.

Since the kernel $G$ contains $\mathinner{\\!\left\lVert r\right\rVert}$ in
the denominator, the integral on the right hand side of equation (1) is weakly
singular, whereas the integral on the left side is strongly singular because
the kernel $K$ contains $\mathinner{\\!\left\lVert r\right\rVert}^{2}$ in the
denominator [26]. However, a weakly singular representation of the integral on
the left may be obtained by _singularity subtraction_ [22]; subtracting and
adding
$\left[\int_{\gamma}K(\chi,\xi)\operatorname{d\\!}\gamma\right]\\!\bar{u}(\chi)$,
and using the identity,
$\int_{\gamma}K(\chi,\xi)\operatorname{d\\!}\gamma\\!=\\!-4\pi I$ [26, 22], we
obtain

$\bar{u}(\chi)+\frac{\lambda-1}{8\pi}\\!\int_{\gamma}K(\chi,\xi)(\bar{u}(\xi)-\bar{u}(\chi))\operatorname{d\\!}\gamma=\bar{u}^{\infty}(\chi)-\frac{1}{8\pi\mu}\\!\int_{\gamma}G(\chi,\xi)\bar{f}(\xi)\operatorname{d\\!}\gamma\;\;\forall\chi\\!\in\\!\gamma$
(3)

where now, the integral on the left is also weakly singular, and
$(\bar{u}(\chi)-\bar{u}(\xi))$ may be thought of as “cancelling” one of the
$\mathinner{\\!\left\lVert r\right\rVert}$s in the denominator of $K$ as
$\chi\xrightarrow{}\xi$. This form allows using the method of [18], involving
rotation of the spherical coordinate system to align the north pole with the
source point, $\chi$, the eigenfunction property and the additional theorem of
spherical harmonics, to compute both integrals (see [19] for further details).
Furthermore, when the BIE is linearized in equation (13) below, the resulting
integrals remain weakly singular, and the method of [18] can be used to
compute them as well.

Introducing spherical coordinates, we have

$\begin{gathered}\hat{u}(\varphi_{\chi})+\frac{\lambda-1}{8\pi}\\!\int_{\phi=0}^{2\pi}\\!\int_{\theta=0}^{\pi}K\left(\hat{\xi}(\varphi_{\chi}),\hat{\xi}(\varphi)\right)(\hat{u}(\varphi)-\hat{u}(\varphi_{\chi}))J_{\hat{\xi}}(\varphi)\operatorname{d\\!}\theta\operatorname{d\\!}\phi\\\
=\hat{u}^{\infty}(\varphi_{\chi})-\frac{1}{8\pi\mu}\\!\int_{\phi=0}^{2\pi}\\!\int_{\theta=0}^{\pi}G\left(\hat{\xi}(\varphi_{\chi}),\hat{\xi}(\varphi)\right)\underaccent{\circ}{\hat{f}}(\varphi)\sin\theta\operatorname{d\\!}\theta\operatorname{d\\!}\phi\\\
\;\,\forall\varphi_{\chi}\\!\in\\![0,\pi]\times[0,2\pi]\end{gathered}$ (4)

where $\varphi_{\chi}$ denotes the spherical coordinates that map to the
source point, $\chi$. Note that the force on the membrane in equation (4) is
expressed per unit area of the unit sphere parent geometry,
$\underaccent{\circ}{\hat{f}}$, since this can be computed readily by equation
(10) below. Moreover, since the Jacobian $J_{\hat{\xi}}$ does not appear in
the right hand side integral of equation (4), the linearization in equation
(13) below has a simpler form.

##  3 Summary of elastic thin-shell mechanics results

This section summarizes the computation of the membrane force. The formulas
are completely equivalent to standard ones in the literature [15, 23, 35].
However, they are arranged so that the force can be computed using only the
operations of the gradient of scalar functions and the divergence of vector
fields on the unit sphere, without invoking Christoffel symbols and
derivatives of higher-order tensor fields. The gradient and divergence are
computed in spherical coordinates as

$\nabla\square=\left(\frac{\partial\square}{\partial\theta},\frac{1}{\sin\theta}\frac{\partial\square}{\partial\phi}\right)^{\top};\quad\nabla\cdot\square=\frac{1}{\sin\theta}\left(\frac{\partial}{\partial\theta}(\square_{\theta}\sin\theta)+\frac{\partial}{\partial\phi}\square_{\phi}\right)$
(5)

directly using the gradgc and divgc functions in SPHEREPACK. Note that the
Jacobians in Figure 1 and normal vector fields to the reference and deformed
surfaces can also be computed using the gradient operator,

$\begin{gathered}J_{\breve{\Xi}}=\mathinner{\\!\left\lVert\nabla_{\theta}\hat{\Xi}\times\nabla_{\phi}\hat{\Xi}\right\rVert};\quad\tilde{N}(\varphi)=\frac{1}{J_{\breve{\Xi}}}\left(\nabla_{\theta}\hat{\Xi}\times\nabla_{\phi}\hat{\Xi}\right)\\\
J_{\breve{\xi}}=\mathinner{\\!\left\lVert\nabla_{\theta}\hat{\xi}\times\nabla_{\phi}\hat{\xi}\right\rVert};\quad\tilde{n}(\varphi)=\frac{1}{J_{\breve{\xi}}}\left(\nabla_{\theta}\hat{\xi}\times\nabla_{\phi}\hat{\xi}\right)\end{gathered}$
(6)

The force per unit area of $\mathbb{S}^{2}$ is computed by the following steps
(details can be found in [8]).

  1. 1.

Kinematics: The strain, $\epsilon$, and curvature change, $\kappa$, are
computed as

$\epsilon=\frac{1}{2}A^{-1}(a-A);\quad\kappa=A^{-1}(b-B)$ (7)

where $A=\nabla\hat{\Xi}^{\top}\nabla\hat{\Xi}$ and
$a=\nabla\hat{\xi}^{\top}\nabla\hat{\xi}$ are matrices closely related to the
metric tensor in the reference and deformed configurations; similarly,
$B=\nabla\hat{\Xi}^{\top}\nabla\tilde{N}$ and
$b=\nabla\hat{\xi}^{\top}\nabla\tilde{n}$ are related to the respective
curvatures.

  2. 2.

Constitutive equations: The membrane in-plane force and moment matrices,
$\mathfrak{n}$ and $\mathfrak{m}$, are then obtained using constitutive
equations,

$\displaystyle\mathfrak{n}$
$\displaystyle=E_{\text{S}}I+\left(\frac{E_{\text{D}}}{2}\log(\det(C))-E_{\text{S}}\right)C^{-\top}$
(8) $\displaystyle\mathfrak{m}$ $\displaystyle=E_{\text{B}}\kappa^{\top}$

where $C=2\epsilon+I$, and following the notation of [25, 35], $E_{\text{D}}$
and $E_{\text{S}}$ are in-plane dilatational and shear moduli, and
$E_{\text{B}}$ is a bending modulus. A neo-Hookean model is used for the in-
plane forces, and a linear model for the moment.

  3. 3.

Equilibrium: The membrane out-of-plane shear, $\tilde{n}\mathfrak{q}$ is
computed as

$\tilde{n}\mathfrak{q}=\frac{1}{J_{\breve{\Xi}}J_{\breve{\xi}}}\left[\begin{array}[]{c|c}\nabla_{\phi}\hat{\xi}\times\mathcal{P}\mathfrak{r}&\mathcal{P}\mathfrak{r}\times\nabla_{\theta}\hat{\xi}\end{array}\right]$
(9)

where the normal projection, $\mathcal{P}=I-\tilde{n}\tilde{n}^{\top}$, and
$\mathfrak{r}=-\nabla\cdot\left(\nabla\hat{\xi}A^{-1}\mathfrak{m}J_{\breve{\Xi}}\right)$.
Finally, the membrane force is calculated using

$\underaccent{\circ}{\hat{f}}(\varphi)=-\nabla\cdot\left[\left(\nabla\hat{\xi}A^{-1}\mathfrak{n}+\nabla\tilde{n}A^{-1}\mathfrak{m}+\tilde{n}\mathfrak{q}\right)J_{\breve{\Xi}}\right]$
(10)

##  4 Galerkin approach and implicit ordinary differential equation (IODE)
structure

The residual of equation (4), which we call
$\text{res}(\hat{\xi},\hat{u};\chi)$, may be computed for any source point,
$\chi$. The BIE is satisfied if this residual is zero for all source points.
Instead of enforcing the BIE at selected source points, the residual is
computed using all the sample points as source points, and the spherical
harmonic transform of the residual is set to zero. This amounts to a Galerkin
approach with the test functions being the spherical harmonic basis functions
(see [19] for details). The process is depicted in Figure 2, and may be
thought of as the equation,

$g(\hat{\xi},\hat{u},t)=0$ (11)

Figure 2: Implicit ODE structure of BIE

Since the velocity field, $\hat{u}\\!=\\!\dot{\hat{\xi}}$, is the rate of
change of configuration, this is an implicit ordinary differential equation
(IODE).

Numerical solutions to such IODE could be computed by means of a number of
solvers such as ode15i in the MATLAB suite [31], DASSL [6], DASPK [7] and IDAS
in the SUNDIALS suite [21] that use multi-step methods based on backward
difference formulas (BDF), and RADAU5 [20] based on an implicit Runge-Kutta
method. All such solvers require the derivatives, or Jacobians,
$\operatorname{D\\!}_{\,1}\\!g$ and $\operatorname{D\\!}_{\,2}\\!g$, to use in
different variations of Newton’s method at each time-step. Computing these
Jacobian matrices is cumbersome and computationally expensive. However, as
outlined below in Section 5, the directional derivatives (or matrix-vector
products),
$\operatorname{D\\!}_{\,1}\\!g(\hat{\xi},\hat{u},t)\\!\cdot\\!\delta\hat{\xi}$
and
$\operatorname{D\\!}_{\,2}\\!g(\hat{\xi},\hat{u},t)\\!\cdot\\!\delta\hat{u}$
along $\delta\hat{\xi}$ and $\delta\hat{u}$ can be readily computed. The
solvers DASPK and IDAS allow a function that returns these matrix-vector
products. We use IDAS with its MATLAB interface. This solver requires two
functions representing the IODE,

  1. 1.

A function that given the state $(\hat{\xi},\hat{u})$ and time, computes the
residual $g(\hat{\xi},\hat{u},t)$; this is implemented as shown in Figure 2.

  2. 2.

A function that given the state, $(\hat{\xi},\hat{u})$, time, directions
$(\delta\hat{\xi},\delta\hat{u})$, and a parameter, $\Lambda$, computes the
sum of the directional derivatives,
$\operatorname{D\\!}_{\,1}\\!g(\hat{\xi},\hat{u},t)\\!\cdot\\!\delta\hat{\xi}+\Lambda\operatorname{D\\!}_{\,2}\\!g(\hat{\xi},\hat{u},t)\\!\cdot\\!\delta\hat{u}$.
This is discussed in the next section.

##  5 Directional derivatives

This section is on the computation of the directional derivatives of the IODE
residual needed for implicit numerical solution. The directional derivatives,
$\operatorname{D\\!}_{\,1}\\!\text{res}(\hat{\xi},\hat{u};\varphi_{\chi})\\!\cdot\\!\delta\hat{\xi}$
and
$\operatorname{D\\!}_{\,2}\\!\text{res}(\hat{\xi},\hat{u};\varphi_{\chi})\\!\cdot\\!\delta\hat{u}$
residual of equation (11) are calculated, and required derivative of the IODE
is obtained after a transformation as shown in Figure 3.

Figure 3: Structure of function to compute directional derivatives (Jacobian
matrix-vector products) of the BIE IODE

The second of these directional derivatives is simple, since the BIE (4) is
linear in the velocity field, and can be computed in the same manner as the
BIE itself.

$\begin{gathered}\operatorname{D\\!}_{\,2}\\!\text{res}(\hat{\xi},\hat{u};\varphi_{\chi})\\!\cdot\\!\delta\hat{u}=\delta\hat{u}(\varphi_{\chi})\\\
+\frac{\lambda-1}{8\pi}\\!\int_{\phi=0}^{2\pi}\\!\int_{\theta=0}^{\pi}K\left(\hat{\xi}(\varphi_{\chi}),\hat{\xi}(\varphi)\right)(\delta\hat{u}(\varphi)-\delta\hat{u}(\varphi_{\chi}))J_{\hat{\xi}}(\varphi)\operatorname{d\\!}\theta\operatorname{d\\!}\phi\end{gathered}$
(12)

The derivative in a direction of configuration change requires more careful
consideration.

$\displaystyle\operatorname{D\\!}_{\,1}\\!\text{res}(\hat{\xi},\hat{u};\varphi_{\chi})\\!\cdot\\!\delta\hat{\xi}$
$\displaystyle=\frac{\lambda-1}{8\pi}\\!\int_{\phi=0}^{2\pi}\\!\int_{\theta=0}^{\pi}\left[\delta\mathcal{K}\Delta\hat{u}J_{\breve{\xi}}+\mathcal{K}\Delta\hat{u}\delta
J_{\breve{\xi}}\right]\sin\theta\operatorname{d\\!}\theta\operatorname{d\\!}\phi$
(13)
$\displaystyle-\delta\hat{u}^{\infty}(\varphi_{\chi})+\frac{1}{8\pi\mu}\\!\int_{\phi=0}^{2\pi}\\!\int_{\theta=0}^{\pi}\left[\delta\mathcal{G}\underaccent{\circ}{\hat{f}}+\mathcal{G}\delta\underaccent{\circ}{\hat{f}}\right]\sin\theta\operatorname{d\\!}\theta\operatorname{d\\!}\phi$

where we have used the abbreviations,
$\Delta\hat{u}\\!=\\!\hat{u}(\varphi)-\hat{u}(\varphi_{\chi})$,
$\mathcal{K}\\!=\\!K(\hat{\xi}(\varphi_{\chi}),\hat{\xi}(\varphi))$ and
$\mathcal{G}\\!=\\!G(\hat{\xi}(\varphi_{\chi}),\hat{\xi}(\varphi))$, and
$\delta\square$ to represent the directional derivative,
$\operatorname{D\\!}\square(\hat{\xi})\\!\cdot\\!\delta\hat{\xi}$. In the
following, we define
$\nu=\nabla_{\theta}\hat{\xi}\times\nabla_{\phi}\hat{\xi}$. The directional
derivatives in equation (13), obtained by straightforward but tedious
calculation, are given by

$\begin{gathered}\delta
r=\delta\hat{\xi}(\varphi_{\chi})-\delta\hat{\xi}(\varphi);\quad\delta\tilde{r}=\frac{1}{\mathinner{\\!\left\lVert
r\right\rVert}}(I-\tilde{r}\tilde{r}^{\top})\delta r\\\
\delta\nu=\nabla_{\theta}\delta\hat{\xi}(\varphi)\times\nabla_{\phi}\hat{\xi}(\varphi)+\nabla_{\theta}\hat{\xi}(\varphi)\times\nabla_{\phi}\delta\hat{\xi}(\varphi)\\\
\delta
J_{\breve{\xi}}=\tilde{n}^{\top}\delta\nu;\;\;\delta\tilde{n}=\frac{1}{J_{\breve{\xi}}}\mathcal{P}\delta\nu;\;\;\delta\mathcal{G}\underaccent{\circ}{\hat{f}}=\frac{1}{\mathinner{\\!\left\lVert
r\right\rVert}}\left[-(\tilde{r}^{\top}\delta
r)\mathcal{G}\underaccent{\circ}{\hat{f}}+(\tilde{r}^{\top}\underaccent{\circ}{\hat{f}})\delta\tilde{r}+(\delta\tilde{r}^{\top}\underaccent{\circ}{\hat{f}})\tilde{r}\right]\\\
\begin{aligned}
\delta\mathcal{K}\Delta\hat{u}=\frac{6}{\mathinner{\\!\left\lVert
r\right\rVert}^{2}}&\left[(\tilde{r}^{\top}\tilde{n})(\tilde{r}^{\top}\Delta\hat{u})\left(\delta\tilde{r}-\frac{2}{\mathinner{\\!\left\lVert
r\right\rVert}}(\tilde{r}^{\top}\delta r)\tilde{r}\right)\right.\\\
&\left.+\left((\tilde{r}^{\top}\Delta\hat{u})(\delta\tilde{r}^{\top}\tilde{n}+\tilde{r}^{\top}\delta\tilde{n})+(\tilde{r}^{\top}\tilde{n})(\delta\tilde{r}^{\top}\Delta\hat{u})\right)\tilde{r}\vphantom{\left(\frac{2}{\mathinner{\\!\left\lVert
r\right\rVert}}\right)}\right]\end{aligned}\end{gathered}$ (14)

The directional derivative of the membrane force is obtained by linearizing
the shell mechanics relations from Section 3.

$\displaystyle\delta\epsilon$
$\displaystyle=\frac{1}{2}A^{-1}\left(\nabla\delta\hat{\xi}^{\top}\nabla\hat{\xi}+\nabla\hat{\xi}^{\top}\nabla\delta\hat{\xi}\right)$
(15) $\displaystyle\delta C$ $\displaystyle=2\delta\epsilon$
$\displaystyle\delta\kappa$
$\displaystyle=A^{-1}\left(\nabla\delta\hat{\xi}^{\top}\nabla\tilde{n}+\nabla\hat{\xi}^{\top}\nabla\delta\tilde{n}\right)$
$\displaystyle\delta\mathfrak{n}$
$\displaystyle=\left[\left(\frac{E_{\text{D}}}{2}\operatorname{tr}(C^{-1}\delta
C)\right)-\left(\frac{E_{\text{D}}}{2}\log(\det(C))-E_{\text{S}}\right)\delta
C\right]C^{-\top}$ $\displaystyle\delta\mathfrak{m}$
$\displaystyle=E_{\text{B}}\delta\kappa^{\top}$
$\displaystyle\delta\mathfrak{r}$
$\displaystyle=-\nabla\cdot\left[\left(\nabla\delta\hat{\xi}A^{-1}\mathfrak{m}+\nabla\hat{\xi}A^{-1}\delta\mathfrak{m}\right)J_{\breve{\Xi}}\right]$
$\displaystyle\delta(\mathcal{P}\mathfrak{r})$
$\displaystyle=\mathcal{P}\delta\mathfrak{r}-(\tilde{n}^{\top}\mathfrak{r})\delta\tilde{n}-(\delta\tilde{n}^{\top}\mathfrak{r})\tilde{n}$
$\displaystyle\delta(\tilde{n}\mathfrak{q})$ $\displaystyle=-\frac{\delta
J_{\breve{\xi}}}{J_{\breve{\xi}}}(\tilde{n}\mathfrak{q})+\frac{1}{J_{\breve{\Xi}}J_{\breve{\xi}}}\left[\begin{array}[]{c|c}\\!\nabla_{\phi}\delta\hat{\xi}\\!\times\\!\mathcal{P}\mathfrak{r}\\!+\\!\nabla_{\phi}\hat{\xi}\\!\times\\!\delta(\mathcal{P}\mathfrak{r})&\\!\mathcal{P}\mathfrak{r}\\!\times\\!\nabla_{\theta}\delta\hat{\xi}\\!+\\!\delta(\mathcal{P}\mathfrak{r})\\!\times\\!\nabla_{\theta}\hat{\xi}\end{array}\right]$
$\displaystyle\delta\underaccent{\circ}{\hat{f}}(\varphi)$
$\displaystyle=-\nabla\\!\cdot\\!\left[\left(\nabla\delta\hat{\xi}A^{-1}\mathfrak{n}+\nabla\hat{\xi}A^{-1}\delta\mathfrak{n}+\nabla\delta\tilde{n}A^{-1}\mathfrak{m}+\nabla\tilde{n}A^{-1}\delta\mathfrak{m}+\delta(\tilde{n}\mathfrak{q})\right)J_{\breve{\Xi}}\right]$

Equation (13) is itself a BIE with only weakly singular terms, and the
integrals can again be computed using the method of [18]. The weakly singular
representation of BIE (4) facilitates the linearization in this convenient
form.

##  6 Numerical examples

Four numerical examples are presented to demonstrate the performance of the
proposed fully implicit method. Various aspects of performance such as (i)
time steps used to achieve the specified tolerances, (ii) order of backward
difference approximation used, (iii) number of function and derivative calls
needed, (iv) conservation of different quantities, and (v) comparison with the
explicit time scheme, as well as the dependence of these on spatial
discretization are illustrated. In each example, the motion of a red blood
cell suspended in an ambient flow is considered. The four examples are

  1. (a)

Shear flow ($k\\!=\\!100$\text{\,}{\mathrm{s}}^{-1}$$): An ambient shear flow
parallel to the $x$ direction with a shear strain rate,
$k=100$\text{\,}{\mathrm{s}}^{-1}$$ (as in [24, 33]), so that ambient velocity
field, $u=(kz,0,0)^{\top}$, with an undeformed red cell placed initially
making an angle of $45^{\circ}$ with the $x$ axis. Capillary number,
$\text{Ca}=\frac{\mu ka}{E_{s}}=0.08$.

  2. (b)

Shear flow ($k\\!=\\!500$\text{\,}{\mathrm{s}}^{-1}$$): Same as (a), but with
a greater shear rate, $k\\!=\\!500$\text{\,}{\mathrm{s}}^{-1}$$ (as in [24,
33]), $\text{Ca}=0.4$.

  3. (c)

Parabolic flow: A parabolic ambient flow along the $x$ direction,
$u_{x}=A(B-(y^{2}+z^{2}))$, where A = 1
($35.46\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{-1}\text{\,}{\mathrm{s}}^{-1}$) and B = 2.3
($18.3\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$) are parameters that
control the curvature and peak velocity of the flow, $\text{Ca}=\frac{\mu
Aa^{4}}{E_{B}}=14.95$.

  4. (d)

Relaxation: A highly deformed cell is suspended in ambient fluid at rest until
it relaxes to its undeformed shape. A deformed shape from (b) is taken as the
initial configuration.

Following common practice in the literature [24, 25, 35], non-dimensional
parameters are used. The equivalent cell radius,
$a\approx$2.82\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ is used as
characteristic length, the inverse of a shear rate $k=100$
$\text{\,}{\mathrm{s}}^{-1}$ as characteristic time, and plasma viscosity
times this shear rate $\mu k=$0.12\text{\,}\mathrm{Pa}$$ as characteristic
stress, where
$\mu=$1.2\text{\times}{10}^{-3}\text{\,}\mathrm{Pa}\cdot\mathrm{s}$$. With
this nondimensionalization, the cell membrane shear modulus,
$E_{\text{S}}=12.4$ $(4.2\times 10^{-6}\text{N}/\text{m})$ and bending
modulus, $E_{\text{B}}=0.0669$ $(1.8\times 10^{-19}\text{N}\cdot\text{m}$).
The membrane dilatational modulus is set as a penalty parameter
$E_{\text{D}}=200$ $(6.8\times 10^{-5}\text{N}/\text{m})$ to control membrane
area change while preserving numerical conditioning. Viscosity ratio of the
intracellular plasma to extracellular plasma, $\lambda=5$. The ability of the
membrane to deform under viscous stress is gauged by capillary number,
$\text{Ca}=\frac{\mu ka}{E_{s}}$, based on shear modulus for case (a) and (b),
and $\text{Ca}=\frac{\mu Aa^{4}}{E_{B}}$, based on bending modulus for case
(c) [16]. The initial undeformed shape of the cell is a biconcave disk, which
following [24, 25] is described by

$\displaystyle\hat{\Xi}(\varphi)=\left(\alpha\sin{\theta}\cos{\phi};\alpha\sin{\theta}\sin{\phi};\frac{\alpha}{2}(0.207+2.003\sin^{2}\theta-1.123\sin^{4}\theta)\cos{\theta}\right)^{\top}$
(16)

where $\alpha\\!=\\!1.386$ is the non-dimensional maximum radius of the
biconcave disk.

The relative tolerance for truncation error is set to $10^{-4}$ and absolute
tolerance to $10^{-6}$ for the IDAS solver for all cases. A maximum time step
restriction is set for each example to minimize failed steps.

Figure 4 shows deformation snapshots of the evolving cell for the different
flow types with spherical harmonic degree $N=24$. The color indicates the norm
of the principal membrane stress. Figure 5 presents the profile snapshots of
the cell cross section in the $xz$ plane corresponding to the deformed shapes
in Figure 4. The behaviors seen in these figures are well documented in the
literature [24, 11]. From Figures 4(a, b) and 5(a, b), it can be seen that the
cell rotates clockwise, and is stretched or compressed depending on the angle
of inclination. Compared to case (a), the cell in case (b) has larger
deformation due to higher shear rate, and the period of rotation is much
higher as well. A star-shape marker serves as tracker of the cell rotating.
The reversed S-shape of the cell mid-plane profile shown in Figure 5(b) is
also reported by Pozrikidis [24]. In Figures 4(c) and 5(c), the cell is
stretched at the middle as it translates in the $x$ direction under parabolic
flow. Due to the initial inclination, the upper half of the cell moves faster
and drags the lower half forward. In the relaxation case, the cell recovers
from the highly deformed shape to reach the stress-free state in ambient
fluid, as shown in Figure 4(d) and 5(d).

|  Flow type |  Deformation snapshots
---|---|---
|  |
(a) |  Shear flow k = 100 $\text{\,}{\mathrm{s}}^{-1}$ |
|  |
(b) |  Shear flow k = 500 $\text{\,}{\mathrm{s}}^{-1}$ |
|  |
(c) |  Parabolic flow |
|  |
(d) |  Relaxation |
Figure 4: Evolution of cell deformation ($N$=24). Colors indicate the norm of the principal membrane stress. |  Flow type |  Profile snapshots
---|---|---
|  |
(a) |  Shear flow k = 100 $\text{\,}{\mathrm{s}}^{-1}$ |
|  |
(b) |  Shear flow k = 500 $\text{\,}{\mathrm{s}}^{-1}$ |
|  |
(c) |  Parabolic flow |
|  |
(d) |  Relaxation |
Figure 5: Cell mid-plane profiles corresponding to Figure 4

In the following, various numerical performance aspects of the implicit time
integration scheme are discussed:

  1. (i)

Time step size: The time-step sizes used by the implicit scheme and the
explicit scheme for different spatial resolutions and flow types are reported
in Table 1. Time-step size for implicit method is nearly independent of the
spherical harmonic degree, and only depends on the flow type and capillary
number. A similar observation was made in [33] for simulation of a vesicle in
simple shear flow with a semi-implicit scheme.

| Explicit | Implicit
---|---|---
Flow type | (a) | (b) | (c) | (a) | (b) | (c) | (d)
$N=8$ | 1.13E-02 | 1.16E-02 | 1.13E-02 | 7.34E-02 | 1.60E-02 | 1.22E-01 | 5.14E-02
16 | 5.07E-03 | 5.63E-03 | 5.38E-03 | 5.10E-02 | 2.26E-02 | 1.17E-01 | 6.55E-02
24 | 3.07E-03 | 3.07E-03 | 3.46E-03 | 6.60E-02 | 2.45E-02 | 8.93E-02 | 5.90E-02

Table 1: Steady state time-step sizes used by the explicit and implicit time
integration schemes for different flow types and spherical harmonic degrees

(a)

(b)

Figure 6: IDAS solver statistics for example (a) for $N$ = 24 (solid line),
$N$ = 8 (dotted line) and with explicit method for $N$ = 8 (dash-dot line);
The first subplot shows evolution of time-step size with simulation time —
after adapting the time step, the solver settles at a time step of $6.60\times
10^{-2}$ for $N$ = 24, the second subplot shows the the order of backward
difference approximation, the third and fourth subplots show numbers of
residual and Jacobian-vector product evaluations for implicit method and
matrix-vector product for explicit method. The first and fourth subplots also
show corresponding data for explicit integration.

(a)

(b)

Figure 7: IDAS solver statistics for example (c) for $N$ = 24 (solid line),
$N$ = 8 (dotted line) and with explicit method for $N$ = 8 (dash-dot line);
The first subplot shows evolution of time-step size with simulation time —
after adapting the time step, the solver settles at a time step of $8.93\times
10^{-2}$ for $N$ = 24, the second subplot shows the the order of backward
difference approximation, the third and fourth subplots show numbers of
residual and Jacobian-vector product evaluations for implicit method and
matrix-vector product for explicit method. The first and fourth subplots also
show corresponding data for explicit integration.

IDAS solver statistics are provided for the shear flow case (a) in Figure 6
and for the parabolic flow case (c) in Figure 7. In each figure, solid lines
show data for $N\\!=\\!24$, and dotted line for $N\\!=\\!8$. The data with
explicit method for $N\\!=\\!8$ is also reported as dash-dot line in both
first and last subplots for comparison. The first subplot shows evolution of
time step size during the simulation. The solver adapts the time step based on
local truncation error criteria and on convergence of Newton iterations. For
each example, we set a maximum step size to minimize the number of failed
steps. The solver settles at a time step of $6.60\\!\times\\!10^{-2}$ in case
(a), and $8.93\\!\times\\!10^{-2}$ in case (c), more than ten times what is
possible with an explicit method based on stability considerations alone.

  2. (ii)

Order of backward difference approximation: The second subplot in each of
Figures 6 and Figures 7 shows how the order of backward difference
approximation employed by the solver. The order is adapted between 1 and 5,
again based on local truncation error test and convergence. The order of
approximation is at least 3, and most often is the highest possible order 5,
indicating the smoothness of the solution.

  3. (iii)

Number of function and derivative calls needed: The third and fourth subplots
in Figures 6 and Figures 7 report the numbers of residual function evaluations
(Figure 2) and Jacobian matrix-vector multiplications (Figure 3).
Approximately ten Jacobian multiplications are needed for every residual
evaluation. Since fewer time steps are used, the total number of matrix-vector
products is comparable to the number of GMRES iterations needed in an explicit
method. The number of matrix-vector products can potentially be reduced by
suitably preconditioning the system.

  4. (iv)

Conservation of different quantities: Besides local truncation error, as other
measures of accuracy, we consider conservation of membrane surface area and
cell volume, and equilibrium of membrane forces. Figure 8 shows membrane
surface area change and volume conservation error for the shear case (a).
Volume conservation error decreases dramatically with increasing spherical
harmonic degree due to the accuracy of representation of the BIE itself. The
area changes is within $1\%$, indicating the extent to which the penalty
parameter, $E_{\text{D}}$, is effective in enforcing membrane inextensibility.
The largest area changes occur at the times when cell is stretched most by the
shear flow when the angle of inclination of the cell is $\pm 45^{\circ}$.
Figure 9 demonstrates equilibrium of forces and moments on the entire cell.

  5. (v)

Comparison with the explicit integration Due to the variabilities in how time
step size can be adapted, tolerances for nonlinear and linear solvers,
truncation error control etc., it is not possible to make an exact
quantitative comparison between the performances of implicit and explicit time
integration schemes. Qualitatively, however, it is clear that steady state
time-step size used by the implicit time integration scheme depends on the
imposed flow type and is not affected by the spherical harmonics order (see
Table 1). On the other hand, the steady time-step size used by the explicit
method decreases with increasing spherical harmonic degree and, interestingly,
is not affected by the flow strength and type. The time step dictated by
stability in an explicit method is approximately 1/10th of the time step size
used by IDAS for truncation error control. Comparison of computational effort
between the IDAS and an explicit method is made using the number of matrix-
vector multiplications (although a Jacobian matrix-vector multiplication is
more expensive than the residual matrix-vector multiplication). The number of
matrix-vector multiplication performed by IDAS is approximately 1/7th of the
number of GMRES iterations in an explicit method as seen in the fourth
subplots of Figures 6 and 7.

Figure 8: Membrane surface area change and volume conservation error for shear
flow example (a) with $k=100$ $\text{\,}{\mathrm{s}}^{-1}$ for different
spherical harmonic degrees $N$. Figure 9: Summation of membrane forces and
moments for shear flow example (a) with $k=100$ $\text{\,}{\mathrm{s}}^{-1}$
for $N$ = 24 (solid line) and $N$ = 8 (dotted line).

##  7 Concluding remarks

An implicit time integration approach has been presented in connection with a
spectral boundary integral method for computation of red blood cell flow. The
simple problem of cell flow in an ambient fluid (rather than in a vessel) has
been considered to describe of the method. The flow is represented by a
boundary integral equation (BIE) expressed in weakly singular form. Cell
configurations are represented by smooth deformations of the unit sphere. The
cell configuration and velocity field are discretized using spherical harmonic
basis functions. The cell membrane is modeled using elastic thin-shell
mechanics, formulas from which are written in a manner that requires only the
operations of gradient of scalar functions and divergence of vector fields on
the unit sphere in spherical coordinates. This BIE representing the flow has
the form of an implicit ordinary differential equation (IODE) in the cell
configuration, and solutions are computed using an off-the-shelf IODE solver
based on backward differences formulas. The solver adapts the time step and
order of backward difference approximation to control local truncation error
and to regulate convergence of Newton’s method. Newton’s method employed in
implicit time integration necessitates computation of Jacobians of the IODE
with respect to the configuration and the velocity field. Rather than assemble
such Jacobians, Jacobian-vector products, i.e. directional derivatives are
calculated. This calculation, which itself is the evaluation of a different
BIE, is aided by the expression of the flow BIE in weakly singular format.
Numerical examples have been presented to demonstrate the working of the
implicit approach. They illustrate that larger time steps can be used compared
to explicit methods, where time-step size is limited by stability
considerations. The number of Jacobian-vector products is comparable to the
number of GMRES iterations of the flow BIE in explicit methods. Source code
for the implementation of the implicit method is included online supplemental
material that accompany this note.

Much further development is needed to employ this implicit method in broader
problems of interest. Suitable preconditions must be constructed for the
directional derivative BIEs to reduce the number of Jacobian-vector products
further, and hence computational effort. Fast summation methods must also be
implemented for the derivative BIEs. Most importantly, the method must be
extended for the presence of other geometries such as vessels and slits. While
these do not alter the mathematical formulation significantly, they raise the
possibility of interpenetration when large time steps are used, In such cases,
interpenetration must be flagged to the solver. Closest-point calculations
used in evaluation of near-singular integrals can potentially be adapted to
detect interpenetration. These are topics of current inquiry.

## References

  * [1] John C Adams and Paul N Swarztrauber.  SPHEREPACK 3.0: A model development facility.  Monthly Weather Review, 127(8):1872–1878, 1999.
  * [2] Uri M Ascher and Linda R Petzold.  Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations.  SIAM: Society for Industrial and Applied Mathematics, 1998.
  * [3] C. Bae, F. Sachs, and P. A. Gottlieb.  The mechanosensitive ion channel PIEZO1 is inhibited by the peptide GsMTx4.  Biochemistry, 50(29):6295–300, 2011.
  * [4] Chilman Bae, Radhakrishnan Gnanasambandam, Chris Nicolai, Frederick Sachs, and Philip A. Gottlieb.  Xerocytosis is caused by mutations that alter the kinetics of the mechanosensitive channel PIEZO1.  Proceedings of the National Academy of Sciences, 110(12):E1162–E1168, 2013.
  * [5] G. Boedec, M. Leonetti, and M. Jaeger.  3d vesicle dynamics simulations with a linearly triangulated surface.  Journal of Computational Physics, 230(4):1020–1034, 2011.
  * [6] Kathryn E. Brenan, Stephen L. Campbell, and Linda R. Petzold.  Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations.  SIAM, Philadelphia, 1996.
  * [7] Peter N. Brown, Alan C. Hindmarsh, and Linda R. Petzold.  Using Krylov methods in the solution of large-scale differential-algebraic systems.  SIAM Journal of Scientific Computing, 15(6):1467–1488, 1994\.
  * [8] Pei-Chuan Chao.  Spectral and implicit methods for computational cell mechanics.  PhD thesis, University at Buffalo, 2020.
  * [9] Ivan Cimràk and Ivetà Jančigova.  Computational Blood Cell Mechanics: Road Towards Models and Biomedical Applications.  CRC Press, 2018.
  * [10] E. Cinar, S. Zhou, J. DeCourcey, Y. Wang, R. E. Waugh, and J. Wan.  Piezo1 regulates mechanotransductive release of ATP from human RBCs.  Proceedins of the National Academy of Sciences, 112(38):11783–11788, 2015.
  * [11] Daniel Cordasco and Prosenjit Bagchi.  On the shape memory of red blood cells.  Physics of Fluids, 29(4):041901, 2017.
  * [12] Jens G. Danielczok, Emmanuel Terriac, Laura Hertz, Polina Petkova-Kirova, Franziska Lautenschläger, Matthias W. Laschke, and Lars Kaestner.  Red blood cell passage of small capillaries is associated with transient Ca${}^{\text{2+}}$-mediated adaptations.  Frontiers in Physiology, 8:979, 2017.
  * [13] Panagiotis Dimitrakopoulos.  Interfacial dynamics in stokes flow via a three-dimensional fully-implicit interfacial spectral boundary element algorithm.  Journal of Computational Physics, 225(1):408–426, 2007.
  * [14] Sai K. Doddi and Prosenjit Bagchi.  Three-dimensional computational modeling of multiple deformable cells flowing in microvessels.  Physical Review E, 79:046318, Apr 2009.
  * [15] Wilhelm Flügge.  Tensor analysis and continuum mechanics.  Springer-Verlag, Berlin, New York, 1972.
  * [16] Jonathan B. Freund.  The flow of red blood cells through a narrow spleen-like slit.  Physics of Fluids, 25(11):110807, 2013.
  * [17] Jonathan B Freund.  Numerical simulation of flowing blood cells.  Annual Review of Fluid Mechanics, 46:67–95, 2014.
  * [18] M. Ganesh and I. G. Graham.  A high-order algorithm for obstacle scattering in three dimensions.  Journal of Computational Physics, 198(1):211–242, 2004.
  * [19] Ali Gürbüz.  Simulation of Pressure-Driven Red Blood Cell Flow through Constrictions by Boundary Integral Methods.  PhD thesis, University at Buffalo, 2020.
  * [20] E. Hairer and G. Wanner.  Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems.  Springer, New York, 1991.
  * [21] Alan C Hindmarsh, Peter N Brown, Keith E Grant, Steven L Lee, Radu Serban, Dan E Shumaker, and Carol S Woodward.  SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers.  ACM Transactions on Mathematical Software (TOMS), 31(3):363–396, 2005.
  * [22] Yijun Liu and TJ Rudolphi.  Some identities for fundamental solutions and their applications to weakly-singular boundary element formulations.  Engineering Analysis with Boundary Elements, 8(6):301–311, 1991\.
  * [23] Frithiof I. Niordson.  Shell theory.  North-Holland series in applied mathematics and mechanics. North-Holland, Amsterdam, New York, 1985.
  * [24] C Pozrikidis.  Numerical simulation of the flow-induced deformation of red blood cells.  Annals of Biomedical Engineering, 31(10):1194–1205, 2003.
  * [25] C Pozrikidis.  Axisymmetric motion of a file of red blood cells through capillaries.  Physics of Fluids, 17(3):031503, 2005.
  * [26] Constantine Pozrikidis.  Boundary Integral and Singularity Methods for Linearized Viscous Flow.  Cambridge University Press, 1992.
  * [27] Bryan Quaife and George Biros.  Adaptive time stepping for vesicle suspensions.  Journal of Computational Physics, 306:478–499, 2016.
  * [28] Abtin Rahimian, Shravan K Veerapaneni, Denis Zorin, and George Biros.  Boundary integral method for the flow of vesicles with viscosity contrast in three dimensions.  Journal of Computational Physics, 298:766–786, 2015.
  * [29] S Ramanujan and C Pozrikidis.  Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: large deformations and the effect of fluid viscosities.  Journal of Fluid Mechanics, 361:117–143, 1998.
  * [30] Timothy W. Secomb.  Blood flow in the microcirculation.  Annual Review of Fluid Mechanics, 49(1):443–461, 2017.
  * [31] The Mathworks.  MATLAB version 9.6.0.1072779 (R2019a).  The Mathworks, Inc., Natick, Massachusetts, 2019.
  * [32] Shravan K Veerapaneni, Denis Gueyffier, George Biros, and Denis Zorin.  A numerical method for simulating the dynamics of 3D axisymmetric vesicles suspended in viscous flows.  Journal of Computational Physics, 228(19):7233–7249, 2009.
  * [33] Shravan K Veerapaneni, Abtin Rahimian, George Biros, and Denis Zorin.  A fast algorithm for simulating vesicle flows in three dimensions.  Journal of Computational Physics, 230(14):5610–5634, 2011.
  * [34] Ryan Zarychanski, Vincent P. Schulz, Brett L. Houston, Yelena Maksimova, Donald S. Houston, Brian Smith, Jesse Rinehart, and Patrick G. Gallagher.  Mutations in the mechanotransduction protein PIEZO1 are associated with hereditary xerocytosis.  Blood, 120(9):1908–1915, 2012.
  * [35] Hong Zhao, Amir HG Isfahani, Luke N Olson, and Jonathan B Freund.  A spectral boundary integral method for flowing blood cells.  Journal of Computational Physics, 229(10):3726–3744, 2010.
  * [36] Hong Zhao and Eric SG Shaqfeh.  The dynamics of a vesicle in simple shear flow.  Journal of Fluid Mechanics, 674:578–604, 2011.