# Probing Shocked Ejecta in SN 1987A with XRISM$-$Resolve: the effects of the
gate valve closed

Vincenzo Sapienza Dipartimento di Fisica e Chimica E. Segrè, Università degli
Studi di Palermo, Piazza del Parlamento 1, 90134, Palermo, Italy  INAF-
Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, 90134, Palermo,
Italy  Department of Physics, Graduate School of Science, The University of
Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  Marco Miceli
Dipartimento di Fisica e Chimica E. Segrè, Università degli Studi di Palermo,
Piazza del Parlamento 1, 90134, Palermo, Italy  INAF-Osservatorio Astronomico
di Palermo, Piazza del Parlamento 1, 90134, Palermo, Italy  Aya Bamba
Department of Physics, Graduate School of Science, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  Research Center for the Early
Universe, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,
Tokyo 113-0033, Japan  Trans-Scale Quantum Science Institute, The University
of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  Salvatore Orlando INAF-
Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, 90134, Palermo,
Italy  Shiu-Hang Lee Department of Astronomy, Kyoto University Oiwake-cho,
Kitashirakawa, Sakyo-ku, Kyoto 606-8502, Japan  Kavli Institute for the
Physics and Mathematics of the Universe (WPI), The University of Tokyo,
Kashiwa 277-8583, Japan  Shigehiro Nagataki RIKEN Interdisciplinary
Theoretical and Mathematical Sciences Program (iTHEMS), 2-1 Hirosawa, Wako,
Saitama 351-0198, Japan  Astrophysical Big Bang Laboratory (ABBL), RIKEN
Cluster for Pioneering Research, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan
Astrophysical Big Bang Group (ABBG), Okinawa Institute of Science and
Technology Graduate University, 1919-1 Tancha, Onna-son, Kunigami-gun, Okinawa
904-0495, Japan  Masaomi Ono RIKEN Interdisciplinary Theoretical and
Mathematical Sciences Program (iTHEMS), 2-1 Hirosawa, Wako, Saitama 351-0198,
Japan  Astrophysical Big Bang Laboratory (ABBL), RIKEN Cluster for Pioneering
Research, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan  Institute of Astronomy
and Astrophysics, Academia Sinica, Taipei 10617, Taiwan  Satoru Katsuda
Graduate School of Science and Engineering, Saitama University, 255 Simo-
Ohkubo, Sakura-ku, Saitama city, Saitama, 338-8570, Japan  Koji Mori  Faculty
of Engineering, University of Miyazaki, 1-1 Gakuen Kibanadai Nishi, Miyazaki,
Miyazaki 889-2192, Japan  Japan Aerospace Exploration Agency, Institute of
Space and Astronautical Science, 3-1-1 Yoshino-dai, Chuo-ku, Sagamihara,
Kanagawa 252-5210, Japan  Makoto Sawada Department of Physics, Rikkyo
University, 3-34-1 Nishi Ikebukuro, Toshima-ku, Tokyo 171-8501, Japan
Yukikatsu Terada Graduate School of Science and Engineering, Saitama
University, 255 Simo-Ohkubo, Sakura-ku, Saitama city, Saitama, 338-8570, Japan
Roberta Giuffrida Dipartimento di Fisica e Chimica E. Segrè, Università degli
Studi di Palermo, Piazza del Parlamento 1, 90134, Palermo, Italy  INAF-
Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, 90134, Palermo,
Italy  Fabrizio Bocchino INAF-Osservatorio Astronomico di Palermo, Piazza del
Parlamento 1, 90134, Palermo, Italy

###### Abstract

Supernova (SN) 1987A is widely regarded as an excellent candidate for
leveraging the capabilities of the freshly launched XRISM satellite. Recent
researches indicate that the X-ray emission from SN 1987A will increasingly
originate from its ejecta in the years to come. In a previous study, we
thoroughly examined the proficiency of XRISM-Resolve in identifying signatures
of shocked ejecta in SN 1987A, synthesizing the XRISM-Resolve spectrum based
on a state-of-the-art magneto-hydrodynamic simulation. However, following the
satellite’s launch, a technical issue arose with the XRISM instrument’s gate
valve, which failed to open, thereby affecting observations with the Resolve
spectrometer. Here, we update our analysis, reevaluating our diagnostic
approach under the assumption that the gate valve remains closed. We find
that, even with the reduced instrumental capabilities, it will be possible to
pinpoint the ejecta contribution through the study of the line profiles in the
XRISM-Resolve spectrum of SN 1987A.

##  1 Introduction

Supernova (SN) 1987A provides a unique opportunity to study the evolution of a
SN into its remnant across the electromagnetic spectrum (McCray & Fransson,
2016). The supernova remnant (SNR) is evolving in a clumpy equatorial ring
within a diffuse $\mathrm{H}\,$II region (Sugerman et al., 2005).

X-ray observations are crucial for studying the interaction between the shock
front, the circumstellar medium (CSM), and the ejecta. Since the initial
detection (Dotani et al. 1987; but see also Beuermann et al. 1994), the X-ray
emission increased for about 25 years, confirming the shock’s encounter with
the equatorial ring (Borkowski et al. 1997,Park et al. 2005, Haberl et al.
2006, Zhekov et al. 2009, Maggi et al. 2012). Recent investigations in the
past few years, hinted that the X-ray emission from the CSM is decreasing,
leaving space to the X-ray emission from the ejecta (Frank et al. 2016, Sun et
al. 2021, Ravi et al. 2021, Maitra et al. 2022).

3D hydrodynamic and magneto-hydrodynamic simulations have replicated the X-ray
light curves, spectra and remnant morphology, as well as other observed
properties of SN 1987A. Initially dominated by the shocked $\mathrm{H}\,$II
region, the X-ray emission later became influenced by the shocked equatorial
ring, and is now driven by the outer ejecta heated by the reverse shock
(Orlando et al., 2015, 2019, 2020).

A novel data analysis approach by Miceli et al. (2019) using Chandra data
compared with 3D simulations provided insights into post-shock temperatures.
The XRISM mission, launched on September 7, 2023, will offer high-resolution
X-ray spectroscopy, allowing in-depth studies of SN 1987A. Using our in-house
tool, we synthesized the XRISM-Resolve X-ray spectrum of SN 1987A from MHD
simulations (Orlando et al., 2020; Ono et al., 2020), successfully tracing new
ejecta diagnostics for future observations (Sapienza et al., 2024).

However, the XRISM Telescope safety gate valve for the Resolve instrument
failed to open after launch. As a consequence, the effective area of XRISM-
Resolve results degraded and severely compromised under 2 keV. Here we update
the analysis in Sapienza et al. (2024) by adopting the updated response files
(Section 2) and investigate the possibility of detecting the ejecta in SN
1987A with the gate valve closed (Section 3).

##  2 Results

Following the same procedure described in Sapienza et al. (2024), we produced
the synthetic spectrum of SN 1987A for the year 2024, using the new RMF and
ARF files that take into account the effects of the gate valve closed,
available in the HEASARC XRISM
website111https://heasarc.gsfc.nasa.gov/docs/xrism/proposals/responses.html.

Figure 1: Upper panel:Synthetic XRISM-Resolve spectrum of SN 1987A for the
year 2024 (assuming the gate valve closed) with an exposure time of 160 ks
including bulk motion Doppler broadening and binned using the Kaastra &
Bleeker (2016) optimal binning algorithm. The synthetic XRISM-Resolve spectral
model is superimposed with a black dashed line. The red curve shows the
contribution of the CSM (ER and $\mathrm{H}\,$II) to the spectral model, while
the blue curve shows the contribution of the ejecta. Lower left panel: Close-
up view of the spectrum shown in upper panel in the $2.38-2.65$ keV band with
best-fit model and residuals. Lower right panel: Close-up view of the spectral
model shown in upper panel in the $2.38-2.65$ keV band.

Upper panel of Figure 1 shows the updated version of the synthetic spectrum,
for an exposure time of 160 ks, which is the amended exposure time for the PV
phase observation of SN 1987A, increased with respect to the original 100 ks
to compensate the loss of effective area. Unfortunately, the effective area
below 2 keV is severely compromised, hampering us to carry on the diagnostic
procedure for the Mg and Si emission lines. However, the quality of the
spectrum remains adequate for applying the same diagnostics to the S lines. In
particular, it will be possible to reveal the ejecta through their enhanced
line widths (with respect to the narrow lines stemming from the shocked
ambient medium) associated with their large velocities. We then fitted the S
line profiles adopting only 2 broad Gaussian components (due to the reduced
statistic), for the S XV w-z and for the S XVI Ly$\alpha$, modeling the
continuum emission with a power law and accounting for the interstellar and
intergalactic absorption. Similarly to Sapienza et al. (2024)
$\sigma_{2}=\sigma_{1}\cdot E_{2}/E_{1}$, being $\sigma_{1,2}$ and $E_{1,2}$
the width and the centroid of the ejecta components in the lines of He-like
and H-like species. Figure 1 (lower-left panel) shows a close-up view of the
synthetic spectrum highlighting the S lines with the corresponding best fit
model and residuals. The different components are detected with high
significance (their normalization being larger than zero at more than the
5$\sigma$ confidence level) and clearly show a large expansion velocity
($3300_{-700}^{+800}$ km s-1), which nicely accounts for the ejecta bulk
motion.

##  3 Discussion and Conclusions

We synthesized the updated PV XRISM-Resolve X-ray spectrum for 2024, including
the effects of the gate valve closed. Our analysis shows that the synthetic
spectra still provides valuable diagnostics to measure the ejecta velocity
from the broadening of sulfur emission lines, highlighting the continued
utility of XRISM in studying SN 1987A.

Our results showcase the adaptability and resilience of our analytical methods
and the importance of ongoing X-ray observations and simulations in unveiling
the complex dynamics of supernova remnants. Future observations, leveraging
both current and upcoming missions, are expected to provide deeper insights
into the evolving nature of SN 1987A.

## References

  * Beuermann et al. (1994) Beuermann, K., Brandt, S., & Pietsch, W. 1994, A&A, 281, L45
  * Borkowski et al. (1997) Borkowski, K. J., Blondin, J. M., & McCray, R. 1997, ApJ, 477, 281, doi: 10.1086/303691
  * Dotani et al. (1987) Dotani, T., Hayashida, K., Inoue, H., et al. 1987, Nature, 330, 230, doi: 10.1038/330230a0
  * Frank et al. (2016) Frank, K. A., Zhekov, S. A., Park, S., et al. 2016, ApJ, 829, 40, doi: 10.3847/0004-637X/829/1/40
  * Haberl et al. (2006) Haberl, F., Geppert, U., Aschenbach, B., & Hasinger, G. 2006, A&A, 460, 811, doi: 10.1051/0004-6361:20066198
  * Kaastra & Bleeker (2016) Kaastra, J. S., & Bleeker, J. A. M. 2016, A&A, 587, A151, doi: 10.1051/0004-6361/201527395
  * Maggi et al. (2012) Maggi, P., Haberl, F., Sturm, R., & Dewey, D. 2012, A&A, 548, L3, doi: 10.1051/0004-6361/201220595
  * Maitra et al. (2022) Maitra, C., Haberl, F., Sasaki, M., et al. 2022, A&A, 661, A30, doi: 10.1051/0004-6361/202141104
  * McCray & Fransson (2016) McCray, R., & Fransson, C. 2016, ARA&A, 54, 19, doi: 10.1146/annurev-astro-082615-105405
  * Miceli et al. (2019) Miceli, M., Orlando, S., Burrows, D. N., et al. 2019, Nature Astronomy, 3, 236, doi: 10.1038/s41550-018-0677-8
  * Ono et al. (2020) Ono, M., Nagataki, S., Ferrand, G., et al. 2020, ApJ, 888, 111, doi: 10.3847/1538-4357/ab5dba
  * Orlando et al. (2015) Orlando, S., Miceli, M., Pumo, M. L., & Bocchino, F. 2015, ApJ, 810, 168, doi: 10.1088/0004-637X/810/2/168
  * Orlando et al. (2019) Orlando, S., Miceli, M., Petruk, O., et al. 2019, A&A, 622, A73, doi: 10.1051/0004-6361/201834487
  * Orlando et al. (2020) Orlando, S., Ono, M., Nagataki, S., et al. 2020, A&A, 636, A22, doi: 10.1051/0004-6361/201936718
  * Park et al. (2005) Park, S., Zhekov, S. A., Burrows, D. N., & McCray, R. 2005, ApJ, 634, L73, doi: 10.1086/498848
  * Ravi et al. (2021) Ravi, A. P., Park, S., Zhekov, S. A., et al. 2021, ApJ, 922, 140, doi: 10.3847/1538-4357/ac249a
  * Sapienza et al. (2024) Sapienza, V., Miceli, M., Bamba, A., et al. 2024, ApJ, 961, L9, doi: 10.3847/2041-8213/ad16e3
  * Sugerman et al. (2005) Sugerman, B. E. K., Crotts, A. P. S., Kunkel, W. E., Heathcote, S. R., & Lawrence, S. S. 2005, ApJS, 159, 60, doi: 10.1086/430408
  * Sun et al. (2021) Sun, L., Vink, J., Chen, Y., et al. 2021, ApJ, 916, 41, doi: 10.3847/1538-4357/ac033d
  * Zhekov et al. (2009) Zhekov, S. A., McCray, R., Dewey, D., et al. 2009, ApJ, 692, 1190, doi: 10.1088/0004-637X/692/2/1190


#  The efficacy of antiviral drug, HIV viral load and the immune response

Mesfin Asfaw Taye  West Los Angeles College, Science Division
9000 Overland Ave, Culver City, CA 90230, USA

###### Abstract

Developing antiviral drugs is an exigent task since viruses mutate to overcome
the effect of antiviral drugs. As a result, the efficacy of most antiviral
drugs is short-lived. To include this effect, we modify the Neumann and Dahari
model. Considering the fact that the efficacy of the antiviral drug varies in
time, the differential equations introduced in the previous model systems are
rewritten to study the correlation between the viral load and antiviral drug.
The effect of antiviral drug that either prevents infection or stops the
production of a virus is explored. First, the efficacy of the drug is
considered to decreases monotonously as time progresses. In this case, our
result depicts that when the efficacy of the drug is low, the viral load
decreases and increases back in time revealing the effect of the antiviral
drugs is short-lived. On the other hand, for the antiviral drug with high
efficacy, the viral load, as well as the number of infected cells,
monotonously decreases while the number of uninfected cells increases. The
dependence of the critical drug efficacy on time is also explored. Moreover,
the correlation between viral load, the antiviral drug, and CTL response is
also explored. In this case, not only the dependence for the basic
reproduction ratio on the model parameters is explored but also we analyze the
critical drug efficacy as a function of time. We show that the term related to
the basic reproduction ratio increases when the CTL response step up. A simple
analytically solvable mathematical model is also presented to analyze the
correlation between viral load and antiviral drugs.

###### pacs:

Valid PACS appear here

##  I Introduction

Viruses are tiny particles that occupy the world and have property between
living and non-living things. As they are not capable of reproducing, they
rely on the host cells to replicate themselves. To gain access, the virus
first binds and intrudes into the host cells. Once the virus is inside the
cell, it releases its genetic materials into the host cells. It then starts
manipulating the cell to multiply its viral genome. Once the viral protein is
produced and assembled, the new virus leaves the cell in search of other host
cells. Some viruses can also stay in the host cells for a long time as a
latent or chronic state. The genetic information of a virus is stored either
in form of RNA or DNA. Depending on the type of virus, the host cell type also
varies. For instance, the Human Immunodeficiency Virus (HIV) (see Fig. 1 muu1
) directly affects Lymphocytes. Lymphocytes can be categorized into two main
categories: the B and T cells. The B cells directly kill the virus by
producing a specific antibody. T cells on the other hand can be categorized as
killer cells (CD 8) and helper cells (CD 4). Contrary to CD 8, CD 4 gives only
warning so that the cells such as CD 8 and B cells can directly kill the virus
mes5 ; muu2 . Although most of these viruses contain DNA as genetic material,
retroviruses such as HIV store their genetic materials as RNA. These viruses
translate their RNA into DNA using an enzyme called reverse transcriptase
during their life cycle. In the case of HIV, once HIV infects the patient, a
higher viral load follows for the first few weeks, and then its replication
becomes steady for several years. As a result, the CD 4 (which is the host
cell for HIV) decreases considerably. When the CD 4 cells are below a certain
threshold, the patient develops AIDS.

Figure 1:  (Color online) Schematic diagram for HIV virion muu1 .

To tackle the spread of virulent viruses such as HIV, discovering potent
antiviral drugs is vital. However, developing antiviral drugs is an exigent
task since viruses mutate to overcome the effect of antiviral drugs because of
this only a few antiviral drugs are currently available. Most of these drugs
are developed to cure HIV and herpes virus. The fact that viruses are obligate
parasites of the cells makes drug discovery complicated since the drug’s
adverse effects directly affect the host cells. Many medically important
viruses are also virulent and hence they cannot be propagated or tested via
animal models. This in turn forces the drug discovery to be lengthy. Moreover,
unlike other antimicrobial drugs, antiviral drugs have to be 100 percent
potent to completely avoid drug resistance. In other words, if the drug
partially inhibits the replication of the virus, through time, the number of
resistant viruses will dominate the cell culture. All of the above factors
significantly hinder drug discovery. Furthermore, even a potent antiviral drug
does not guarantee a cure if the acute infection is already established.

Understanding the dynamics of the virus in vivo or vitro is crucial since
viral diseases are the main global health concern. For instance, recent
outbreaks of viral diseases such as COVID-19 not only cost trillion dollars
but also killed more than 217,721 people in the USA alone. To control such a
global pandemic, developing an effective therapeutic strategy is vital.
Particularly, in the case of virulent viruses, mathematical modeling along
with the antiviral drug helps to understand the dynamics of the virus in vivo
mes2 . The pioneering mathematical models on the Human Immunodeficiency virus
depicted in the works mes1 ; mes3 ; mes4 ; mu1 ; mu2 ; mu3 ; mu4 ; mu5 ; mu6 ;
mu7  shed light regarding the host-virus correlation. Latter these model
systems are modified by Neumann $et.$ $al.$ mes1 ; mu8  to study the dynamics
of HCV during treatment. To study the dependence of uninfected cells, infected
cells, and virus load on model parameters, Neumann proposed three differential
equations. More recently, to explore the observed HCV RNA profile during
treatment, Dahari $et.$ $al.$ mes1 ; mu9  extended the original Neumann model.
Their findings disclose that critical drug efficacy plays a critical role.
When the efficacy is greater than the critical value, the HCV will be cleared.
On the contrary, when the efficacy of the drug is below the critical
threshold, the virus keeps infecting the target cells.

As discussed before, the effect of antiviral drugs is short-lived since the
virus mutates during the course of treatment. To include this effect, we
modify Neumanny and Dahariy models. Considering the fact that the efficacy of
the antiviral drug decreases in time, we rewrite the three differential
equations introduced in the previous model systems. The mathematical model
presented in this work analyzes the effect of an antiviral drug that either
prevents infection ($e_{k}$) or stops the production of virus ($e_{p})$.
First, we consider a case where the efficacy of the drug decreases to zero as
time progresses and we then discuss the case where the efficacy of the drug
decreases to a constant value as time evolves maintaining the relation
$e_{P}=e_{p}^{\prime}(1+e^{-rt})/m$ and $e_{r}=e_{r}^{\prime}(1+e^{-rt})/m$.
Here $r$, $e_{k}^{\prime}$ and $e_{p}^{\prime}$ measure the ability of
antiviral drug to overcome drug resistance. When $r$ tends to increase, the
efficacy of the drug decreases. The results obtained in this work depict that
for large $r$, the viral load decreases and increases back as the antiviral
drug is administered showing the effect of antiviral drugs is short-lived. On
the other hand, for small $r$, the viral load, as well as the number of
infected cells monotonously decreases while the host cell increases. The
dependence of the critical drug efficacy on time is also explored.

The correlation between viral load, antiviral therapy, and cytotoxic
lymphocyte immune response (CTL) is also explored. Not only the dependence for
the basic reproduction ratio on the model parameters is explored but also we
find the critical drug efficacy as a function of time. The basic reproduction
ratio increases when the CTL response decline. When the viral load inclines,
the CTL response step up. We also present a simple analytically solvable
mathematical model to address the correlation between drug resistance and
antiviral drugs.

The rest of the paper is organized as follows: in Section II, we explore the
correlation between antiviral treatment and viral load. In Section III the
relation between viral load, antiviral therapy, and the CTL immune response is
examined. A simple analytically solvable mathematical model that addresses the
correlation between drug resistance and viral load is presented in section IV.
Section V deals with summary and conclusion.

##  II  The relation between antiviral drug and virus load

In the last few decades, mathematical modeling along with antiviral drugs
helps to develop a therapeutic strategy. The first model that describes the
dynamics of host cells $x$, virus load $v$, and infected cells $y$ as a
function of time $t$ was introduced in the works mu1 ; mu2 ; mu3 ; mu4 ; mu5 .
Accordingly, the dynamics of the host cell, infected cell, and virus is
governed by

$\displaystyle{\dot{x}}$ $\displaystyle=$ $\displaystyle\lambda-dx-\beta xv,$
$\displaystyle{\dot{y}}$ $\displaystyle=$ $\displaystyle\beta xv-ay,$
$\displaystyle{\dot{v}}$ $\displaystyle=$ $\displaystyle ky-uv.$ (1)

The host cells are produced at rate of $\lambda$ and die naturally at a
constant rate $d$ with a half-life of $x_{t_{1\over 2}}={ln(2)\over d}$. The
target cells become infected at a rate of $\beta$ and die at a rate of $a$
with a corresponding half-life of $y_{t_{1\over 2}}={ln(2)\over a}$. On the
other hand, the viruses reproduce at a rate of $k$ and die with a rate $u$
with a half-life of $v_{t_{1\over 2}}={ln(2)\over u}$ mu10 ; mu11 ; mu12 ;
mu14 . In this model, only the interaction between the host cells and viruses
is considered neglecting other cellular activities. The host cells, the
infected cells, and the viruses have a lifespan of $1/d$, $1/a$, and $1/u$,
respectively. During the lifespan of a cell, one infected cell produces
$N=u/a$ viruses on average muu2 ; mu12 ; mu14 .

The capacity for the virus to spread can be determined via the basic
reproductive ratio

$\displaystyle R_{0}={\lambda\beta k\over adu}.$ (2)

Whenever $R_{0}>1$, the virus spreads while when $R_{0}<1$ the virus will be
cleared by the host immune system muu2 .

To examine the dependence of uninfected cells, infected cells and virus load
on the system parameters during antiviral treatment, the above model system
(Eq. (1)) was modified by Neumann $et.$ $al.$ mu8  and Dahari $et.$ $al.$ mu9
. The modified mathematical model presented in those works analyzes the effect
of antiviral drugs that either prevents infection of new cells ($e_{k}$) or
stops production of the virion ($e_{p})$. In this case, the above equation can
be remodified to include the effect of antiviral drugs as

$\displaystyle{\dot{x}}$ $\displaystyle=$ $\displaystyle\lambda-
dx-(1-e_{k})\beta xv$ $\displaystyle{\dot{y}}$ $\displaystyle=$
$\displaystyle(1-e_{k})\beta xv-ay$ $\displaystyle{\dot{v}}$ $\displaystyle=$
$\displaystyle(1-e_{p})ky-uv$ (3)

where the terms $e_{k}$ and $e_{p}$ are used when the antivirus blocks
infection and virion production, respectively. For instance, $e_{p}=0.8$
indicates that the drug has efficacy in blocking virus production by $80$
percent. The antiviral drug such as protease inhibitor inhibits the infected
cell from producing the right gag protein as a result the virus becomes
noninfectious. A drug such as a reverse transcriptase inhibitor prohibits the
infection of new cells.

Moreover, the results obtained in the last few decades depict that, usually
HIV patient shows high viral load in the first few weeks of infection. As a
result, the viral load becomes the highest then it starts declining for a few
weeks. The viruses then keep replicating for many years until the patient
develops AIDS. Since virus replication is prone to errors, the virus often
develops drug resistance; HIV mutates to become drug-resistant. Particularly
when the antiviral drug is administered individually, the ability of the virus
to develop drug resistance steps up. However, a triple-drug therapy which
includes one protease inhibitor combined with two reverse transcriptase
inhibitors helps to reduce the viral load for many years muu2 .

Since the antiviral drugs are sensitive to time, to include this effect, next,
we will modify the Neumann and Dahari model.

Case one.—As discussed before, the effect of antiviral drugs is short-lived
since the virus mutates once the drug is administrated. To include this
effect, we modify the above equation by assuming that
$e_{P}=e_{p}^{\prime}e^{-rt}$ and $e_{r}=e_{r}^{\prime}e^{-rt}$. The efficacy
of the drugs declines exponentially as time progresses. The decaying rate
aggravates when $r$ tends to increase. Hence let us rewrite Eq. (3) as

$\displaystyle{\dot{x}}$ $\displaystyle=$ $\displaystyle\lambda-
dx-(1-e_{k}^{\prime}e^{-rt})\beta xv$ $\displaystyle{\dot{y}}$
$\displaystyle=$ $\displaystyle(1-e_{k}^{\prime}e^{-rt})\beta xv-ay$
$\displaystyle{\dot{v}}$ $\displaystyle=$
$\displaystyle(1-e_{p}^{\prime}e^{-rt})ky-uv.$ (4)

The term related to the reproductive ratio is given as

$\displaystyle R_{1}={\lambda\beta k\over
adu}(1-e_{k}^{\prime})(1-e_{p}^{\prime}).$ (5)

When $R_{1}<1$, the antivirus drug is capable of clearing the virus and if
$R_{1}>1$, the virus tends to spread. At steady state,

$\displaystyle\overline{x}$ $\displaystyle=$ $\displaystyle{au\over\beta k}$
$\displaystyle\overline{y}$ $\displaystyle=$ $\displaystyle{\beta\lambda
k-adu\over a\beta k}$ $\displaystyle\overline{v}$ $\displaystyle=$
$\displaystyle{\beta\lambda k-adu\over a\beta u}.$ (6)

As one can note that, when $t\to\infty$, $e_{p}\to 0$ and $e_{k}\to 0$. At
steady state, only one newly infected cell arises from one infected cell muu2
.

Figure 2:  (Color online) (a) The number of host cells $x$ as a function of
time (days). (b) The number of infected cells as function of time (days). (c)
The virus load as a function of the time (days). In the figure, we fix
$\lambda=10^{7}$, $d=0.1$, $a=0.5$, $\beta=2X10^{-9}$, $k=1000.0$, $u=5.0$,
$e_{p}=0.5$, $e_{k}=0.5$ and $r=0.06$.

Let us next explore how the number of host cells $x$, the number of infected
cells $y$, and the viral load $v$ behave as a function of time by exploiting
Eq. (4) numerically. From now on, all of the physiological parameters are
considered to vary per unit time (days). Figure 1 depicts the plot of the
number of host cells $x$, the number of infected cells $y$ and the number of
virus as function of time (days) for parameter choice of $\lambda=10^{7}$,
$d=0.1$, $a=0.5$, $\beta=2X10^{-9}$, $k=1000.0$, $u=5.0$, $e_{p}=0.5$,
$e_{k}=0.5$ and $r=0.06$. The figure depicts that in the presence of an
antiviral drug, the number of $CD_{4}$ cells increases and attains a maximum
value. The cell numbers then decrease and saturate to a constant value. The
number of infected cells decreases and saturates to a constant value. On the
other hand, the viral load decreases as the antiviral takes an effect.
However, this effect is short-lived since the the viral load increases back as
the viruses develop drug resistance.

Figure 3:  (Color online) (a) The number of host cells $x$ as a function of
time (days). (b) The number of infected cells as function of time (days). (c)
The plasma virus load as a function of time (days). In the figure, we fix
$\lambda=10^{7}$, $d=0.1$, $a=0.5$, $\beta=2X10^{-9}$, $k=1000.0$, $u=5.0$,
$e_{p}=0.9$, $e_{k}=0.9$ and $r=0.0001$.

When $r$ is small, the ability of the antiviral drug to overcome drug
resistance increases. As depicted in Fig. (2), for very small $r$, the host
cells increase in time, and at a steady state, the cells saturate to a
constant value. On the contrary, the infected cells as well as the plasma
virus load monotonously decrease as time progresses. The figure is plotted by
fixing $\lambda=10^{7}$, $d=0.1$, $a=0.5$, $\beta=2X10^{-9}$, $k=1000$,
$u=5.0$, $e_{p}=0.9$, $e_{k}=0.9$ and $r=0.0001$. These all results indicate
that when combined drugs are administrated, the viral load is significantly
reduced depending on the initial viral load.

Figure 4:  (Color online) (a) The number of host cells $x$ as a function of
days. (b) The number of infected cells as function of days. (c) The virus load
as a function of the days. In the figure, we fix $\lambda=10^{7}$, $d=0.1$,
$a=0.5$, $\beta=2X10^{-9}$, $k=1000.0$, $u=5.0$, $e_{p}=0.9$, $e_{k}=0.9$ and
$r=20.0$.

Figure 3 is plotted by fixing $\lambda=10^{7}$, $d=0.1$, $a=0.5$,
$\beta=2X10^{-9}$, $k=1000$, $u=5.0$, $e_{p}=0.9$, $e_{k}=0.9$ and $r=20.0$.
The figure exhibits that for large $r$, the $CD_{4}$ cells decrease and
exhibit a local minima. As time progresses, the number of cells increases and
saturates to a constant value. On the other hand, the number of infected cells
$y$ and the viral load $v$ decreases and saturates to considerably large value
as time progresses. Figure 3 also depicts that when the drug is unable to
control the infection in a very short period of time, the number of drug
resistant viruses steps up.

Case two.— In the previous case, the efficacy of the drug is considered to
decrease monotonously as time progresses. In this section, the efficacy of the
drug is assumed to decrease to a constant value as time increases maintaining
the relation $e_{P}=e_{p}^{\prime}(1+e^{-rt})/m$ and
$e_{r}=e_{r}^{\prime}(1+e^{-rt})/m$. The dynamics of host cells, infected
cells, and viral load is governed by the equation

$\displaystyle{\dot{x}}$ $\displaystyle=$ $\displaystyle\lambda-
dx-(1-e_{r}^{\prime}(1+e^{-rt})/m)\beta xv$ $\displaystyle{\dot{y}}$
$\displaystyle=$ $\displaystyle(1-e_{r}^{\prime}(1+e^{-rt})/m)\beta xv-ay$
$\displaystyle{\dot{v}}$ $\displaystyle=$
$\displaystyle(1-e_{p}^{\prime}(1+e^{-rt})/m)ky-uv.$ (7)

After some algebra, the term related to the basic reproductive ratio reduces
to

$\displaystyle R_{1}={\lambda\beta k\over
adum^{2}}\left(m-2e_{k}^{\prime}\right)\left(m-2e_{p}^{\prime}\right).$ (8)

As one can see from Eq. (8) that as $m$ steps up, the drug losses its potency
and as a result $R_{1}$ increases. When $R_{1}<1$, the antivirus drug
treatment is successful and this occurs for large values of $e_{p}$ and
$e_{k}$. When $R_{1}>1$, the virus overcomes the antivirus treatment.

At equilibrium, one finds

$\displaystyle\overline{x}$ $\displaystyle=$ $\displaystyle{aum^{2}\over\beta
k(m-e_{r}^{\prime})(m-e_{p}^{\prime})}$ $\displaystyle=$
$\displaystyle{\lambda\over
dR_{0}}{m^{2}\over(m-e_{r}^{\prime})(m-e_{p}^{\prime})}$
$\displaystyle\overline{y}$ $\displaystyle=$ $\displaystyle{\lambda\over
a}+{dum^{2}\over\beta k(e_{r}^{\prime}-m)(m-e_{p}^{\prime})}$ $\displaystyle=$
$\displaystyle\left(R_{0}-{m^{2}\over(m-e_{r}^{\prime})(m-e_{p}^{\prime})}\right){du\over\beta
k}$ $\displaystyle\overline{v}$ $\displaystyle=$ $\displaystyle{dm\over(\beta
e_{k}^{\prime}-\beta m})+{\lambda k(m-e_{p}^{\prime})\over(amu)}$ (9)
$\displaystyle=$ $\displaystyle\left({R_{0}(m-e-p^{\prime})\over
m}-{m\over(m-e_{k}^{\prime})}\right){d\over\beta}.$

The case $R_{0}\gg 1$ indicates that the equilibrium abundance of the
uninfected cells is much less than the number of uninfected cells before
treatment. When the drug is successful, (large values of $e_{p}^{\prime}$ or
$e_{r}^{\prime}$), the equilibrium abundance of the uninfected cells
increases. On the contrary, for a highly cytopathic virus ( $R_{1}\gg 1$ ),
the number of infected cells, as well as the viral load steps up. When
$e_{p}^{\prime}$ and $e_{r}^{\prime}$ increase, the equilibrium abundance of
infected cells as well as viral load decreases.

In general for large $R_{0}$, Eq. (9) converges to

$\displaystyle\overline{y}$ $\displaystyle=$ $\displaystyle{\lambda\over a}$
$\displaystyle\overline{v}$ $\displaystyle=$ $\displaystyle{\lambda
k(m-e_{p}^{\prime})\over(amu)}.$ (10)

Clearly $\overline{v}$ decreases as $e_{p}^{\prime}$ and $e_{r}^{\prime}$
increase.

The overall efficacy can be written as mu9

$\displaystyle 1-e=\left(1-{e_{r}^{\prime}\left(1+e^{-rt}\right)\over
m}\right)\left(1-{e_{p}^{\prime}\left(1+e^{-rt}\right)\over m}\right)$ (11)

where $0<e_{r}^{\prime}<1$ and $0<e_{p}^{\prime}<1$. At steady state
$1-e=\left(1-{e_{r}^{\prime}\over m}\right)\left(1-{e_{p}^{\prime}\over
m}\right)$. The transcritical bifurcation point (at steady state) can be
analyzed via Eq. (9) and after some algebra we find

$\displaystyle 1-e=\left(1-{e_{r}^{\prime}\over
m}\right)\left(1-{e_{p}^{\prime}\over m}\right)={adum\over\lambda\beta
k}={x_{1}\over x_{0}}$ (12)

where $x_{0}={\lambda\over d}$ denotes the number of uninfected host cells
before infection and $x_{1}={aum\over\beta k}$ designates the number of
uninfected cells in the chronic case. This implies the critical efficacy is
given as $e_{c}=1-{x_{1}\over x_{0}}=1-{adum\over\lambda\beta k}$.

To write the overall efficacy as a function of time, for simplicity, let us
further assume that $e_{r}^{\prime}=e_{p}^{\prime}$. In this case, Eq. (12)
can be rewritten as

$\displaystyle{e_{r}^{\prime}\over m}=1\pm\sqrt{{adum\over\lambda\beta k}}$
(13)

and hence

$\displaystyle 1-e=\left(1-\left(1\pm\sqrt{{adum\over\lambda\beta
k}}\right)\left(1+e^{-rt}\right)\right)^{2}.$ (14)

From Eq. (14), one finds

$\displaystyle e_{c}=1-\left(1-\left(1\pm\sqrt{{adum\over\lambda\beta
k}}\right)\left(1+e^{-rt}\right)\right)^{2}.$ (15)

The critical efficacy serves as an alternative way of determining whether
antiviral treatment is successful or not. When $e>e_{c}$, the antiviral clears
the infection and if $e<e_{c}$, the virus replicates.

##  III  The correlation between antiviral drug, immune response and virus
load

The basic mathematical model that specifies the relation between the immune
response, antiviral drug, and viral load is given by

$\displaystyle{\dot{x}}$ $\displaystyle=$ $\displaystyle\lambda-
dx-(1-e_{r}^{\prime}(1+e^{-rt})/m)\beta xv$ $\displaystyle{\dot{y}}$
$\displaystyle=$ $\displaystyle(1-e_{r}^{\prime}(1+e^{-rt})/m)\beta xv-ay-pyz$
$\displaystyle{\dot{v}}$ $\displaystyle=$
$\displaystyle(1-e_{p}^{\prime}(1+e^{-rt})/m)ky-uv$ $\displaystyle{\dot{z}}$
$\displaystyle=$ $\displaystyle c-bz.$ (16)

Once again the terms $x$, $y$, and $v$ denote the uninfected cells, infected
cells, and the viral load. The term $z$ denotes the CTL response and the CTL
die at a rate of $b$ and produced at a rate of $c$. The term CTL is defined as
cytotoxic lymphocytes that have responsibility for killing the infected cells.
The term related to the basic reproduction rate is given as

$\displaystyle R_{1}={\lambda\beta k\over(a+{cp\over
b})dum^{2}}(m-2e_{k}^{\prime})(m-2e_{p}^{\prime}).$ (17)

It vital to see that when $R_{0}>0$ the virus becomes successful to persist an
infection which triggers an immune responses $z$. As long as the coordination
between the immune response and the antiviral drug treatment is strong enough,
the virus will be cleared $R_{1}<1$. As it can be clearly seen from Eq. (17),
when the CTL response step up, $R_{1}$ declines as expected.

The equilibrium abundance of the host cells, infected cells, viral load, and
CTL response can be given as

$\displaystyle\overline{x}$ $\displaystyle=$
$\displaystyle{m^{2}u(ab+cp)\over\beta
bk(m-e_{k}^{\prime})(m-e_{p}^{\prime})}$ $\displaystyle=$
$\displaystyle{\lambda\over
adR_{0}}{(ab+cp)m^{2}\over(m-e_{r}^{\prime})(m-e_{p}^{\prime})}$
$\displaystyle\overline{y}$ $\displaystyle=$
$\displaystyle{b\lambda\over(ab+cp)}+{dum^{2}\over\beta k(e_{r}-m)(m-e_{p})}$
$\displaystyle=$
$\displaystyle\left({abR_{0}\over(ab+cp)}-{m^{2}\over(m-e_{r}^{\prime})(m-e_{p}^{\prime})}\right){du\over\beta
k}$ $\displaystyle\overline{v}$ $\displaystyle=$ $\displaystyle{dm\over(\beta
e_{k}-\beta m})+{b\lambda k(m-e_{p})\over mu(ab+cp)}$ $\displaystyle=$
$\displaystyle\left({bR_{0}(m-e-p^{\prime})\over
m(ab+cd)}-{m\over(m-e_{k}^{\prime})}\right){d\over\beta}$
$\displaystyle\overline{z}$ $\displaystyle=$ $\displaystyle{c\over b}.$ (18)

Exploiting Eq. (18), one can deduce when $R_{0}\gg 1$, the equilibrium
abundance of the uninfected cells becomes much lower in comparison to the
number of cells before treatment. The equilibrium abundance of the uninfected
cells steps up when there is CTL response (when $c$ and $p$ increase) or when
the antiviral treatment is successful (when $e_{p}^{\prime}$ and
$e_{r}^{\prime}$ increase). On the contrary, $\overline{y}$ and $\overline{v}$
decline whenever there is a strong CTL response or when the antiviral
treatment is successful.

In Fig. 5a , we plot the phase diagram for a regime $R_{1}<1$ (shaded region)
in the phase space of $e_{k}$ and $e_{p}$. In Figure 5b, the phase diagram for
a regime $R_{1}<1$ (shaded region) in the phase space of $m$ and $e_{p}=e_{k}$
is plotted. In the figure, we fix $\lambda=10^{7}$, $d=0.1$, $a=0.5$,
$\beta=2X10^{-9}$, $k=1000$, $u=5$, $m=2$, $r=0.0001$, $p=1.0$, $b=0.5$ and
$c=2.0$.

Furthermore, the transcritical bifurcation point can be analyzed via Eq. (18)
and after some algebra (at steady state) we find

$\displaystyle
1-e=\left(1-e_{r}^{\prime}/m\right)\left(1-e_{p}^{\prime}/m\right)={(ab+cp)dum\over\lambda\beta
kb}={x_{1}\over x_{0}}$ (19)

where $x_{0}={\lambda\over d}$ denotes the number of uninfected host cells
before infection and $x_{1}={(ab+cp)um\over\beta kb}$ designates the number of
uninfected cells in a chronic case. This implies the critical efficacy is
given as $e_{c}=1-{x_{1}\over x_{0}}=1-{(ab+cp)umd\over\beta kb\lambda}$.

Figure 5:  (Color online) (a) The phase diagram for a regime $R_{0}<1$ in the
phase space of $e_{k}$ and $e_{p}$. (b) The phase diagram for a regime
$R_{0}<1$ in the phase space of $m$ and $e_{p}=e_{k}$. In the figure, we fix
$\lambda=10^{7}$, $d=0.1$, $a=0.5$, $\beta=2X10^{-9}$, $k=1000$, $u=5$, $m=2$,
$r=0.0001$, $p=1.0$, $b=0.5$ and $c=2.0$.

Assuming $e_{r}^{\prime}=e_{p}^{\prime}$, Eq. (19) can be rewritten as

$\displaystyle e_{r}^{\prime}/m=1\pm\sqrt{{(ab+cp)dum\over\lambda\beta kb}}$
(20)

and the overall efficacy as a function of time is given by

$\displaystyle\left(1-e\right)=\left(1-\left(1\pm\sqrt{{(ab+cp)dum\over\lambda\beta
kb}}\right)\left(1+e^{-rt}\right)\right)^{2}.$ (21)

From Eq. (21), one finds

$\displaystyle
e_{c}=1-\left(1-\left(1\pm\sqrt{{\left(ab+cp\right)dum\over\lambda\beta
kb}}\right)\left(1+e^{-rt}\right)\right)^{2}.$ (22)

Once again, the infection will be cleared when $e>e_{c}$, and the virus
replicates as long as $e<e_{c}$.

Figure 6:  (Color online) (a) The number of host cells $x$ as a function of
days. (b) The number of infected cells as function of days (c) The virus load
as a function of the days. In the figure, we fix $\lambda=10^{7}$, $d=0.1$,
$a=0.5$, $\beta=2X10^{-9}$, $k=1000$, $u=5.0$, $e_{p}^{\prime}=0.9$,
$e_{k}^{\prime}=0.9$, $r=0.1$, $p=10.0$, $b=0.5$ and $c=1.0$.

In order to get a clear insight, let us explore Eq. (16). In Fig. 6, the
number of host cells $x$, the number of infected cells and the viral load as a
function of time is plotted. In the figure, we fix $\lambda=10^{7}$, $d=0.1$,
$a=0.5$, $\beta=2X10^{-9}$, $k=1000$, $u=5.0$, $e_{p}^{\prime}=0.9$,
$e_{k}^{\prime}=0.9$, $r=0.1$, $p=10.0$, $b=0.5$ and $c=1.0$. For such
parameter choice, $R_{0}\gg 1.0$ and $R_{1}\ll 1.0$ revealing that initially,
the virus establishes an infection but latter the antiviral drug and CTL
response collaborate to clear the infection. As a result, the number of host
cells increases while the viral load as well as the infected cells decreases.
Since the virus is responsible for initiating the CTL response, as the viral
load declines, the CTL response step down.

##  IV The dynamics of mutant virus in the presence of antiviral drugs

As discussed before, the fact that viruses are an obligate parasite of the
cells forces drug discovery to be complicated as the drug’s adverse effects
directly affect the host cells. Many medically important viruses are also
virulent and hence they cannot be propagated or tested via animal models. This
in turn forces the drug discovery to be lengthy. Moreover, unlike other
antimicrobial drugs, antiviral drugs have to be 100 percent potent to
completely avoid drug resistance. In other words, if the drug partially
inhibits the replication of the virus, through time, the number of the
resistant virus will dominate the cell culture.

To discuss the dynamics of the mutant virus, let us assume that the virus
mutates (when a single drug is administered ) by changing one base every
$10^{4}$ viruses. If $10^{11}$ viruses are produced per day, then this results
in $10^{7}$ mutant viruses. On the contrary, when two antiviral drugs are
administrated, $10^{3}$ mutant viruses will be produced. In the case of triple
drug therapy, no mutant virus is produced. To account for this effect, let us
remodify the rate of mutant virus production per day as

$\displaystyle k=10^{11-4s},{s=1,2,3}$ (23)

where the variable $s=1,~{}2,~{}3$ corresponds to a single, double, and triple
drug therapy, respectively.

Here we assume only the rate of mutant virus production $k$ determines the
dynamics and $10^{11}$ viruses are produced per day. To get an instructive
analytical solution regarding the relation between antiviral drug and viral
load, let us solve the differential equation

$\displaystyle{\dot{v}}$ $\displaystyle=$ $\displaystyle k-uv$ (24)

neglecting the effect of uninfected and infected host cells. Here the mutant
virus produced at rate of $k$ and die with the rate of $u$. The solution for
the above equation is given as

$\displaystyle v(t)={e^{-ut}(-k+ke^{ut}+uN)\over u}.$ (25)

Whenever ${k\over u}>1$ the virus spreads and when ${k\over u}<1$, the
antiviral is capable of eliminating the virus.

Exploiting Eq. (25) one can comprehend that, in the case of single therapy,
the virus load decreases during the course of treatment. As time progresses,
the viral load increases back due to the emergence of drug resistance (see
Fig. 8a). In the case of double drug therapy, as shown in Fig. 8b, the viral
load decreases but relapses back as time progresses. When triple drugs are
administered, the viral replication becomes suppressed as depicted in Fig 8c.
The readers should understand that the triple drug therapy does not guarantee
a cure. If the patient halts his or her therapy, the viral replications will
resume because of the latent and chronic infected cells. At steady state, we
get

$\displaystyle\overline{v}={k\over u}$ (26)

At equilibrium, the viral load spikes as $k$ increases and it decreases as $u$
steps up.

Figure 7:  (Color online) The virus load as a function of time. In the figure,
we fix $u=1.0$. (a) Single drug therapy $s=1$. (b) Double drug therapy $s=2$.
(c) Triple drug therapy $s=3$.

##  V Summary and conclusion

Developing antiviral drugs is challenging but an urgent task since outbreaks
of viral diseases not only killed several people but also cost trillion
dollars worldwide. The discovery of new antiviral drugs together with emerging
mathematical models helps to understand the dynamics of the virus in vivo. For
instance, the pioneering mathematical models on HIV shown in the works mu1 ;
mu2 ; mu3 ; mu4 ; mu5 ; mu6 ; mu7  disclose the host-virus correlation.
Moreover, to study the correlation between, antiviral drugs and viral load, an
elegant mathematical model was presented in the works mu8 ; mu9 .

Due to the emergence of drug resistance, the efficiency of antiviral drugs is
short-lived. To study this effect, we numerically study the dynamics of the
host cells and viral load in the presence of an antiviral drug that either
prevents infection ($e_{k}$) or stops the production of virus ($e_{p}$). For
the drug whose efficacy depends on time, we show that when the efficacy of the
drug is low, the viral load decreases and increases back in time revealing the
effect of the antiviral drugs is short-lived. On the contrary, for the
antiviral drug with high efficacy, the viral load, as well as the number of
infected cells, monotonously decreases while the number of uninfected cells
increases. The dynamics of critical drug efficacy on time is also explored.
Furthermore, the correlation between viral load, an antiviral drug, and CTL
response is also explored. Not only the dependence for the basic reproduction
ratio on the model parameters is explored but also we analyze the critical
drug efficacy as a function of time. The term related to the basic
reproduction ratio increases when the CTL response step up. A simple
analytically solvable mathematical model to analyze the correlation between
viral load and antiviral drugs is also presented.

In conclusion, in this work, we present a simple model which not only serves
as a basic tool for better understanding of viral dynamics in vivo and vitro
but also helps in developing an effective therapeutic strategy.

## Acknowledgment

I would like to thank Mulu Zebene and Blaynesh Bezabih for the constant
encouragement.

## References

  * (1) This file is made available by MechESR under the Creative Commons CC 1.0 Universal Public Domain Dedication.Images (Wikimedia Commons).
  * (2) G. Doitsh and W. C. Greene, Cell Host and Microbe 19, 280 (2016).
  * (3) Nowak, M. A., and May, R. (2001). Virus Dynamics: Mathematical Principles of Immunology and Virology. Oxford University Press.
  * (4) S. Wang, Y. Pan, Q. Wang, H. Miao, A. N. Brown and L. Rong, Mathematical Biosciences 328, 108438 (2020).
  * (5) C. Zitzmann and L. Kaderali, Forntiers in Microbiology 9, 1546 (2018).
  * (6) Y. Wang, J. Liu and L. Liu, Advances in Difference Equations 1, 225 (2018).
  * (7) S. S. Chen, C. Y. Cheng and Y. Takeuchi, Journal of Mathematical Analysis and Applications 442, 642 (2016).
  * (8) A. S. Perelson, D. E. Kirschner and R. De Boer , Math. Biosci 114, 81 (1993).
  * (9) A. S. Perelson, , A. U. Neumann, M. Markowitz, J. M Leonard and D. D Ho , Science 271, 1582 (1996).
  * (10) A. S. Perelson, P. Essunger, Y. Cao, M. Vesanen, A. Hurley and K. Saksela, Nature 387, 188 (1997).
  * (11) Ho, D. D., Neumann, A. U., Perelson, A. S., Chen, W., Leonard, J. M., and Markowitz,M. (1995). Rapid turnover of plasma virions and CD4 lymphocytes in HIV-1 infection. Nature 373, 123126. doi: 10.1038/373123a.
  * (12) Bonhoeffer, S., May, R. M., Shaw, G. M., and Nowak, M. A. (1997). Virus dynamics and drug therapy. Proc. Natl. Acad. Sci. U.S.A. 94, 69716976. doi: 10.1073/pnas.94.13.6971.
  * (13) Stafford, M. A., Corey, L., Cao, Y., Daar, E. S., Ho, D. D., and Perelson, A. S. (2000). Modeling plasma virus concentration during primary HIV infection. J. Theor. Biol. 203, 285301. doi: 10.1006/jtbi.2000.1076
  * (14) Wei, X., Ghosh, S. K., Taylor, M. E., Johnson, V. A., Emini, E. A., Deutsch, P., et al. (1995).Viral dynamics in human immunodeficiency virus type 1 infection. Nature 373, 117122. doi: 10.1038/373117a0
  * (15) Neumann, A. U. (1998). Hepatitis C viral dynamics in vivo and the antiviral efficacy of interferon- therapy. Science 282, 103107. doi: 10.1126/science.282.5386.103
  * (16) Dahari, H., Lo, A., Ribeiro, R. M., and Perelson, A. S. (2007a). Modeling hepatitis C virus dynamics: Liver regeneration and critical drug efficacy. J. Theor. Biol. 247, 371381. doi: 10.1016/j.jtbi.2007.03.006
  * (17) Nowak, M. A., and Bangham, C. R. (1996). Population dynamics of immune responses to persistent viruses. Science 272, 7479.
  * (18) Nowak,M. A., Bonhoeffer, S., Hill, A.M., Boehme, R., Thomas, H. C., andMcDade, H. (1996). Viral dynamics in hepatitis B virus infection. Proc. Natl. Acad. Sci. U.S.A. 93, 43984402. doi: 10.1073/pnas.93.9.4398
  * (19) Perelson, A. S. (2002). Modelling viral and immune system dynamics. Nat. Rev. Immunol. 2, 2836. doi: 10.1038/nri700
  * (20) Wodarz, D., and Nowak,M. A. (2002).Mathematical models of HIV pathogenesis and treatment. BioEssays 24, 11781187. doi: 10.1002/bies.10196


# Massively-Parallel Lagrangian Particle Code and Applications

Shaohua Yuan  Mario Zepeda Aguilar  Nizar Naitlho  Roman Samulyak
<EMAIL_ADDRESS>Department of Applied Mathematics and
Statistics, Stony Brook University, Stony Brook, New York 11794, USA.
Computational Science Initiative, Brookhaven National Laboratory, Upton, New
York 11973, USA.

###### Abstract

Massively-parallel, distributed-memory algorithms for the Lagrangian particle
hydrodynamic method [R. Samulyak, X. Wang, H.-C. Chen, Lagrangian particle
method for compressible fluid dynamics, J. Comput. Phys., 362 (2018), 1-19]
have been developed, verified, and implemented. The key component of parallel
algorithms is a particle management module that includes a parallel
construction of octree databases, dynamic adaptation and refinement of
octrees, and particle migration between parallel subdomains. The particle
management module is based on the p4est (parallel forest of k-trees) library.
The massively-parallel Lagrangian particle code has been applied to a variety
of fundamental science and applied problems. A summary of Lagrangian particle
code applications to the injection of impurities into thermonuclear fusion
devices and to the simulation of supersonic hydrogen jets in support of laser-
plasma wakefield acceleration projects has also been presented.

###### keywords:

Lagrangian particle method , Multiphase flows , Parallel algorithms , Parallel
k-tree

††journal: Mechanics Research Communications

###### Contents

  1. 1 Introduction
  2. 2 Lagrangian Particle models and algorithms
    1. 2.1 Main governing equations
    2. 2.2 Massively parallel algorithms
    3. 2.3 Verification and parallel scalability
  3. 3 Applications
    1. 3.1  Summary of thermonuclear fusion applications
    2. 3.2 Supersonic jets for LPWA
  4. 4 Summary and conclusions

##  1 Introduction

Particle-based Lagrangian methods [1, 2, 3] possess a number of unique
features that make them the ideal choice for numerical simulations of certain
classes of problems. By discretizing the medium (compressible fluid, gas,
plasma etc.) with particles, Lagrangian methods achieve a continuum, high-
level of adaptivity: particles move with the flow increasing their density in
compressed regions and automatically provide refinement / coarsening of
numerical resolution where it is needed. Such a continuous adjustment of
resolution eliminates the presence of interfaces between fine and coarse mesh
patches, present in adaptive mesh refinement, that may create a source of
numerical artifacts. Lagrangian particle methods are capable of working with
extremely non-uniform distributions of matter (as typical in astrophysics and
some high-energy-density applications) as they do not need to resolve ”empty”
space between sparse regions of matter. Particle methods can resolve material
boundaries or interfaces of high geometric complexity and they exactly
conserve certain properties of matter in numerical simulations. Like all
numerical methods, they should be used in their best applicability range. For
example, they would under-perform compared to certain spectral methods if a
very high resolution is necessary in a rectangular domain with periodic
boundary conditions.

Smooth particle hydrodynamics (SPH) was the 1st Lagrangian particle method
suitable for 3D computational fluid dynamics. The traditional SPH method [4]
was often limited by a major drawback: its very poor accuracy of discrete
differential operators. It is widely accepted [1, 5, 6] that the traditional
SPH discretization has zeroth-order convergence for widely used kernels.
Numerous methods have been developed to improve the SPH discretization of
differential operators. They include the Reproducing Kernel Particle Methods
(RKPM) [7] and the Normalized Smoothing Function (NSF) methods [8], the
Conservative Reproducing Kernel Smoothed Particle Hydrodynamics (CRKSPH) [2],
and moving-least-squares WENO-SPH schemes [3].

While most of the previous works present modifications and improvements of the
SPH while still preserving the main ideas of SPH and its smooth kernel, the
Lagrangian particle (LP) method [9] significantly deviates from them. It
implements an unsplit Lax–Wendroff type scheme with a limiter. In the location
of points selected by a state estimator, the Lax–Wendroff scheme is blended
with a first-order upwind method to control numerical dispersion. In the LP
method, approximations of spatial derivatives are obtained by employing a
local polynomial fit known also as the generalized finite difference (GFD)
method [10], which is always accurate to a prescribed order. The LP method is
second-order accurate in space and time and it is extendable to higher order
accuracy. In addition to compressible hydrodynamics, the LP code implements
equations of low magnetic Reynolds number magnetohydrodynamics [11] and
granular flow equations [12].

The 1st version of the Lagrangian particle code [9] was written for multi-core
processors that significantly limited its applicability to large scale
scientific simulations. In this paper, we describe a parallel Lagrangian
particle algorithm and the corresponding object-oriented code for massively-
parallel, distributed memory supercomputers. The massively - parallel
Lagrangian particle code is being used for a variety of fundamental science
and applied problems, in particular for the simulation of impurities in
thermonuclear fusion devices. In the last Section, we briefly summarize the
thermonuclear fusion applications and describe in more detail hydrodynamic
simulations for the nozzle design optimization for supersonic hydrogen jets.
Such simulations support experimental work at Brookhaven National Laboratory
on laser-plasma interaction and laser-driven plasma wakefield acceleration.

##  2 Lagrangian Particle models and algorithms

###  2.1 Main governing equations

The core Lagrangian particle (LP) method [9] solves the system of compressible
3D Euler equations in Lagrangian coordinates:

$\displaystyle\frac{d\rho}{dt}=-\rho\nabla\mathbf{u},$ (1)
$\displaystyle\rho\frac{d\mathbf{u}}{dt}=-\nabla P,$ (2)
$\displaystyle\rho\frac{de}{dt}=-P\nabla\cdot\mathbf{u},$ (3) $\displaystyle
P=P(\rho,e),$ (4)

where $d/dt$ is the Lagrangian time derivative, $\mathbf{u}$, $\rho$ and $e$
are the velocity, density and specific internal energy, respectively, and $P$
is the pressure. In addition, various extensions of the system (1 \- 4),
including magnetohydrodynamic equations in the low Magnetic Reynolds number
approximation [11] and granular flow equations [12], have been implemented.
The LP method also implements various models for the equation of state (4)
describing gases, plasmas, and liquids.

The LP method can be summarized as follows. The medium is discretized by
Lagrangian particles, i.e. each Lagrangian particle represents a fluid parcel
characterized by density, pressure, internal energy and other quantities. The
number density of Lagrangain particles is proportional to the mass density of
the medium at each point in space. Differential operators are discretized by
employing a local polynomial fit known also as the generalized finite
difference (GFD) method [10]. The main component of the GFD method, finding
neighborhoods for every Lagrangian particle in a highly non-uniform system, is
implemented via octree data structure construction and search [13]. A second-
order accurate in space and time discretization of the system is obtained
based on the Lax-Wendroff method. The second order method, while used most of
the time in numerical simulations, is prone to numerical dispersion in the
location of discontinuities and high gradients of physics states which leads
to unphysical oscillations of states. To eliminate numerical dispersion, the
2nd order method is blended with a 1st order accurate, diffusive upwind method
via a limiter in locations detected by a state estimator. After updating
states of the Lagrangian particles, they are moved to the next location in
space by the velocity field. These main components of the Lagrangian particle
method are described below.

The main idea of the 1st order upwind method is easy to introduce for a 1D
system of Euler equations which can be written in the following quasi-linear
form using $U=[V$ $u$ $P]^{T}$ as the state vector:

$U_{t}+A(U)U_{x}=0,$ (5) $U=\left(\begin{array}[]{c}V\\\ u\\\
P\end{array}\right),\quad A(U)=V\left(\begin{array}[]{ccc}0&-1&0\\\ 0&0&1\\\
0&K&0\end{array}\right),$ (6)

where

$K=\left(P+\frac{\partial e}{\partial V}\right)\left/\frac{\partial
e}{\partial P}\right.,$ (7)

and $V=1/\rho$ is the specific volume. Performing diagonalization of $A$ as
$A=R\Lambda R^{-1}$, equations (5) and (6) become

$U_{t}+R\Lambda R^{-1}U_{x}=0,$ $R^{-1}U_{t}+\Lambda R^{-1}U_{x}=0,$ (8)

where

$R^{-1}=\left(\begin{array}[]{ccc}1&0&\frac{1}{K}\\\
0&-\frac{1}{2\sqrt{K}}&-\frac{1}{2K}\\\
0&\frac{1}{2\sqrt{K}}&-\frac{1}{2K}\end{array}\right),\,\,R=\left(\begin{array}[]{ccc}1&1&1\\\
0&-\sqrt{K}&\sqrt{K}\\\ 0&-K&-K\end{array}\right),$
$\Lambda=V\left(\begin{array}[]{ccc}0&&\\\ &\sqrt{K}&\\\
&&-\sqrt{K}\end{array}\right)$ (9)

This gives rise to three advection equations for components of $R^{-1}U$ which
are solved by placing numerical stencils in the upwind direction. In terms of
the original state quantities, the upwind equations are

$\displaystyle V_{t}$ $\displaystyle=$
$\displaystyle\frac{V}{2}\left(u_{xr}+u_{xl}\right)-\frac{V}{2\sqrt{K}}\left(P_{xr}-P_{xl}\right),$
(10) $\displaystyle u_{t}$ $\displaystyle=$
$\displaystyle\frac{V\sqrt{K}}{2}\left(u_{xr}-u_{xl}\right)-\frac{V}{2}\left(P_{xr}+P_{xl}\right),$
(11) $\displaystyle P_{t}$ $\displaystyle=$
$\displaystyle-\frac{VK}{2}\left(u_{xr}+u_{xl}\right)+\frac{V\sqrt{K}}{2}\left(P_{xr}-P_{xl}\right).$
(12)

The subscripts $l$ and $r$ to the spatial derivatives in equations (10 \- 12)
indicate that the corresponding terms, in the discrete form, will be computed
using one-sided left (l) or right (r) derivatives. Extension of these
equations to the 3D case can be found in [9].

The first-order upwind scheme is stable but diffusive. We now present a
second-order, unspilt, Lax-Wendroff-type method based on symmetric in space
stencils for Lagrangian particles. As before, we consider a 1D Lagrangian
formulation of Euler equation for polytropic gas. Differentiating (5) - (7)
with respect to $t$ and exchanging the order of the material and spatial
derivatives, we obtain

$\displaystyle V_{tt}+VV_{x}P_{x}+V^{2}P_{xx}=0$ $\displaystyle u_{tt}-\gamma
Vu_{x}P_{x}-\gamma VPu_{xx}=0$ (13) $\displaystyle
P_{tt}-(\gamma^{2}+\gamma)Pu_{x}^{2}-\gamma V_{x}PP_{x}-\gamma VPP_{xx}=0$

Notice that the exchange of the order of the material and spatial derivatives
introduces an additional term

$\frac{Df_{x}}{Dt}=\left(\frac{Df}{Dt}\right)_{x}-u_{x}f_{x}.$ (14)

Using the Taylor series expansion, we obtain the semi-discrete equations

$\frac{V^{n+1}-V^{n}}{\Delta t}=Vu_{x}+\frac{1}{2}\Delta
t(-VV_{x}P_{x}-V^{2}P_{xx})+O(\Delta t^{2})\\\ $ $\frac{u^{n+1}-u^{n}}{\Delta
t}=-VP_{x}+\frac{1}{2}\Delta t(\gamma VPu_{xx}-Vu_{x}P_{x})+O(\Delta t^{2})\\\
$ $\frac{P^{n+1}-P^{n}}{\Delta t}=-\gamma Pu_{x}+\frac{1}{2}\Delta
t(\gamma^{2}Pu_{x}^{2}+\gamma V_{x}PP_{x}+\gamma VPP_{xx}+$ $\gamma
Pu_{x}^{2})+O(\Delta t^{2})$

The derivation in two- and three-dimensional spaces is tedious but
straightforward. The corresponding equations are given in [9]

Numerical approximations of spatial differential operators at locations of
Lagrangian particles are obtained using the Generalized Finite Differences
[10]. Suppose we need to calculate derivatives of a function $U$ at the
location of particle 0 in a 2D domain. In the vicinity of particle $0$, the
function value in the location of particle $i$ can be expressed by the Taylor
series as

$U_{i}=U_{0}+h_{i}\left.\frac{\partial U}{\partial
x}\right|_{0}+k_{i}\left.\frac{\partial U}{\partial y}\right|_{0}+$ (15)
$\frac{1}{2}\left(h^{2}_{i}\left.\frac{\partial^{2}U}{\partial
x^{2}}\right|_{0}+k^{2}_{i}\left.\frac{\partial^{2}U}{\partial
y^{2}}\right|_{0}+2h_{i}k_{i}\left.\frac{\partial^{2}U}{\partial x\partial
y}\right|_{0}\right)+\ldots,$

where $U_{i}$ and $U_{0}$ are the corresponding function values in the
location of particles $i$ and $0$, $h_{i}=x_{i}-x_{0}$, $k_{i}=y_{i}-y_{0}$,
and the derivatives are calculated in the location of particle $0$. Let’s
truncate the Taylor series to a second-order polynomial. In order to
approximate values of the derivatives, we perform a local polynomial fitting
using states on $m>=5$ particles in the vicinity of particle $0$. The
following linear system $Ax=b$

$\begin{bmatrix}h_{1}&k_{1}&\frac{1}{2}h^{2}_{1}&\frac{1}{2}k^{2}_{1}&h_{1}k_{1}\\\
h_{2}&k_{2}&\frac{1}{2}h^{2}_{2}&\frac{1}{2}k^{2}_{2}&h_{2}k_{2}\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
h_{m}&k_{m}&\frac{1}{2}h^{2}_{m}&\frac{1}{2}k^{2}_{m}&h_{m}k_{m}\end{bmatrix}\left[\begin{array}[]{c}U_{x}\\\
U_{y}\\\ U_{xx}\\\ U_{yy}\\\
U_{xy}\end{array}\right]=\left[\begin{array}[]{c}U_{1}-U_{0}\\\ U_{2}-U_{0}\\\
\vdots\\\ U_{m}-U_{0}\end{array}\right],$ (16)

is usually overdetermined. An optimal solution to (16) is a solution $x$ that
minimizes the $L_{2}$ norm of the residual, i.e.,

$min\|Ax-b\|_{2},$ (17)

and the QR decomposition with column pivoting is employed to obtain $x$.

For the second-order method, particle neighborhoods are intended to be
symmetric while including particles from all 8 octants. For the upwind method,
the left and right derivatives are approximated using the left and right semi-
spheres in the corresponding spatial dimension. The 2nd-order method is
blended with the upwind method as

$s_{i}^{n+1}=s_{i}^{n}+\Delta
t\left(UW_{i}-\phi(r_{i})\left(UW_{i}-LW_{i}\right)\right),$ (18)

where $s_{i}^{n}$ is a hydrodynamic state at point $i$ at time $t^{n}$,
$UW_{i}$ and $LW_{i}$ schematically denote the upwind and Lax-Wendroff method
at point $i$, respectively, and $\phi(r_{i})$ is a limiter function based on
values of a state estimator. The state estimator takes into consideration the
quotient of the left $s_{xl}$ and right $s_{xr}$ derivatives and is defined as

$r_{i}=\begin{cases}\infty&s_{xl}\cdot s_{xr}\leq 0\\\
\frac{1}{2}\left(\frac{s_{xr}}{s_{xl}}+\frac{s_{xl}}{s_{xr}}\right)&elsewhere\end{cases}$

The limiter controls the total variation of numerical solutions and reduces
the numerical order of accuracy at extremum points and high solution
gradients. The specific form of the limiter function is often problem
dependent. In some simulations, we use the following empirically designed
limiter function

$\phi(r)=\begin{cases}1.25e^{(1-\sqrt{r})}&r<\frac{3}{2}\\\
1.25e^{-\frac{9}{2}(\sqrt{r}-1)^{2}}&r>\frac{3}{2}\end{cases}$

that effectively eliminates the dispersion while optimizing the convergence
rate.

After the update of states of each Lagrangian particle, particles are advanced
in space by a mixture of the forward and backward Euler schemes:

$\displaystyle\frac{x^{n+1}-x^{n}}{\Delta t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(u^{n}+u^{n+1}\right)$ (19)

Other details of the method can be found in [9].

###  2.2 Massively parallel algorithms

In this Section, we describe main algorithms included in the time step loop
and their implementation for massively parallel supercomputers. The fluid
particles are generated during the initialization. The LP algorithm for each
time step can be presented as the following steps for a 3D problem.

1\. Boundary. Generate boundary particles for given boundary conditions. For
example, mirror particles are used for solid boundary conditions. Inflow and
outflow boundaries can also be set up at this stage by properly generating or
deleting particles.

2\. Octree. Construct an octree for the present distribution of fluid and
boundary particles. Adaptively refine or coarsen the octree until the particle
number in an octant does not exceed a prescribed number.

3\. Neighbor searching. Find neighbors for each particles in six directions
(top, down, front, back, left, right). If the free boundary condition is
applied and neighbors are not complete in certain directions, ghost particles
will be generated to complete the neighborhood.

4\. Stencil. Given the neighbors of each particles, generate stencils for both
the upwind and Lax-Wendroff methods. For the upwind method, neighbors are
organized into six lists in order to compute spatial derivatives in six
directions. Stencils are close to spherical in shape for the Lax-Wendroff
method. A simple sorting of particles in the distance ascending order may lead
to the selection of most particles from one direction while missing
information in other directions. Thus, the second principle of ordering
neighbor lists is that the top members of each list should contain neighbors
from all directions.

5\. Solver. Use GFD to compute spatial derivatives. The optimal number of
neighbors for the GFD solver is obtained from numerical experiments. 8 and 54
neighboring particles are used in upwind and Lax-Wendroff schemes to compute
spatial derivatives. Then we compute time derivatives based on the numerical
discretization.

6\. Particle mover. Update states of particles and move them to the next
location in space.

The major octree algorithms such as building, refining and searching that
significantly affect the code performance are parallelized using the p4est
library. p4est [14] enables a dynamic management of a collection of adaptive
octrees on distributed memory supercomputers. It has the functionality of
building, refining, coarsening, 2:1 balancing and partitioning on
computational domains composed of multiple connected two-dimensional quadtrees
or three-dimensional octrees, referred to as a forest of octrees. Here we
briefly introduce three features of p4est which greatly contribute to the
massively parallel LP code.

1\. Space-filling curve. Octants within an octree can be assigned a natural
ordering by traversing across all leaves. This results in a one-dimensional
sequence corresponding to a space-filling Z-shaped curve in the geometric
domain. p4est extends this concept to a forest of octrees by connecting the
space-filling curve between octrees, thus generating a total ordering of all
octants in the domain (see Figure 1).

2\. Parallel partition. The objective of the partition algorithm is to
equidistribute the computational work uniformly among processes, which is
necessary to ensure parallel scalability of applications. In addition to
uniformly distributing octants, p4est allows for a user-defined weight
function that returns a non-negative integer weight for each octant and
creates a partition that is evenly distributed by weight. After the partition,
the relative order of octants with respect to processes stays the same as the
order in the space-filling curve.

Figure 1: A 2D example with two octrees($k_{0}$ and $k_{1}$). A left-to-right
traversal of leaves through two octrees creates a space-filling z-curve(black
line). In this example the forest is partitioned among three
processes(represented by red, green and blue) and this partition divides the
geometric domain into three segments with almost equal ($\pm 1$) number of
octants (see [14]).

3\. Ghost. p4est is parallelized using the Message Passing Interface (MPI),
which enables its applications to run on distributed memory supercomputers.
For many computational applications, it is useful to know information on
remote octants. Ghost algorithm collects one layer of neighboring octants on
remote processes, providing the complete parallel neighborhood information.
This feature reduces the communication pattern in the application to simple
point-to-point transfers of numerical information.

The massively parallel LP code stores the particle data in a single array.
Elements in the array contain all information about particles, such as
coordinates, velocity, hydrodynamic states and supporting quantities. Not only
does this make the code easier to maintain and extend, but also it simplifies
the work to communicate particle data among processes as the particle data is
stored on a contiguous memory block. The particle data is uniformly
distributed among available processes. Particles on local processes are sorted
by following the order of their owner octants in the space-filling curve.
Particles in the same octants are not imposed with any specific orders.

The LP code has been parallelized using MPI. The macro-structure of the
octrees and some global information, such as boundary conditions and initial
states, are shared among processes. The micro-structure of octrees, i.e., the
division of each octree into octants, particle data and some other highly
memory-demanding data structures designed for certain problems are dynamically
managed and distributed among all processes. In order to achieve balanced
computational loads across processes, we assign weights in the p4est Partition
algorithm as the number of computational particles in octants. As a result,
local numbers of octants vary among processes if the particle distribution is
non-uniform, but computational particles are evenly distributed among
processes. This partition optimizes the load balance as the neighbor searching
and stencil operations dominate the computational time.

p4est allows computational domains to be composed of multiple connected
octrees, referred to a forest of octrees. Neighboring octrees can be
arbitrarily rotated against each other, communicating through non-local
connections. This flexible encoding of octrees enables implementations of a
wide range of possible computational domains and periodic boundary conditions
without any special treatment to stencils.

The LP code strictly follows the idea of objected-oriented design principles,
dividing the algorithm into major tasks and implementing them with classes. LP
supports algorithm extensions using polymorphic classes and public interface.
It consists of several main classes. Initializer reads the input information
from a formatted input file and saves data into class member variables. These
variables are used as inputs to member functions of several other classes to
generate initial particle distributions, particle states, EOS models, boundary
conditions, octree structures, data analysis etc. Global_data stores the
particle data, EOS model, boundary model, octree data structures and
implements particle initialization routines, particle migration, particle
neighbor searching and ghost particle generation algorithms. LP_solver
implements the first order upwind scheme and the second order Lax-Wendroff
scheme for 2D and 3D problems. Octree_manager implements wrapper functions for
building, adapting, balancing, and partitioning octrees using p4est library
functions. Particle_viewer outputs particle data stored on local processes in
the VTK format. p4est-based implementation of another Lagrangian particle
based code AP-Cloud (Adaptive Particle-in-Cloud) for optimal solutions of
Vlasov-Poisson problems can be found in [15].

###  2.3 Verification and parallel scalability

The Lagrangian Particle code has been extensively verified and validated using
both classical hydrodynamic problems [9] and applied problems of plasma
physics relevant to thermonuclear fusion [11, 16]. Since the main novelty of
the algorithmic part of this paper is the parallelization of the Lagrangian
particle code for distributed memory supercomputers, it is sufficient to show
that the parallel code gives the same accuracy as the original code for some
typical hydrodynamic problems. We use the classical Yee vortex problem to
demonstrate the numerical accuracy and convergence.

The Yee vortex [17] is an isentropic 2D vortex in a free-stream flow with the
velocity $(u_{\infty},v_{\infty})=(1,0)$, density $\rho_{\infty}=1$, and
pressure $p_{\infty}=1$. The initial conditions for the Yee vortex are as
follows

$\begin{split}&\rho=T^{\frac{1}{\gamma-1}}=[1-\frac{(\gamma-1)\beta^{2}}{8\gamma\pi^{2}}e^{1-r^{2}}]^{\frac{1}{\gamma-1}},\\\
&u=1-\frac{\beta}{2\pi}e\frac{1-r^{2}}{2}(y-y_{0}),\\\
&v=\frac{\beta}{2\pi}e^{\frac{1-r^{2}}{2}}(x-x_{0}),\\\
&P=\rho^{\gamma},\end{split}$ (20)

where $(x_{0},y_{0})$ are coordinates of the vortex center,
$r^{2}=(x-x_{0})(x-x_{0})+(y-y_{0})(y-y_{0})$, and $\beta$ is the vortex
strength. Notice that the vortex extends to the entire 2D space, but it can be
truncated for all practical purposes as the velocity rapidly decays to the
free-stream velocity. The exact solution with given initial states is just a
passive convection of the vortex moving with the free-stream velocity.
Therefore, it provides a good accuracy measure for numerical schemes for the
nonlinear Euler equations. For the numerical test, we select $\beta=5$ and
$\gamma=1.4$. Simulations that used a perfect hexagonal initial placement of
particles and small random perturbations of initial particle coordinates
produced very similar results in terms of numerical errors. The norms of
numerical errors of pressure and velocity for the initial hexagonal placement
of particles are given in Tables 1 and 2.The convergence order is higher than
second in both $L_{2}$ and $L_{\infty}$ norms. Figure 2 shows radial scatter
plots for different numerical resolutions and visually demonstrates the
numerical convergence. While the problem is two-dimensional, it has also been
performed with the full-3D code by assuming constant states in the third
dimension, thus testing the performance of data structure and functions
relevant to octrees. 3D results were used for the parallel scalability
analysis below.

dx | $L_{2}$ | $L_{\infty}$ | ratio | order
---|---|---|---|---
0.4 | 0.02508 | 0.1488 | - | -
0.2 | 0.00549 | 0.0366 | 4.57 | 2.19
0.1 | 0.0011 | 0.0063 | 4.99 | 2.32
0.05 | 0.00029 | 0.0017 | 3.77 | 1.91
Table 1: Errors of pressure at $t=8$ and convergence rates for 2D Yee vortex problem. dx | $L_{2}$ | $L_{\infty}$ | ratio | order
---|---|---|---|---
0.4 | 0.067 | 0.369 | - | -
0.2 | 0.01789 | 0.1035 | 3.75 | 1.91
0.1 | 0.00324 | 0.00241 | 5.52 | 2.47
0.05 | 0.00056 | 0.00047 | 5.78 | 2.53
Table 2: Errors of velocity at $t=8$ and convergence rates for 2D Yee vortex
problem.

(a) Pressure

(b) Velocity

Figure 2: Radial scatter plots of particles states of Yee vortex simulation at
$t=8$

We have studied parallel scalability of the code based on the Yee vortex
problem in 3D using the Seawulf cluster located at Stony Brook University. We
tested the strong scaling by running the Yee vortex simulation with different
number of CPU cores while keeping the total number of computational particles
constant (1,060,000 particles were used in the strong scalability tests).
Figure 3 shows the strong scalability of the LP code.

Figure 3: Strong scaling for massively-parallel LP code with using 1,060,000
particles. The largest speedup (22.28) was reaches at 240 CPU cores, which was
55.7% of ideal speedup.

The weak scaling was measured by running the code with different number of CPU
cores and using a proportionally scaled number of particles such that the
number of particle per core remained constant. The results of weak scaling
tests are demonstrated in Figure 4.

Figure 4: Weak scaling for massively-parallel LP code. The largest speedup
(22.89) was reached at 192 CPU cores, which was 71.5% of ideal speedup.

The tests demonstrated satisfactory scalability. As expected, the weak scaling
outperformed the strong one. Since only the weak scaling is important for
scientific computing in our area (which requires solving problems of ever-
increasing size in a reasonable amount of time), achieving more than 70% of
ideal speedup and a linear scaling is sufficient for practical applications.
It is important to note that the saturation of scaling has not been reached as
we limited the test to 200 CPU cores, a typical allocation for our
applications runs. We have reasons to expect that the weak scaling beyond 200
cores will remain close to linear.

##  3 Applications

###  3.1  Summary of thermonuclear fusion applications

The Lagrangian particle code has been broadly applied to the simulation of
impurities in hot thermonuclear fusion plasmas. Thermonuclear fusion devices
called tokamaks explore the magnetic confinement of hot deuterium-tritium
plasma for achieving conditions of thermonuclear fusion burn. The
International Thermonuclear Experimental Reactor (ITER) [www.iter.org] is an
international nuclear fusion research and engineering project that builds a
new generation tokamak in Cadarache, France. Upon completion, the ITER tokamak
is expected to achieve the fusion energy gain of at least 10 times larger
compared to the energy deposited into the plasma. The ITER tokamak requires
the injection of cryogenic material into the plasma for the purpose of fueling
and plasma disruption mitigation.

As a cryogenic pellet is injected into a thermonuclear plasma, it is hit by
hot plasma electrons that rapidly ablate the pellet surface leading to the
formation of a cold and dense ablation cloud. This ablated cloud absorbs most
of the incoming plasma electrons. The fraction of electrons that reach the
pellet surface is responsible for the pellet ablation rate. In our recent
works, 2D axisymmetric [18] and full 3D simulations [11, 16] of ablating
pellets in tokamaks have been performed, and the influence of tokamak plasma
conditions, magnetic field and 3D plasma drift effects on their ablation rates
have been studied and compared with the available experimental data. The
pellet ablation problem is intrinsically multiscale. The resolution of large
gradients of plasma states occurring at sub-millimeter length scales near the
pellet surface is critically important for computing pellet ablation rates. At
the same time, large, device-scale expansion of the ablated material along
tokamaks change the background plasma properties and affect the distribution
of plasma heat sources for the pellet ablation. The adaptive property of the
Lagrangian particle method made it possible to accurately resolve small scale
processes near the pellet surface while propagating the ablated material over
distances of dozens of meters.

This section intended to emphasize adaptive multi-physics capabilities of the
Lagrangian particle code. Detailes of the plasma physics simulations
summarized above go beyond the scope of this paper and can be found in the
referenced papers.

###  3.2 Supersonic jets for LPWA

In this section, we present applications of the Lagrangian Particle code to
the problem of nozzle design optimization for laser plasma wakefield
acceleration (LPWA) experiments being conducted at Brookhaven National
Laboratory. LPWA is potentially a high-impact research area. By replacing
conventional, large-size accelerator structures with electromagnetic wakes
created in plasma by lasers or high energy charged particles, LPWA is expected
to increase the accelerating field gradients by orders of magnitude, leading
to the reduction of accelerator sizes from the scale of tens of kilometers to
the meter scale.

Broad studies of the plasma wakefield acceleration of electrons are being
conducted at the Brookhaven National Laboratory (BNL) based on the powerful
CO2 laser facility. Our research group supports this theoretical and
experimental effort by numerical simulations. The main processes of interest,
the CO2 laser interaction with hydrogen gas, ionization of hydrogen, formation
of plasma wakes and the subsequent acceleration of electrons are described by
electromagnetic system of equations for relativistic particles and plasma. The
main simulation code is SPACE [15], a parallel, 3D electromagnetic Particle-
in-Cell code with atomic physics support, developed in our research group.
Recent simulations results are described in [19, 20].

There is one component of the laser plasma wakefield experiments that is in
the range of classical hydrodynamics. The problem is to obtain a proper
distribution of hydrogen gas in vacuum before its interaction with the laser.
Various designs have been considered for this tasks, including a hydrogen
chamber with special windows for the laser beam. While various pro- and
contra- of such chambers are being discussed, the BNL experiment currently
operates with a simple solution - a hydrogen jet injected into the interaction
region through a nozzle connected to a high-pressure reservoir. While being
the simplest solution from the engineering point of view, the nozzle jet does
not provide completely uniform initial gas distribution, which is desirable
for the experiment. The purpose of numerical simulations is to optimize the
process and obtain a nozzle that provides the most uniform hydrogen density
distribution of specific requirements. A schematic of the BNL laser plasma
wakefield experiment is shown Figure 5.

Figure 5: Schematic of experiment on laser-plasma wakefield acceleration. CO2
and YAG denote laser beams generated by the CO2 and solid state lasers,
respectively, propagating across the hydrogen jet.

A nozzle is typically attached to a high-pressure reservoir containing the
working gas. Computationally, we model the presence of the reservoir by
special boundary conditions at the nozzle inlet. Three hypotheses are used in
order to calculate the pressure, density and velocity of the particles
entering the nozzle.

  * 1.

The pressure of newly injected particles into the nozzle is equal to the
average of the pressure of the particles in a small neighborhood of the nozzle
inlet for the current timestep. (For the first timestep the pressure of the
reservoir is used.)

  * 2.

The entropy in the reservoir is equal to the entropy in the neighborhood of
the nozzle inlet. In terms of pressure and density, this condition can be
written as:

$\frac{P_{res}}{\rho_{res}^{\gamma}}=\frac{P_{inf}}{\rho_{inf}^{\gamma}},$

where the subscripts ’res’ and ’inf’ denote the reservoir and nozzle inflow
quantities.

  * 3.

The transverse velocity component of particles is negligibly small (assumed
zero), and the longitudinal velocity component is calculated using the
enthalpy condition as

$\frac{1}{2}u_{res}^{2}+\frac{1}{\gamma-1}c^{2}_{res}=\frac{1}{2}u_{inf}^{2}+\frac{1}{\gamma-1}c^{2}_{inf}$

where $c$ is the sound speed, given by $c^{2}=\gamma\frac{P}{\rho}$. Assuming
$u_{res}=0$, the velocity of the inflow particles can be calculated.

Two types of nozzles have been used in LPWA experiments, a cylindrical-
divergent nozzle and a convergent-divergent one. Their schematics and
geometric measurements are given in Figure 6, and typical 3D numerical
simulations are shown in Figures 7(a) and 7(b).

(a) Convergent-divergent nozzle

(b) Cylindrical-divergent nozzle

Figure 6: Schematics and design parameters for convergent-divergent nozzle (a)
and cylindrical-divergent nozzle (b)

(a) Jet from convergent-divergent nozzle

(b) Jet from cylindrical-divergent nozzle

Figure 7: LP simulation of hydrogen jets. Pressure distributions at $t=50\mu
s$ are shown.

The radial profiles of hydrogen density for the jets formed by the
cylindrical-divergent nozzle at various reservoir parameters, used for
determining plasma density after the interaction with laser, are shown in
Figure 8. Similar profiles were also obtained for the convergent-divergent
nozzle.

Figure 8: Radial profiles of hydrogen density at 5 mm form the nozzle exit for
the reservoir temperature 293K and the following values of the reservoir
pressure (bar): 1 (A), 1.27 (B), 1.72 (C), and 4.21 (D).

By optimizing the convergent-divergent nozzle shape, we were able to obtain
reasonably small variations of hydrogen density over the distance of 6 mm as
shown in Figure 9. The optimized nozzle dimensions are as follows: diameter of
the reservoir opening is 1 mm, length of the convergent section is 4 mm,
diameter of the neck is 0.5 mm, length of the divergent section is 8 mm, and
diameter of the nozzle exit is 4.1 mm.

Figure 9: Hydrogen density for the optimized nozzle at several distances from
the nozzle exit.

##  4 Summary and conclusions

In this paper, we described a parallel Lagrangian particle algorithm for
compressible hydrodynamic problems and the corresponding object-oriented code
for massively-parallel distributed memory supercomputers. The Lagrangian
particle (LP) algorithm implements an unsplit Lax – Wendroff type scheme with
a limiter. In the location of points selected by a state estimator, the
Lax–Wendroff scheme is blended with a first-order upwind method to control
numerical dispersion. In the LP method, approximations of spatial derivatives
are obtained by employing a local polynomial fit known also as the generalized
finite difference (GFD) method. The LP method is second-order accurate in
space and time and it is extendable to higher order accuracy.

The particle manager module of the massively-parallel LP code is based on the
p4est library. p4est enables a dynamic management of a collection of adaptive
octrees on distributed memory supercomputers. It has the functionality of
building, refining, coarsening, 2:1 balancing and partitioning of
computational domains composed of multiple connected two-dimensional quadtrees
or three-dimensional octrees, referred to as a forest of octrees. The
Lagrangian particle code has been verified, validated, and tested on
distributed memory supercomputers. It demonstrated good strong and weak
scalability on hundreds of supercomputer cores.

The massively-parallel Lagrangian particle code is applicable to a variety of
fundamental science and applied problems. Its current most important
application is simulation of impurities in thermonuclear fusion devices. We
briefly summarize the thermonuclear fusion applications in this paper and
referred the reader to published works for further details. Hydrodynamic
simulations for the nozzle design optimization for supersonic hydrogen jets
are described in more detail. Such simulations support experimental work at
Brookhaven National Laboratory on CO2 laser-driven plasma wakefield
acceleration.

Acknowledgments

This work was supported in part by the U.S. Department of Energy grants SciDAC
Center for Tokamak Transient Simulations and Grant No. DE-SC0014043.

## References

  * [1] J. Monaghan, Smoothed particle hydrodynamics, Rep. Prog. Phys. 68 (2005) 1703–1759.
  * [2] N. Frontiere, C. Raskin, J. Owen, CRKSPH conservative reproducing kernel smoothed particle hydrodynamics scheme, J. Comput. Phys. 332 (2017) 160–209.
  * [3] D. Avesani, M. Dumbser, A. Bellin, A new class of moving-least-squares WENO-SPH schemes, J. Comput. Phys. 270 (2014) 278–299.
  * [4] J. Monaghan, Smoothed particle hydrodynamics, Annual Review of Astronomy and Astrophysics 30 (1992) 543–574.
  * [5] G. Diltz, Moving-least-squares particle hydrodynamics - I. consistency and stability, Int. J. Num. Methods in Engineering 44 (1999) 1115–1155.
  * [6] P. Hopkins, GIZMO: a new class of accurate, mesh-free hydrodynamic simulation methods, Mon. Not, R. Astron. Soc. 450 (2014).
  * [7] W. K. Liu, J. Sukky, Y. F. Zhang, Reproducing kernel particle methods, International Journal for Numerical Methods in Fluids 20 (1995) 1081–1106.
  * [8] G. Oger, An improved SPH method: Towards higher order convergence, J. Comput. Phys. 225 (2007) 1472–1492.
  * [9] R. Samulyak, X. Wang, H.-C. Chen, Lagrangian particle method for compressible fluid dynamics, J. Comput. Phys. 362 (2018) 1–19.
  * [10] J. Benito, F. Urena, L. Gavete, Influence of several factors in the generalized finite difference method, Appl. Math. Modelling 25 (12) (2001) 1039–1053.
  * [11] R. Samulyak, S. Yuan, N. Naitlho, P. Parks, Lagrangian particle model for 3d simulation of pellets and SPI fragments in tokamaks, Nuclear Fusion 61 (4) (2021) 046007.
  * [12] M. Zepeda, R. Samulyak, Computational model for granular flow based on the lagrangian particle method, Computer Methods in Applied Mechanics and EngineeringSubmitted (2022).
  * [13] D. Meagher, Geometric modeling using octree encoding, Computer Graphics and Image Processing 19 (1982) 129–147.
  * [14] C. Burstedde, L. C. Wilcox, O. Ghattas, P4EST: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees, SIAM J. Sci. Comp. 33 (3) (2011) 1103–1133.
  * [15] K. Yu, P. Kumar, S. Yuan, A. Cheng, R. Samulyak, SPACE: 3d parallel solvers for vlasov-maxwell and vlasov-poisson equations for relativistic plasmas with atomic transformations, Comp. Physics Communications 277 (2022) 108396.
  * [16] E. Hollmann, N. Naitlho, S. Yuan, R. Samulyak, P. Parks, D. Shiraki, J. Herfindal, C. Marini, Measurement and simulation of small cryogenic neon pellet Ne-I 640 nm photon efficiency during ablation in DIII-D plasma, Phys. PlasmasSubmitted (2022).
  * [17] H. C. Yee, M. Vinokur, M. J. Djomehri, Entropy splitting and numerical dissipation, Journal of Computational Physics 162 (1) (2000) 33–81.
  * [18] N. Bosviel, P. Parks, R. Samulyak, Near-field models and simulations of pellet ablation in tokamaks, Physics of Plasmas 28 (1) (2021) 012506.  doi:10.1063/5.0029721.
  * [19] P. Kumar, K. Yu, R. Zgadzaj, L. Amorim, M. Downer, J. Welsh, V. Litvinenko, N. Vafaei-Najafabadi, R. Samulyak, Simulation study of co2 laser-plasma interactions and self-modulated wakefield acceleration, Physics of Plasmas 26 (8) (2019) 083106.
  * [20] P. Kumar, K. Yu, R. Zgadzaj, M. Downer, I. Petrushina, R. Samulyak, V. Litvinenko, N. Vafaei-Najafabadi, Evolution of the self-injection process in long wavelength infrared laser driven lwfa, Physics of Plasmas 28 (1) (2021) 013102.