Present address: ] Physics Division, Argonne National Laboratory, Argonne, IL
60439, USA.

# Elastic scattering measurements for the 10C + 208Pb system at
E${}_{\textnormal{lab}}$ = 66 MeV

R. Linares<EMAIL_ADDRESS>Instituto de Física, Universidade Federal
Fluminense, 24210-340, Niterói, Rio de Janeiro, Brazil  Mandira Sinha
Department of Physics, Bose Institute, 93/1 Acharya Prafulla Chandra Road,
Kolkata-700009, W.B, India  E. N. Cardozo  Instituto de Física, Universidade
Federal Fluminense, 24210-340, Niterói, Rio de Janeiro, Brazil  Instituto de
Física, Universidade de São Paulo, 05508-090 São Paulo, Brazil  V. Guimarães
Instituto de Física, Universidade de São Paulo, 05508-090 São Paulo, Brazil
G. V. Rogachev  Department of Physics and Astronomy, Texas A & M University,
77843 TX  Cyclotron Institute, Texas A & M University, 77843 TX  Nuclear
Solutions Institute, Texas A & M University, 77843 TX  J. Hooker  Department
of Physics and Astronomy, Texas A & M University, 77843 TX  Cyclotron
Institute, Texas A & M University, 77843 TX  E. Koshchiy  T. Ahn  Cyclotron
Institute, Texas A & M University, 77843 TX  C. Hunt  Department of Physics
and Astronomy, Texas A & M University, 77843 TX  Cyclotron Institute, Texas A
& M University, 77843 TX  H. Jayatissa  [  Department of Physics and
Astronomy, Texas A & M University, 77843 TX  Cyclotron Institute, Texas A & M
University, 77843 TX  S. Upadhyayula  Department of Physics and Astronomy,
Texas A & M University, 77843 TX  Cyclotron Institute, Texas A & M University,
77843 TX  B. Roeder  A. Saastomoinen  Cyclotron Institute, Texas A & M
University, 77843 TX  J. Lubian  Instituto de Física, Universidade Federal
Fluminense, 24210-340, Niterói, Rio de Janeiro, Brazil  M. Rodríguez-Gallardo
Departamento de Física Atómica, Molecular y Nuclear, Facultad de Física,
Universidad de Sevilla, Apartado 1065, E-41080 Sevilla, Spain
and Instituto Carlos I de Física Teórica y Computacional, Universidad de
Sevilla, Spain  J. Casal  Dipartimento di Fisica e Astronomia “G. Galilei”,
Università degli studi di Padova, and INFN - Sezione di Padova, Via Marzolo 8,
I-35131 Padova, Italy  K. C. C. Pires  Instituto de Física, Universidade de
São Paulo, 05508-090 São Paulo, Brazil  M. Assunção  Departamento Física,
Universidade Federal de São Paulo, Campus Diadema, 09913-030, Diadema, São
Paulo, Brazil  Y. Penionzhkevich  S. Lukyanov  Joint Institute for Nuclear
Research, Dubna. Moscow region, 141980

###### Abstract

Background: The influence of halo structure of 6He, 8B, 11Be and 11Li nuclei
in several mechanisms such as direct reactions and fusion is already
established, although not completely understood. The influence of the 10C
Brunnian structure is less known.

Purpose: To investigate the influence of the cluster configuration of 10C on
the elastic scattering at an energy close to the Coulomb barrier.

Methods: We present experimental data for the elastic scattering of the
10C+208Pb system at $E_{\rm lab}$ = 66 MeV. The data are compared to the
three- and the four-body continuum-discretized coupled-channels calculations
assuming 9B+$p$, 6Be+$\alpha$ and 8Be+$p$+$p$ configurations.

Results: The experimental angular distribution of the cross sections shows the
suppression of the Fresnel peak that is reasonably well reproduced by the
continuum-discretized coupled-channels calculations. However, the calculations
underestimate the cross sections at backward angles. Couplings to continuum
states represent a small effect.

Conclusions: The cluster configurations of 10C assumed in the present work are
able to describe some of the features of the data. In order to explain the
data at backward angles, experimental data for the breakup and an extension of
theoretical formalism towards a four-body cluster seem to be in need to
reproduce the measured angular distribution.

##  I Introduction

The investigation of nuclear structure and reaction mechanisms induced by
radioactive nuclei is a subject of considerable interest. Detailed reviews on
this subject are presented in refs. Kolata, J. J.  _et al._ (2016); Canto  _et
al._ (2015a); Canto, L. F.  _et al._ (2020). Some of the weakly-bound nuclei
are called exotic as they exhibit low binding energies and halo structure, in
which the valence nucleon(s) orbits a compact inert core and therefore forms
an extended matter distribution Tanihata  _et al._ (1985, 2013). These
features of the exotic nuclei manifest themselves in nuclear collisions mainly
at energies around the Coulomb barrier. At these energies, a strong interplay
between elastic scattering and other processes, such as breakup and transfer
reactions, is expected. For instance, the low binding energy produces a
decoupling between the valence nucleon(s) and the core during the collision
which might cause suppression of the Fresnel peak (Coulomb rainbow) in the
angular distributions Kolata, J. J.  _et al._ (2016).

Also, the decoupling to the continuum of the single-particle or cluster motion
of the valence nucleon(s) with respect to the core has motivated the
application of reaction models such as continuum-discretized coupled-channels
(CDCC) Austern  _et al._ (1987), in which the coupling between the elastic and
the breakup channel is explicitly taken into account. The CDCC calculations
have been applied to study the elastic scattering of exotic neutron-rich
nuclei, such as 6He, 11Li and 11Be, on several medium to heavy target nuclei
and at energies around the Coulomb barrier, Morcelle  _et al._ (2014); Di
Pietro  _et al._ (2004); de Faria  _et al._ (2010); Sánchez-Benítez  _et al._
(2008); Acosta  _et al._ (2011); Aguilera  _et al._ (2011); Di Pietro  _et
al._ (2010, 2012); Acosta  _et al._ (2009); Pesudo  _et al._ (2017); Mazzocco
_et al._ (2006). In some of these works, the 11Be has been described as 10Be
core and a valence neutron while 6He has been described as a 4He core and a
paired di-neutron system. This is referred as three-body CDCC (3b-CDCC)
framework. An extension of this framework allows to treat the projectile as a
core plust two valence particles configuration, simply referred as the four-
body CDCC (4b-CDCC) calculations Rodríguez-Gallardo  _et al._ (2009). The
4b-CDCC formalism has been applied to study reactions with 6He Morcelle  _et
al._ (2014), 9Be Casal  _et al._ (2015); Arazi  _et al._ (2018) and 11Li
Cubero  _et al._ (2012); Fernández-García  _et al._ (2013), in which the 6He
and 11Li are treated as a core+$n$+$n$ and 9Be as a $\alpha$+$\alpha$+$n$
configurations. These cluster configurations are quite important to describe
the elastic scattering data.

Moving towards the proton-rich side, one of the most studied nucleus is 8B.
This nucleus has a very weakly bound valence proton with separation energy
$S_{p}=0.138$ MeV. The centrifugal barrier between the valence proton, in the
$p$-orbital, and the core inhibits the wider radial extension. Therefore the
halo structure is relatively less pronounced when compared to the neutron-halo
nuclei. Elastic scattering for the 8B+58Ni system was investigated in Refs.
Aguilera  _et al._ (2009); Lubian  _et al._ (2009); Keeley  _et al._ (2010).
The observed large reaction cross sections has been interpreted as an
indication of halo configuration for this nucleus. This is also important to
describe the 8B $\rightarrow$ 7Be + $p$ breakup cross sections at sub-Coulomb
energy Guimaraes  _et al._ (2000).

The effect of the breakup channel on the elastic scattering can be studied in
collisions with heavier targets due to the increasing predominance of the
long-range Coulomb interaction over the nuclear interaction, that induces a
strong dipole response. In Ref. Yang  _et al._ (2013) it has been reported
measurements of the elastic scattering in the 8B + 208Pb system at
E${}_{\textnormal{lab}}$ = 170 MeV, that corresponds to about three times the
Coulomb barrier height. The effect of breakup on the elastic scattering is
small at this bombarding energy but becomes relevant at energies around the
Coulomb barrier, as reported in Refs. Rangel  _et al._ (2016); Yang  _et al._
(2016). Recently, Mazzocco et al., reported a measurement of the elastic
scattering for the same 8B + 208Pb system at E${}_{\textnormal{lab}}$ = 50 MeV
Mazzocco  _et al._ (2019). The data indicate a strong suppression of the
Fresnel peak and a large reaction cross section. However, the 3b-CDCC
calculations do not provide a full description of the data.

The 10C is another interesting proton-rich nucleus due to its Brunnian (super-
Borromean) structure where the four interconnected rings are associated to the
interactions between $\alpha$+$\alpha$+$p$+$p$ Curtis  _et al._ (2008). The
removal of one particle breaks the remaining system apart. This nucleus can
decay by three possible channels: 8Be + $p$ \+ $p$ , 9B + $p$ and 6Be+$\alpha$
with binding energies of 3.821, 4.006 and 5.101 MeV, respectively. Since 8Be,
9B and 6Be are unbound systems, the final state is the $\alpha+\alpha+p+p$
system. Breakup measurement of the 10C on 12C at 33 MeV/A indicates that the
final state is predominantly fed by the 9B+$p$ and 8Be+$p$+$p$ channels Curtis
_et al._ (2008).

Measurementes of the elastic scattering has been recently performed for the
10C+58Ni at E${}_{\textnormal{lab}}$ = 35.3 MeV Guimarães  _et al._ (2019).
Couplings to the first excited state in 10C ($2_{1}^{+}$) has been quite
important in the description of cross sections at backward angles, indicating
the importance of the large deformation for this nucleus. Analysis within the
3b-CDCC, assuming the 9B + $p$ and 6Be+$\alpha$ configurations, also indicates
a strong cluster configuration in the 10C. Elastic measurements of 10C on
208Pb target were performed at 226 and 256 MeV Yang  _et al._ (2014), well
above the Coulomb barrier, and the angular distributions exhibit a typical
Fresnel-like shape.

In this work we present new experimental data for the elastic scattering of
10C on 208Pb target measured at E${}_{\textnormal{lab}}=66$ MeV, which
corresponds to an energy close to the Coulomb barrier (V${}_{\textnormal{B}}$
= 61.3 MeV in the laboratory framework). The theoretical analyses are carried
out within the optical model, to determine the reaction cross section, and the
coupled channel, 3b- and 4b-CDCC formalisms to assess the contribution of
individual reaction channels. The paper is organized as follows: the
experimental details and the theoretical analysis are discussed in sections II
and III, respectively. The conclusions are given in section IV.

##  II Experimental details

The experiment was performed at the Cyclotron Institute of Texas A&M
University. The radioactive 10C beam was produced by the charge-exchange
p(10B,10C)n reaction. The primary 9.6 MeV/u 10B beam was delivered by the K500
superconducting cyclotron and impinged on a H2 gas cell cooled to 77 K. The
produced 10C ions were selected by the Momentum Achromatic Recoil Spectrometer
(MARS) and directed to the scattering chamber Tribble  _et al._ (2002). The
production rate of the radioactive beam was monitored as the ratio between the
thin plastic scintillators counting, at the entrance of the scattering
chamber, and the primary beam current collect at the gas cell. The average
production rate was $5\times 10^{3}$ pps for a 300 nA incident 10B+5 beam. The
obtained purity for 10C beam was better than 95%.

The schematic diagram of the experimental arrangement in the scattering
chamber is shown in Fig. 1. The target consisted of two isotopically enriched
208Pb foils, 2.37 mg/cm2 thick, evaporated on a 0.08 mg/cm2 thick carbon
backing. These targets were mounted together on the target frame. The
detection system consisted of a double-sided silicon strip detectors (DSSSD)
$\Delta$E-E1 telescope placed at forward scattering angles and a single DSSSD
(E2) at backward scattering angles. The E1 and E2 detectors have an active
area of 97.2 mm $\times$ 97.2 mm, 500 $\mu$m thick and with front and back
sides segmented in 128 strips. Both detectors were placed 100 mm away from the
target position. The $\Delta$E detector of the telescope has an active area of
49.5 mm $\times$ 49.5 mm, 38 $\mu$m thick, and had front and back sides
segmented in 16 strips. This detector was placed 50 mm away from the target
position. These detectors were mounted in a diamond-like geometry to allow for
a wide angular coverage. In this configuration, the telescope covered the
scattering angles between $\theta_{\rm lab}=15\degree$ and $89\degree$ whereas
the E2 covered the angular range from $\theta_{\rm lab}=87\degree$ to
$161\degree$.

An extra $16\times 16$ DSSSD detector was placed at $0\degree$, downstream the
beam and at the back of the chamber to monitor the radioactive 10C beam. To
avoid damages in this monitor due to direct exposure to the secondary beam, a
screen with four 2-millimeter holes, positioned at the vertices of a 20 mm
square, was mounted in front of the detector. The relative intensities
registered at these positions were used to determine the central position of
the secondary beam as well as the 10C beam spot size, estimated to be 6 mm.

Figure 1: (Color online) Sketch of the detection system for the measurements,
composed by one $\Delta E-E$ telescope, placed at forward angles, and a
E-detector at backward angles. The E-detectors are $128\times 128$ DSSSD and
the $\Delta E$ is a $16\times 16$ thin DSSSD. A third DSSSD (monitor) was
placed at $0\degree$ to monitor the beam axis.

The two-dimensional particle identification spectrum, $\mathit{\Delta
E-E\textsubscript{\rm res.}}$, obtained with the telescope is shown in Fig. 2.
The elastically scattered 10C ions are well separated from the 10B ions. The
small amount of 10B ions went through the MARS beam line due to sparks in the
Wien filter during the acquisition runs. The observed light particles
($\alpha$ and protons) are from fusion-evaporation and break-up reactions.
However, according to calculations of the fusion-evaporation process performed
with the the pace code Tarasov and Bazin (2008); Gavron (1980), there is no
significant contribution of p and $\alpha$ particles in the energy region of
interest of the elastic scattering at backward angles (detector E2).

Figure 2: (Color online) The $\Delta E-E$ plot for particle identification. A
small fraction of 10B is observed due to sparks in the Wien filter of MARS
separator during the acquisition runs.

The total energy spectrum (E${}_{\textnormal{total}}$ = $\Delta E+E$) at 3
different scattering angles are shown in Fig. 3. The achieved energy
resolution is about 2.9 MeV (Full Width at Half Maximum). This energy
resolution does not allow a complete separation of the elastic peak from the
inelastic scattering leading to the $3^{-}_{1}$ state in 208Pb (at 2.61 MeV).
At the very backward angles ($\theta_{\textnormal{lab.}}>120\degree$), yields
of the elastic contribution were determined considering 2-Gaussian curve fits,
being the second Gaussian associate to the inelastic scattering contributions
of the $2^{+}$ state in 10C and $3^{-}$ in 208Pb Therefore, in principle, this
is a quasi-elastic scattering measurement. Also, according to calculations,
the contribution of inelastic scattering is found to be quite small, as it
will be shown in the next sections. Also, the spectra at backward angles show
almost no background in the elastic peak region, which corroborate with the
absence of evaporated $\alpha$ in the region of interest. Considering the
energy losses in the scintillators ($14\mu$m thick), upstream the scattering
chamber, and in the Pb+C target, the averaged energy is 65.7 MeV. In
particular, the estimated energy loss in the target is 2.1 MeV.

Figure 3: (Color online) Example of total energy spectrum for 10C particles
scattered by 208Pb at (a) $30\degree\leq\theta_{\textnormal{lab}}\leq
35\degree$, (b) $105\degree\leq\theta_{\textnormal{lab}}\leq 115\degree$ and
(c) $135\degree\leq\theta_{\textnormal{lab}}\leq 140\degree$. In panels (a)
and (b), the elastic peak is represented by a single gaussian curve. In panel
(c), data points were adjusted to a two-gaussian curve with the same standard
deviation. The green dotted and violet dashed-dot curves correspond to elastic
and the contribution from $2^{+}$ state in 10C and $3^{-}$ in 208Pb,
respectively.

The yields of 10C scattered particles were integrated in angular ranges of
$5\degree$ or $10\degree$, depending on the counting statistics. The
scattering angle of each detected particle was determined from the coordinates
of the strips triggered in the E-detector. The angular resolution is estimated
to be about $3\degree$ (in the laboratory framework), where the beam spot size
at the target, the relative distances between target holder and E-detectors
and the strip width have been taken into account. The solid angles were
determined using an alpha source placed at the target position and assuming an
uniform emission. The final overall normalization was performed by considering
the cross sections for $\theta_{\rm lab}<45\degree$ as Rutherford cross
section.

The experimental angular distribution of the elastic cross sections is shown
in Fig. 4. The error bars in the experimental cross sections correspond to the
statistic uncertainties only. Systematic uncertainties due to solid angles and
the beam spot size at the target introduce about 10% uncertainty to the cross
sections. The angular distribution for the elastic scattering of 12C + 208Pb
system at Elab = 64.9 MeV Santra  _et al._ (2001) is also plotted in the
figure (blue triangle points). The comparison of both angular distributions
highlights the suppression of the Fresnel peak for the 10C projectile. The
solid curves correspond to optical model calculations with fitting parameters
and are described in the next section.

Figure 4: (Color online) Elastic scattering angular distribution for the 10C +
208Pb system at E${}_{\textnormal{lab}}$ = 66 MeV (black circle points). The
points in blue triangles correspond to data for 12C + 208Pb system at 64.9 MeV
from Ref. Santra  _et al._ (2001). The solid curves were obtained by fitting
the parameters of the optical potential.

##  III Theoretical Analysis

###  III.1 Optical model calculation

In a first approach, the angular distribution of the elastic scattering is
analyzed within the optical model (OM) using a Woods-Saxon (WS) shape for the
real and imaginary potentials. The depths, radii and diffuseness parameters
were adjusted to minimize the $\chi^{2}$ parameter between data and the OM
cross section. The parameters obtained in the fitting procedure for the 12C +
208Pb data at 64.9 MeV Santra  _et al._ (2001) were adopted as initial
guessing parameters for the 10C + 208Pb system. The best values obtained by
the fitting procedure ($\chi^{2}=48.2$) are shown in Table 1 along with the
corresponding reaction cross section, $\sigma_{R}=753$ mb.

Table 1: Parameters of the optical potentials obtained by fitting of the elastic scattering in the 12C + 208Pb data at Elab = 64.9 MeV and 10C + 208Pb data at Elab = 66 MeV. Parameters $V$ and $W$ are in MeV. The $r_{V}$ and $r_{W}$ are reduced radii, using the convention $R_{i}=r_{i}\times[A_{p}^{1/3}+A_{t}^{1/3}]$ fm and $r_{C}=1.3$ fm. The $a_{V}$ and $a_{W}$ are the diffuseness, in fm. The reaction cross sections ($\sigma_{R}$) obtained with the optical potentials are also presented. projectile | $V$ | $r_{V}$ | $a_{V}$ | $W$ | $r_{W}$ | $a_{W}$ |  $\sigma_{R}$ (mb)
---|---|---|---|---|---|---|---
12C | 83.7 | 1.21 | 0.63 | 33.7 | 1.37 | 0.24 | 382
10C | 82.2 | 1.19 | 0.12 | 17.6 | 1.60 | 0.14 | 753

To compare the $\sigma_{R}$ for different projectiles on 208Pb it is necessary
to adopt a reduction method which removes the geometrical and charged
properties of the systems. There are three different reduction methods for
reaction cross section available in the literature, which, following the
nomenclature used by Canto et al. Canto  _et al._ (2015b), are referred as:
traditional, simplified traditional Gomes  _et al._ (2005), and fusion
(reaction) function Canto  _et al._ (2009a); Shorto  _et al._ (2009). Although
none of this methods have acceptable performance for a wide variation of
systems, they can be reliable for reaction cross sections of similar systems
with same target Canto  _et al._ (2015b); Canto, L. F.  _et al._ (2020). Here
we adopt the fusion (reaction) function method. In this case, the reduced
energy ($E_{\rm Red}$) and reduced reaction cross section ($\sigma_{\rm Red}$)
are defined as follows:

$E_{\rm Red}=\frac{E-V_{\rm B}}{\hbar\omega},\qquad\sigma_{\rm
Red}=\sigma_{\rm R}\left[\frac{2E}{\hbar\omega R_{\rm B}^{2}}\right].$ (1)

In both equations, $R_{\rm B}$, $V_{\rm B}$ and $\hbar\omega$ are the radius,
height and width of a parabolic shape adjusted to the Coulomb barrier. Here we
considered the double folding São Paulo potential Chamon  _et al._ (2002) to
obtain these parameters. For the 10C projectile, the adjusted $R_{\rm B}$,
$V_{\rm B}$ and $\hbar\omega$ are 11.3 fm, 58.8 MeV and 5.2 MeV, respectively.
Also, it is important to correct the averaged energy for the target thickness,
according to Ref. Fisichella  _et al._ (2015). Assuming that energies of
particles that compose the incident beam has a gaussian distribution, with
standard deviation of 0.8 MeV, and the target has an uniform thickness, the
effective energy ($E_{\text{eff}}$) is 66.4 MeV. This value has been
considered to calculate the corresponding $E_{\rm Red}$.

Figure 5: (Color online) Reduced reaction cross section for the several
projectiles on 208Pb target. The universal fusion function (UFF) is also shown
in the plot as a lower bound to the reduced reaction cross sections.

In Fig. 5, we compared the reduced cross section for the 10C with other
projectiles in 208Pb target: 6,8He Acosta  _et al._ (2011); Sánchez-Benítez
_et al._ (2008); Marquínez-Durán  _et al._ (2016), 6,7,8Li Keeley  _et al._
(1994); Kolata  _et al._ (2002), 9,10Be Woolliscroft  _et al._ (2004); Kolata
_et al._ (2004), 8B Mazzocco  _et al._ (2019) and 12C Santra  _et al._ (2001).
The dashed curve in the figure corresponds to the Universal Fusion Function
(UFF), introduced in Refs. Canto  _et al._ (2009b, a). This function is
related to the fusion cross sections based on the Wong’s equation. It is,
thus, a lower bound for the reduced reaction cross sections. The reduced
reaction cross section for the 10C, obtained with the optical model, lies a
little above of the weakly-bound 6,8He and 8Li projectiles.

An enhancement of the reduced cross section, compared to the UFF curve,
indicates that other reaction channels, driven by large deformation and/or
cluster configuration, give a sizeable contribution. This is actually the case
for 6He and 8B, for which transfer and/or breakup are known to be relevant
reaction channels. For 10C projectile, inelastic scattering and breakup
reaction could be of importance in the description of the elastic scattering
on 208Pb target, even though the binding energies (for instance, 3.821 MeV for
the 8Be+$p$+$p$) are relatively higher as compared to those for exotic nuclei,
which are typically less than 1 MeV.

In the next sections the coupled-channels (CC) and the continuum-discretized
coupled-channels (CDCC) calculations are explored to understand the effects of
inelastic scattering and cluster structure in the elastic scattering of 10C.

###  III.2 The Coupled-Channels calculations

The coupling to the $2_{1}^{+}$ state in 10C has been show to be quite
important to describe the elastic scattering of the 10C+58Ni system at
$E_{\textnormal{lab.}}=35.3$ MeV Guimarães  _et al._ (2019). Here, we assess
this effect of couplings to the inelastic channel. The coupled-channels
calculations were performed using the São Paulo potential (SPP) Chamon  _et
al._ (2002) for the real part and a short-range Wood-Saxon shape for the
imaginary part of the optical potential. This short range imaginary potential
takes into account the absorption flux due to the fusion process, since this
channel is not explicitly included in the coupling scheme. The parameters for
this potential are $W=50$ MeV, $r_{W}=1.06$ fm and $a_{W}=0.2$ fm, where $W$,
$r_{W}$ and $a_{W}$ are the depth, reduced radius and diffuseness,
respectively. This methodology has been applied with success to describe the
elastic scattering of many systems Canto, L. F.  _et al._ (2014); Deshmukh
_et al._ (2015). In the coupling scheme, the $2_{1}^{+}$ in 10C was treated in
the independent model adopting the transition probability B(E2) from Ref.
Raman  _et al._ (2001). The coupled-equations were solved using the fresco
code Thompson (1988).

The comparison between data and coupled-channels calculations is shown in Fig.
6. The pronounced Fresnel peak present in the calculation is completely absent
in the data. This is an indication that additional imaginary potential at the
surface is necessary to reproduce the data. The coupling to the strongly
deformed $2_{1}^{+}$ state in the projectile (solid orange curve) is
responsible for a shift and small reduction of the Fresnel peak and an
increase of the elastic cross sections at backward angles compared to the no-
coupling calculation (dotted black curve). This effect is also observed in the
10C+58Ni system Guimarães  _et al._ (2019). The inclusion of excited states of
the target in the coupling scheme does not change significantly the results
presented here (deviations less than 3%).

Figure 6: Comparison between data and coupled-channels calculations for the
elastic angular distribution for the 10C + 208Pb system at Elab = 66 MeV. The
no-coupling and the coupled-channel calculations are represented by dotted
black and solid orange curves, respectively.

The effect of transfer channels, such as neutron, proton and alpha, has been
investigated for the 10C+58Ni system Guimarães  _et al._ (2019) and little
influence was observed. Here, we have considered the inclusion of one- and
two-neutron pick-up reaction into the coupling scheme. The effect of such
channels is small and only sizeable at large scattering angles and, therefore,
not shown here.

###  III.3 The three-body continuum-discretized coupled-channels calculations

To investigate the effect of the continuum we performed continuum discretized
coupled-channels calculations (CDCC). As mentioned before, the 10C nucleus may
be described as $\alpha+\alpha+p+p$ cluster configuration. Such configuration
would require a 5-body CDCC calculation (4 bodies in the projectile and the
target nucleus), which is beyond the limits of contemporary reaction models.
However, there is an indication that the breakup of 10C proceeds through
sequential decays with three possible intermediate channels Curtis  _et al._
(2008): 8Be+$p$+$p$, 9B+$p$ and 6Be+$\alpha$ with binding energies of 3.821,
4.006 and 5.101 MeV, respectively. In this section we consider the effect of
the two-body 9B+$p$, and 6Be+$\alpha$ cluster configurations in the CDCC
calculations for the 10C + 208Pb elastic scattering.

Although both 9B and 6Be cores are unbound (9B $\rightarrow$
$\alpha$+$\alpha$+$p$ and 6Be $\rightarrow$ $\alpha$+$p$+$p$), we are
implicitly assuming that these cores remain bound during the reaction and
decay well away from the region of interaction with the target nucleus. This
hypothesis is particularly an over-simplified picture in the case of 6Be, for
which the ground state half-life is 92 keV ($\approx 7\times 10^{-21}$ s) and
comparable to the interaction time (roughly $10^{-21}$ s). However, it can be
useful to assess the influence of two-body configurations in the scattering.

The projectile-target relative-motion wave functions were computed considering
partial waves up to 1000$\hbar$ and radii up to $\mathit{R_{coupl}}$= 500 fm.
The energy bins with relative angular momentum of the fragments up to $4\hbar$
were generated by the superposition of scattering states, and potential
multipoles $\lambda$ up to the 5-pole term were considered. The energy bins
were equally spaced in energy up to the maximum value of 10.0 MeV and width of
3.0 MeV, for the 9B+p cluster configuration, and maximum value of 8.0 MeV and
width of 3.0 MeV for the 6Be+$\alpha$ cluster configuration. These values were
varied until a good convergence was achieved.

To reproduce the structure of 10C in the two-body model, the
${}^{9}\text{B}+p$ and ${}^{6}\text{Be}+\alpha$ potentials are required. In
both cases we have adjusted a Woods-Saxon potential to reproduce the energy
and matter radius as obtained in the three-body model for 10C (to be seen in
the next subsection). In this way we can compare the results from the 3b- and
4b-CDCC. The shape of the adjusted potential is a Woods-Saxon form for central
part with the radius fixed to 2.5 fm. For the ${}^{6}\text{Be}+\alpha$, the
strength and diffuseness of the central part was 85 MeV and 0.8 fm. For the
${}^{9}\text{B}+p$ we also included a derivative Woods-Saxon potential for the
spin-orbit. The strengths were 15 and 48 MeV for the spin-orbit and the
central part, respectively. The diffuseness was 1.2 fm in both parts. The
calculated ground-state matter radius are $r_{\rm m}=2.73$ fm and $r_{\rm
m}=2.61$ fm for the 9B+p and 6Be+$\alpha$ cluster configurations,
respectively, for which it was assumed a radius of 0.87 fm for $p$ and 1.47 fm
for $\alpha$. The SPP was used as the complex optical potential between the
fragments and the 208Pb target. The real and imaginary parts were multiplied
by a strength factor of 1.00 and 0.78, respectively.

The results of the 3-body calculations are presented in Fig. 7. The dotted
black curve corresponds to the calculation with no-couplings, using a
projectile-target double folding potential, as discussed in the previous
subsection. The dashed green and dot-dashed blue curves correspond to a
folding cluster model for the 10C as one-channel 9B+$p$ and 6Be+$\alpha$,
respectively, and considering only the ground state. In this model, couplings
to the continuum states are switched off and the optical potential is defined
as the expectation values of the sum of the fragments-target interactions over
the projectile wave function. This optical potential is complex and the
imaginary part is responsible for the suppression of the Fresnel peak. More
details about this one-channel calculation for 10C can be found in Ref.
Guimarães  _et al._ (2019). The inclusion of the $2_{1}^{+}$ state in 10C and
continuum couplings are represented by solid blue and dot-dot-dashed green
curves in Fig. 7. The difference between 1-channel (no coupling to continuum)
and full couplings is very small for both 3-body configurations. The major
contribution of the full coupling, with respect to the one-channel
calculation, comes from coupling to inelastic channel.

Figure 7: Comparison between data and 3b-CDCC calculations considering the 9B
+ p (green curves) and 6Be + $\alpha$ (blue curves) cluster configurations in
(a) linear scale and (b) log scale. The no-coupling calculation, presented in
the previous subsection, is also included (dotted black curve). See text for
details.

It is worth to highlight that, in this work, excited states are treated in a
different fashion in the CC and CDCC formalisms. In the first, the deformation
of the nuclei, which is the main ingredient to build the coupling potentials,
is obtained within a collective model. In the CDCC formalism, excited states
are built as single-particle states of a valence particle binding to its core.
We also checked the effect of different single-particle configuration in the
9B+$p$ cluster. In this case, the excited state can be build as the valence
particle in the $1p_{1/2}$ or a $1p_{3/2}$. The results for elastic and
inelastic cross sections are indistinguishable between these single-particle
configurations. The 2-body description for the $2_{1}^{+}$ state results in a
less deformed state and, consequently, its effect on the elastic angular
distribution is attenuated compared to the CC calculation (see Fig. 6). In
these cases, the deformation of the cores could be an important feature that
shall be included in this type of calculations.

In Ref. Curtis  _et al._ (2008), it has been suggested that 9B+$p$ is the
dominant configuration even though the ground and excited states in 10C shall
exhibit a mixture of cluster configurations. A full 3b-CDCC calculation, where
both cluster configurations are simultaneously considered in the coupling
scheme, is beyond the current possibilities. Keeping in mind that we are
dealing with pure cluster configurations, the 3b-CDCC calculation considering
the 9B+$p$ configuration gives a slightly better description compared to the
6Be+$\alpha$ channel, but the overall agreement between theoretical curves and
the data is limited.

###  III.4  The four-body continuum-discretized coupled-channels calculations

Going a step further in the CDCC calculations, we have performed four-body
CDCC calculation Rodríguez-Gallardo  _et al._ (2008, 2009), in which the 10C
is described within a 8B+$p$+$p$ configuration. The states of 10C are computed
using a pseudostate approach with the analytical Transformed Harmonic
Oscillator (THO) method presented in Ref. Casal  _et al._ (2013) for three-
body projectiles. This formalism has been applied to reactions induced by the
Borromean nucleus 9Be on different targets Casal  _et al._ (2015); Arazi  _et
al._ (2018), showing the validity of the method to describe reactions induced
by three-body projectiles with more than one charged particle.

To reproduce the structure of 10C in the three-body model,
${}^{8}\text{Be}+p+p$, the ${}^{8}\text{Be}+p$ and $p+p$ potentials are
required. For the last, we have taken the Gogny-Pires-Tourreil (GPT) p-p
potential Gogny  _et al._ (1970) whereas for the ${}^{8}\text{Be}+p$ subsystem
we have adjusted a potential to reproduce the low-energy known states for 9B
Tilley  _et al._ (2004), shown in table 2, for $p$ and $d$ waves. For
$s$-waves we use the same as for $d$-waves. The shape of the adjusted
potential is a Woods-Saxon form for central part and the derivative Woods-
Saxon for the spin-orbit part. The values of the radius and the diffuseness
were fixed to 2.5 fm and 0.65 fm. The spin-orbit strength was 22 MeV and the
strength for the central part was 68.4 MeV for $s$\- and $d$-waves and 38.5
MeV for $p$-waves. Apart from the two-body potentials, we include an effective
three-body force to reproduce the energies of the bound states of 10C. The
three-body force used has a Gaussian form depending on the hyper-radius $\rho$
(the three-body problem is treated in hyperspherical coordinates, see for
example Rodríguez-Gallardo  _et al._ (2008) for details):

$V_{3}(\rho)=v_{3}\exp{\left[\frac{\rho}{R_{3}}\right]}.$ (2)

The radius $R_{3}$ was fixed to 5 fm, meanwhile the strength $v_{3}$ is
adjusted to obtain the correct energies for bound states. This is fixed to
-2.5 MeV for $j^{\pi}=0^{+}$ and 0.33 MeV for $j^{\pi}=2^{+}$. For
$j^{\pi}=1^{-}$, we use the same strength as for $j^{\pi}=0^{+}$. To block the
Pauli forbidden states we use the adiabatic projection method Thompson  _et
al._ (2000).

Table 2: 9B energy states with respect to the 8Be$+p$ threshold. $j^{\pi}$ |  $E_{x}$ (MeV)
---|---
$3/2^{-}$ | 0.1851
$5/2^{+}$ | 2.5649
$1/2^{-}$ | 2.6029

The parameters of the analytical THO basis chosen are fixed to $m=4$, $b=0.7$
and $\gamma=1.2$ (see Ref. Casal  _et al._ (2013)). The maximum hypermomentum
is set to $K_{\rm max}=8$, which has been checked to provide converged results
for reaction calculations at the range of energies considered. The convergence
is also reached using a THO basis with $i_{\rm max}=8$ hyper-radial
excitations. The calculated ground-state energy is $\varepsilon_{B}=-3.798$
MeV (fairly close to the binding energy 3.821 MeV) and rms matter radius
$r_{\rm mat}=2.651$ fm, assuming a radius of 2.5 fm for 8Be.

Then, the 10C$-^{208}$Pb four-body wave functions are expanded in the internal
states of the three-body projectile, leading to a coupled-equations system
that has to be solved. For that, a multipole ($Q$) expansion is performed for
each projectile fragment-target interacting potential. The procedure is
explained in detail in Ref. Rodríguez-Gallardo  _et al._ (2008). The
interaction between each projectile fragment and the target is represented by
a complex optical potential, with both Coulomb and nuclear contributions. For
consistency with the previous subsections, the complex optical potentials
between the fragments and the 208Pb are the SPP with the strength factors
indicated in the previous sections. In this way we can compare the results
between the 3b-CDCC and 4b-CDCC on the same basis. Our model space includes
$j^{\pi}=0^{+},1^{-},2^{+}$ states up to 8 MeV above the breakup threshold,
which ensures convergence of the elastic angular distributions for this
reaction. The coupled equations are solved up to 100 partial waves, including
continuum couplings to all multipole orders, i.e., up to $Q=4$.

The 4b-CDCC result for the elastic cross section is shown in Fig. 8. The one-
channel calculation is represented by the dashed red curve and the inclusion
of couplings to the $2_{1}^{+}$ state and to the continuum (solid red curve)
just provides a small increase in the cross sections at backward angles, that
is mainly due to the inelastic coupling. For the sake of comparison, we also
included the full 3b-CDCC calculation within the 9B+$p$ configuration. The 3b-
and 4b-CDCC results are quite similar and although both of them underpredicts
the cross sections at the backward angles. The agreement with the data in the
region of Fresnel peak is reasonable. We have to consider that a more
realistic description of this nucleus would require a five-body model
($\alpha+\alpha+p+p$). Also, we stress that double-folding SPP, that is
adopted for the complex potentials, can be regarded as energy independent at
energies around the coulomb barrier. The relevant feature for the calculations
is its double-folding characteristic and there were no free adjustable
parameters in these calculations. We have considered other optical potentials
in the CDCC calculations and results are very similar.

Figure 8: Comparison between data and 4b-CDCC calculation considering the 10C
nucleus has a 8Be+p+p cluster. The one-channel and full couplings results are
represented by dashed and solid red curves, respectively. For reference, the
no-coupling (dotted black curve) and the 3b-CDCC assuming the 9B+$p$
configuration (dot-dashed green curve) are also included. Figure 9: Reaction
cross sections as a function of the total angular momentum ($J$) as obtained
from optical model fitting (OM), 3b-CDCC (6Be+$\alpha$ and 9B+$p$) and 4b-CDCC
(8Be+p+p) calculations.

The comparison of reaction cross sections as a function of the total angular
momentum of the scattering ($J$) between the OM, the 3b- and 4b-CDCC
calculations is shown in Fig. 9. The shape of the reaction cross sections
(Fig. 9) as obtained from the OM fitting (dotted black solid), the 9B+$p$
(dashed blue), the 6Be+$\alpha$ (dash-dotted green) and the 8Be+$p$+$p$
cluster (solid red) are all similar. We note that the maximum partial reaction
cross section in the OM is at slightly higher $J$ ($\approx 21\hbar$) than in
the CDCC calculations ($\approx 18\hbar$). The integrated cross sections are
listed in Table 3 and are close to the value estimated in the optical model
fitting ($\sigma_{R}$ = 753 mb). However, the behaviour of the reaction cross
sections as obtained in the OM analysis is different than the ones obtained
with the CDCC calculation.

Table 3: Total reaction cross sections ($\sigma_{\textnormal{total}}$) obtained from the OM and CDCC calculations. model |  $\sigma_{\textnormal{total}}$ (mb)
---|---
OM | 753
${}^{9}\textnormal{B}+p$ | 699
${}^{6}\textnormal{Be}+\alpha$ | 715
${}^{8}\textnormal{Be}+p+p$ | 727

##  IV Conclusions

The angular distribution for the elastic scattering of 10C on 208Pb target has
been measured for the first time at an energy close to the Coulomb barrier.
The experimental angular distribution was analyzed in terms of optical model,
coupled channels and CDCC calculations. The total reaction cross section
derived from the optical model analysis was compared to those from other
projectiles, with similar masses, on the 208Pb target. The reduced reaction
cross section is a little higher than those for 6He and 8Li projectiles in the
same nuclei target. The enhanced reduced reaction cross section, compared to
the UFF curve, and the absence of a Fresnel peak in the data indicate a strong
absorption at the surface.

From the coupled channel calculation, the couplings to the $2^{+}_{1}$ state
in the 10C, treated as a collective state, is an important ingredient to
describe the data at large scattering angles. The cluster configurations for
the 10C has been studied within the 3b-CDCC, assuming the 9B+$\mathit{p}$ and
6Be+$\alpha$ configurations, and within the 4b-CDCC calculation where the 10C
is described by the 8Be+$p$+$p$ configuration. Even though the binding
energies of these configurations are high (more than 3.8 MeV), the
calculations provide a reasonable description of the data, including the
prediction of the suppression of the Fresnel peak but underestimate the data
at backward angles.

It would be interesting to consider a four-body configuration for the 10C
which naturally will embrace clusters explored in this work and could provide
a better description of the $2_{1}^{+}$ state. However, such calculation are
still beyond the present capability of the reaction theory.

## Acknowledgment

The authors would like to thank the following foundations for financial
support: CNPq (Grants 439375/2018-5 and 304961/2017-5), São Paulo Research
Foundation (FAPESP) (Grants 2016/02863-4 and 2016/17612-7), FAPERJ, CAPES
(Grant PNPD 88887.475459/2020-00), INCT-FNA (Instituto Nacional de Ciência e
Tecnologia - Física Nuclear e Aplicações) and Department of Energy, Office of
Science, under Award No. DE-FG02-93ER40773. Also we would like to thank the
operators team at the Cyclotron Institute of Texas A&M University for reliable
operation of the Cyclotron during the data acquisition that generated this
work. This work was partially supported by the Spanish Ministry of Economy and
Competitiveness, the European Regional Development Fund (FEDER), under Project
$\rm N^{o}$ FIS2017-88410-P, by the European Union’s Horizon 2020 research and
innovation program, under Grant Agreement $\rm N^{o}$ 654002 and by SID funds
2019 (Università degli Studi di Padova, Italy) under project No.
CASA_SID19_01.

## References

  * Kolata, J. J.  _et al._ (2016) Kolata, J. J., Guimarães, V., and Aguilera, E. F., Eur. Phys. J. A 52, 123 (2016).
  * Canto  _et al._ (2015a) L. Canto, P. Gomes, R. Donangelo, J. Lubian, and M. Hussein, Physics Reports 596, 1 (2015a), recent developments in fusion and direct reactions with weakly bound nuclei.
  * Canto, L. F.  _et al._ (2020) Canto, L. F., Guimarães, V., Lubian, J., and Hussein, M. S., Eur. Phys. J. A in press (2020).
  * Tanihata  _et al._ (1985) I. Tanihata, H. Hamagaki, O. Hashimoto, Y. Shida, N. Yoshikawa, K. Sugimoto, O. Yamakawa, T. Kobayashi, and N. Takahashi, Phys. Rev. Lett. 55, 2676 (1985).
  * Tanihata  _et al._ (2013) I. Tanihata, H. Savajols, and R. Kanungo, Progress in Particle and Nuclear Physics 68, 215 (2013).
  * Austern  _et al._ (1987) N. Austern, Y. Iseri, M. Kamimura, M. Kawai, G. Rawitscher, and M. Yahiro, Physics Reports 154, 125 (1987).
  * Morcelle  _et al._ (2014) V. Morcelle, K. C. C. Pires, M. Rodríguez-Gallardo, R. Lichtenthäler, A. Lépine-Szily, V. Guimarães, P. N. Faria, D. R. Mendes Junior, A. M. Moro, L. R. Gasques, E. Leistenschneider, R. Pampa Condori, V. Scarduelli, M. C. Morais, A. Barioni, I. C. Zamora, and J. M. B. Shorto, Phys. Lett. B 732, 228 (2014).
  * Di Pietro  _et al._ (2004) A. Di Pietro, P. Figuera, F. Amorini, C. Angulo, G. Cardella, S. Cherubini, T. Davinson, D. Leanza, J. Lu, H. Mahmud, M. Milin, A. Musumarra, A. Ninane, M. Papa, M. G. Pellegriti, R. Raabe, F. Rizzo, C. Ruiz, A. C. Shotter, N. Soic, S. Tudisco, and L. Weissman, Phys. Rev. C 69, 044613 (2004).
  * de Faria  _et al._ (2010) P. N. de Faria, R. Lichtenthäler, K. C. C. Pires, A. M. Moro, A. Lépine-Szily, V. Guimarães, D. R. Mendes, A. Arazi, M. Rodríguez-Gallardo, A. Barioni, V. Morcelle, M. C. Morais, O. Camargo, J. Alcantara Nuñez, and M. Assunção, Phys. Rev. C 81, 044605 (2010).
  * Sánchez-Benítez  _et al._ (2008) A. M. Sánchez-Benítez, D. Escrig, M. A. G. Álvarez, M. V. Andrés, C. Angulo, M. J. G. Borge, J. Cabrera, S. Cherubini, P. Demarelt, J. M. Espino, P. Figuera, M. Freer, J. E. Garí a Ramos, J. Gómez-Camacho, M. Gulino, O. R. Kakuee, I. Martel, C. Metelco, A. M. Moro, F. Pŕez-Bernal, J. Rahigi, K. Rusek, D. Smirnov, O. Tengblad, P. Van Duppen, and V. Ziman, Nucl. Phys. A 803, 30 (2008).
  * Acosta  _et al._ (2011) L. Acosta, A. M. Sánchez-Benítez, M. E. Gómez, I. Martel, F. Pérez-Bernal, F. Pizarro, J. Rodríguez-Quintero, K. Rusek, M. A. G. Alvarez, M. V. Andrés, J. M. Espino, J. P. Fernández-García, J. Gómez-Camacho, A. M. Moro, C. Angulo, J. Cabrera, E. Casarejos, P. Demaret, M. J. G. Borge, D. Escrig, O. Tengblad, S. Cherubini, P. Figuera, M. Gulino, M. Freer, C. Metelko, V. Ziman, R. Raabe, I. Mukha, D. Smirnov, O. R. Kakuee, and J. Rahighi, Phys. Rev. C 84, 044604 (2011).
  * Aguilera  _et al._ (2011) E. F. Aguilera, P. Amador-Valenzuela, E. Martinez-Quiroz, D. Lizcano, P. Rosales, H. García-Martínez, A. Gómez-Camacho, J. J. Kolata, A. Roberts, L. O. Lamm, G. Rogachev, V. Guimarães, F. D. Becchetti, A. Villano, M. Ojaruega, M. Febbraro, Y. Chen, H. Jiang, P. A. De Young, G. F. Peaslee, C. Guess, U. Khadka, J. Brown, J. D. Hinnefeld, L. Acosta, E. S. Rossi Jr, J. F. P. Huiza, and T. L. Belyaeva, Phys. Rev. Lett. 107, 092701 (2011).
  * Di Pietro  _et al._ (2010) A. Di Pietro, G. Randisi, V. Scuderi, L. Acosta, F. Amorini, M. J. G. Borge, P. Figuera, M. Fisichella, L. M. Fraile, J. Gomez-Camacho, H. Jeppesen, M. Lattuada, I. Martel, M. Milin, A. Musumarra, M. Papa, M. G. Pellegriti, F. Perez-Bernal, R. Raabe, F. Rizzo, D. Santonocito, G. Scalia, O. Tengblad, D. Torresi, A. M. Vidal, D. Voulot, F. Wenander, and M. Zadro, Phys. Rev. Lett 105, 022701 (2010).
  * Di Pietro  _et al._ (2012) A. Di Pietro, V. Scuderi, A. M. Moro, L. Acosta, F. Amorini, M. J. G. Borge, P. Figuera, M. Fisichella, L. M. Fraile, J. Gomez-Camacho, H. Jeppesen, M. Lattuada, I. Martel, M. Milin, A. Musumarra, M. Papa, M. G. Pellegriti, F. Perez-Bernal, R. Raabe, G. Randisi, F. Rizzo, G. Scalia, O. Tengblad, D. Torresi, A. M. Vidal, D. Voulot, F. Wenander, and M. Zadro, Phys. Rev. C 85, 054607 (2012).
  * Acosta  _et al._ (2009) L. Acosta, M. Álvarez, M. Andrés, M. J. G. Borge, M. Cortés, J. M. Espino, D. Galaviz, J. Gómez-Camacho, A. Maira, I. Martel, A. M. Moro, I. Mukha, F. Pérez-Bernal, E. Reillo, D. Rodríguez, K. Rusek, A. M. Sánchez-Benítez, , and O. Tengblad, Eur. Phys. J. A 42, 461 (2009).
  * Pesudo  _et al._ (2017) V. Pesudo, M. J. G. Borge, A. M. Moro, J. A. Lay, E. Nácher, J. Gómez-Camacho, O. Tengblad, L. Acosta, M. Alcorta, M. A. G. Alvarez, C. Andreoiu, P. C. Bender, R. Braid, M. Cubero, A. Di Pietro, J. P. Fernández-García, P. Figuera, M. Fisichella, B. R. Fulton, A. B. Garnsworthy, G. Hackman, U. Hager, O. S. Kirsebom, K. Kuhn, M. Lattuada, G. Marquínez-Durán, I. Martel, D. Miller, M. Moukaddam, P. D. O’Malley, A. Perea, M. M. Rajabali, A. M. Sánchez-Benítez, F. Sarazin, V. Scuderi, C. E. Svensson, C. Unsworth, and Z. M. Wang, Phys. Rev. Lett. 118, 152502 (2017).
  * Mazzocco  _et al._ (2006) M. Mazzocco, C. Signorini, M. Romoli, A. De Franesco, M. Di Pietro, E. Vardaci, K. Yoshida, R. Yoshida, A. Bonetti, A. De Rosa, T. Glodariu, A. Guglielmetti, G. Inglima, M. La Commara, B. Martin, D. Pierroutsakou, M. Sandoli, F. Soramel, L. Stroe, R. Kanungo, N. Khai, T. Motobayshi, T. Nomura, T. Ishikawa, H. Ishiyama, S. Jeong, H. Miyatake, M. H. Tanaka, I. Sugai, and Y. Watanabe, Eur. Phys. J. A 28, 295 (2006).
  * Rodríguez-Gallardo  _et al._ (2009) M. Rodríguez-Gallardo, J. M. Arias, J. Gómez-Camacho, A. M. Moro, I. J. Thompson, and J. A. Tostevin, Phys. Rev. C 80, 051601(R) (2009).
  * Casal  _et al._ (2015) J. Casal, M. Rodríguez-Gallardo, and J. M. Arias, Phys. Rev. C 92, 054611 (2015).
  * Arazi  _et al._ (2018) A. Arazi, J. Casal, M. Rodríguez-Gallardo, J. M. Arias, R. Lichtenthäler Filho, D. Abriola, O. A. Capurro, M. A. Cardona, P. F. F. Carnelli, E. de Barbará, J. Fernández Niello, J. M. Figueira, L. Fimiani, D. Hojman, G. V. Martí, D. Martínez Heimman, and A. J. Pacheco, Phys. Rev. C 97, 044609 (2018).
  * Cubero  _et al._ (2012) M. Cubero, J. P. Fernández-García, M. Rodríguez-Gallardo, L. Acosta, M. Alcorta, M. A. G. Alvarez, M. J. G. Borge, L. Buchmann, C. A. Diget, H. A. Falou, B. R. Fulton, H. O. U. Fynbo, D. Galaviz, J. Gómez-Camacho, R. Kanungo, J. A. Lay, M. Madurga, I. Martel, A. M. Moro, I. Mukha, T. Nilsson, A. M. Sánchez-Benítez, A. Shotter, O. Tengblad, and P. Walden, Phys. Rev. Lett. 109, 262701 (2012).
  * Fernández-García  _et al._ (2013) J. P. Fernández-García, M. Cubero, M. Rodríguez-Gallardo, L. Acosta, M. Alcorta, M. A. G. Alvarez, M. J. G. Borge, L. Buchmann, C. A. Diget, H. A. Falou, B. R. Fulton, H. O. U. Fynbo, D. Galaviz, J. Gómez-Camacho, R. Kanungo, J. A. Lay, M. Madurga, I. Martel, A. M. Moro, I. Mukha, T. Nilsson, A. M. Sánchez-Benítez, A. Shotter, O. Tengblad, and P. Walden, Phys. Rev. Lett. 110, 142701 (2013).
  * Aguilera  _et al._ (2009) E. F. Aguilera, E. Martinez-Quiroz, D. Lizcano, A. Gómez-Camacho, J. J. Kolata, L. O. Lamm, V. Guimarães, R. Lichtenthäler, O. Camargo, F. D. Becchetti, H. Jiang, P. A. DeYoung, P. J. Mears, and T. L. Belyaeva, Phys. Rev. C 79, 021601(R) (2009).
  * Lubian  _et al._ (2009) J. Lubian, T. Correa, E. F. Aguilera, L. F. Canto, A. Gomez-Camacho, E. M. Quiroz, and P. R. S. Gomes, Phys. Rev. C 79, 064605 (2009).
  * Keeley  _et al._ (2010) N. Keeley, R. Mackintosh, and C. Beck, Nuclear Physics A 834, 792c (2010), the 10th International Conference on Nucleus-Nucleus Collisions (NN2009).
  * Guimaraes  _et al._ (2000) V. Guimaraes, J. J. Kolata, D. Peterson, P. Santi, R. H. White-Stevens, S. M. Vincent, F. D. Becchetti, M. Y. Lee, T. W. O’Donnell, D. A. Roberts, and J. A. Zimmerman, Phys. Rev. Lett. 84, 1862 (2000).
  * Yang  _et al._ (2013) Y. Y. Yang, J. S. Wang, Q. Wang, D. Pang, J. B. Ma, M. R. Huang, J. L. Han, P. Ma, S. L. Jin, Z. Bai, Q. Hu, L. Jin, J. B. Chen, N. Keeley, K. Rusek, R. Wada, S. Mukherjee, Z. Y. Sun, R. F. Chen, X. Y. Zhang, Z. G. Hu, X. H. Yuan, X. G. Cao, Z. G. Xu, S. W. Xu, C. Zhen, Z. Q. Chen, Z. Chen, S. Z. Chen, C. M. Du, L. M. Duan, F. Fu, B. X. Gou, J. Hu, J. J. He, X. G. Lei, S. L. Li, Y. Li, Q. Y. Lin, L. X. Liu, F. D. Shi, S. W. Tang, G. Xu, X. Xu, L. Y. Zhang, X. H. Zhang, W. Zhang, M. H. Zhao, Z. Y. Guo, Y. H. Zhang, H. S. Xu, and G. Q. Xiao, Phys. Rev. C 87, 044613 (2013).
  * Rangel  _et al._ (2016) J. Rangel, J. Lubian, L. F. Canto, and P. R. S. Gomes, Phys. Rev. C 93, 054610 (2016).
  * Yang  _et al._ (2016) Y. Y. Yang, X. Liu, and D. Y. Pang, Phys. Rev. C 94, 034614 (2016).
  * Mazzocco  _et al._ (2019) M. Mazzocco, N. Keeley, A. Boiano, C. Boiano, M. La Commara, C. Manea, C. Parascandolo, D. Pierroutsakou, C. Signorini, E. Strano, D. Torresi, H. Yamaguchi, D. Kahl, L. Acosta, P. Di Meo, J. P. Fernandez-Garcia, T. Glodariu, J. Grebosz, A. Guglielmetti, Y. Hirayama, N. Imai, H. Ishiyama, N. Iwasa, S. C. Jeong, H. M. Jia, Y. H. Kim, S. Kimura, S. Kubono, G. La Rana, C. J. Lin, P. Lotti, G. Marquínez-Durán, I. Martel, H. Miyatake, M. Mukai, T. Nakao, M. Nicoletto, A. Pakou, K. Rusek, Y. Sakaguchi, A. M. Sánchez-Benítez, T. Sava, O. Sgouros, V. Soukeras, F. Soramel, E. Stiliaris, L. Stroe, T. Teranishi, N. Toniolo, Y. Wakabayashi, Y. X. Watanabe, L. Yang, Y. Y. Yang, and H. Q. Zhang, Phys. Rev. C 100, 024602 (2019).
  * Curtis  _et al._ (2008) N. Curtis, N. L. Achouri, N. I. Ashwood, H. G. Bohlen, W. N. Catford, N. M. Clarke, M. Freer, P. J. Haigh, B. Laurent, N. A. Orr, N. P. Patterson, N. Soić, J. S. Thomas, and V. Ziman, Phys. Rev. C 77, 021301(R) (2008).
  * Guimarães  _et al._ (2019) V. Guimarães, E. N. Cardozo, V. B. Scarduelli, J. Lubian, J. J. Kolata, P. D. O’Malley, D. W. Bardayan, E. F. Aguilera, E. Martinez-Quiroz, D. Lizcano, A. Garcia-Flores, M. Febbraro, C. C. Lawrence, J. Riggins, R. O. Torres-Isea, P. N. de Faria, D. S. Monteiro, E. S. Rossi, and N. N. Deshmukh, Phys. Rev. C 100, 034603 (2019).
  * Yang  _et al._ (2014) Y. Y. Yang, J. S. Wang, Q. Wang, D. Y. Pang, J. B. Ma, M. R. Huang, P. Ma, S. L. Jin, J. L. Han, Z. Bai, L. Jin, J. B. Chen, Q. Hu, R. Wada, S. Mukherjee, Z. Y. Sun, R. F. Chen, X. Y. Zhang, Z. G. Hu, X. H. Yuan, S. W. Xu, S. Z. Chen, X. G. Lei, L. X. Liu, W. H. Ma, S. T. Wang, D. Yan, X. H. Zhang, M. H. Zhao, Y. Zhou, Y. J. Zhou, Z. Y. Guo, Y. H. Zhang, H. S. Xu, and G. Q. Xiao, Phys. Rev. C 90, 014606 (2014).
  * Tribble  _et al._ (2002) R. Tribble, A. Azhari, C. Gagliardi, J. Hardy, A. Mukhamedzhanov, X. Tang, L. Trache, and S. Yennello, Nuclear Physics A 701, 278 (2002), 5th International Conference on Radioactive Nuclear Beams.
  * Tarasov and Bazin (2008) O. Tarasov and D. Bazin, Nucl. Instrum. Methods Phys. Res. Sect. B 266, 4657 (2008), proceedings of the XVth International Conference on Electromagnetic Isotope Separators and Techniques Related to their Applications.
  * Gavron (1980) A. Gavron, Phys. Rev. C 21, 230 (1980).
  * Santra  _et al._ (2001) S. Santra, P. Singh, S. Kailas, A. Chatterjee, A. Shrivastava, and K. Mahata, Phys. Rev. C 64, 024602 (2001).
  * Canto  _et al._ (2015b) L. F. Canto, D. R. M. Junior, P. R. S. Gomes, and J. Lubian, Phys. Rev. C 92, 014626 (2015b).
  * Gomes  _et al._ (2005) P. R. S. Gomes, J. Lubian, I. Padrón, and R. M. Anjos, Phys. Rev. C 71, 017601 (2005).
  * Canto  _et al._ (2009a) L. F. Canto, P. R. S. Gomes, J. Lubian, L. C. Chamon, and E. Crema, Nucl. Phys. A 821, 51 (2009a).
  * Shorto  _et al._ (2009) J. M. B. Shorto, P. R. S. Gomes, J. Lubian, L. F. Canto, S. Mukherjee, and L. C. Chamon, Phys. Letters B 678, 77 (2009).
  * Chamon  _et al._ (2002) L. C. Chamon, B. V. Carlson, L. R. Gasques, D. Pereira, C. De Conti, M. A. G. Alvarez, M. S. Hussein, M. A. Cândido Ribeiro, E. S. Rossi Jr., and C. P. Silva, Phys. Rev. C 66, 014610 (2002).
  * Fisichella  _et al._ (2015) M. Fisichella, A. C. Shotter, A. Di Pietro, P. Figuera, M. Lattuada, C. Marchetta, V. Privitera, L. Romano, C. Ruiz, and M. Zadro, Phys. Rev. C 92, 064611 (2015).
  * Marquínez-Durán  _et al._ (2016) G. Marquínez-Durán, I. Martel, A. M. Sánchez-Benítez, L. Acosta, R. Berjillos, J. Dueñas, K. Rusek, N. Keeley, M. A. G. Álvarez, M. J. G. Borge, A. Chbihi, C. Cruz, M. Cubero, J. P. Fernández-García, B. Fernández-Martínez, J. L. Flores, J. Gómez-Camacho, K. W. Kemper, J. A. Labrador, M. Marqués, A. M. Moro, M. Mazzocco, A. Pakou, V. V. Parkar, N. Patronis, V. Pesudo, D. Pierroutsakou, R. Raabe, R. Silvestri, N. Soic, L. Standyło, I. Strojek, O. Tengblad, R. Wolski, and Z. Abou-Haidar, Phys. Rev. C 94, 064618 (2016).
  * Keeley  _et al._ (1994) N. Keeley, S. J. Bennet, M. N. Clarke, B. R. Fulton, G. Tungate, P. V. Drumm, M. A. Nagarajan, and J. L. Lilley, Nucl. Phys. A571, 326 (1994).
  * Kolata  _et al._ (2002) J. J. Kolata, V. Z. Goldberg, L. O. Lamm, M. G. Marino, C. J. O’Keeffe, G. Rogachev, E. F. Aguilera, H. García-Martínez, E. Martinez-Quiroz, P. Rosales, F. D. Becchetti, T. W. O’Donnell, D. A. Roberts, J. A. Brown, P. A. DeYoung, J. D. Hinnefeld, and S. A. Shaheen, Phys. Rev. C 65, 054616 (2002).
  * Woolliscroft  _et al._ (2004) R. J. Woolliscroft, B. R. Fulton, R. L. Cowin, M. Dasgupta, D. J. Hinde, C. R. Morton, and A. C. Berriman, Phys. Rev. C 69, 044612 (2004).
  * Kolata  _et al._ (2004) J. J. Kolata, E. F. Aguilera, F. D. Becchetti, Y. Chen, P. A. DeYoung, H. García-Martínez, J. D. Hinnefeld, J. H. Lupton, E. Martinez-Quiroz, and G. Peaslee, Phys. Rev. C 69, 047601 (2004).
  * Canto  _et al._ (2009b) L. F. Canto, P. R. S. Gomes, J. Lubian, L. C. Chamon, and E. Crema, J. Phys. G: Nucl. Part. Phys. 36, 015109 (2009b).
  * Canto, L. F.  _et al._ (2014) Canto, L. F., Gomes, P. R. S., Lubian, J., Hussein, M. S., and Lotti, P., Eur. Phys. J. A 50, 89 (2014).
  * Deshmukh  _et al._ (2015) N. N. Deshmukh, V. Guimarães, E. Crema, D. Abriola, A. Arazi, E. de Barbará, O. A. Capurro, M. A. Cardona, J. Gallardo, D. Hojman, G. V. Martí, A. J. Pacheco, D. Rodrígues, Y. Y. Yang, A. N. Deshmukh, D. R. Mendes, V. Morcelle, V. Scarduelli, and D. S. Monteiro, Phys. Rev. C 92, 054615 (2015).
  * Raman  _et al._ (2001) S. Raman, C. Nestor, Jr., and P. Tikkanen, Atomic Data and Nuclear Data Tables 78, 1 (2001).
  * Thompson (1988) I. J. Thompson, Comput. Phys. Rep. 7, 167 (1988).
  * Rodríguez-Gallardo  _et al._ (2008) M. Rodríguez-Gallardo, J. M. Arias, J. Gómez-Camacho, R. C. Johnson, A. M. Moro, I. J. Thompson, and J. A. Tostevin, Phys. Rev. C 77, 064609 (2008).
  * Casal  _et al._ (2013) J. Casal, M. Rodríguez-Gallardo, and J. M. Arias, Phys. Rev. C 88, 014327 (2013).
  * Gogny  _et al._ (1970) D. Gogny, P. Pires, and R. de Tourreil, Phys. Lett. 32B, 591 (1970).
  * Tilley  _et al._ (2004) D. Tilley, J. Kelley, J. Godwin, D. Millener, J. Purcell, C. Sheu, and H. Weller, Nucl. Phys. A 745, 155 (2004).
  * Thompson  _et al._ (2000) I. J. Thompson, B. V. Danilin, V. D. Efros, J. S. Vaagen, J. M. Bang, and M. V. Zhukov, Phys. Rev. C 61, 024318 (2000).


# A hybrid physics-informed neural network for nonlinear partial differential
equation

Chunyue Lv  Lei Wang<EMAIL_ADDRESS><EMAIL_ADDRESS>Chenming Xie
School of Mathematics and Physics, China University of Geosciences, Wuhan
430074, China  Center for Mathematical Sciences, China University of
Geosciences, Wuhan 430074, China  Kindo medical data technology company
limited , Wuhan 430073, China

###### Abstract

The recently developed physics-informed machine learning has made great
progress for solving nonlinear partial differential equations (PDEs), however,
it may fail to provide reasonable approximations to the PDEs with
discontinuous solutions. In this paper, we focus on the discrete time physics-
informed neural network (PINN), and propose a hybrid PINN scheme for the
nonlinear PDEs. In this approach, the local solution structures are classified
as smooth and non-smooth scales by introducing a discontinuity indicator, and
then the automatic differentiation technique is employed for resolving smooth
scales, while an improved weighted essentially non-oscillatory (WENO) scheme
is adopted to capture discontinuities. We then test the present approach by
considering the viscous and inviscid Burgers equations , and it is shown that
compared with original discrete time PINN, the present hybrid approach has a
better performance in approximating the discontinuous solution even at a
relatively larger time step.

###### keywords:

Physics-informed neural network , WENO scheme , Smoothness indicator , Runge-
Kutta method

††journal: XXX

##  1 Introduction

Nonlinear partial differential equations (PDEs) have been the focus of many
studies due to their frequent appearance in many applications in physical,
biological, and social sciences [1]. Because of most nonlinear PDEs do not
have analytical solutions expect for some extremely simple situations, and
therefore, various numerical methods have been developed during the past
several decades, such as the finite-difference [2], finite-volume [2] and
finite-element [3], and the lattice Boltzmann methods [4], to name but a few.

Apart from the traditional numerical approaches, the rise of deep learning in
recent years offers an opportunity to develop surrogate model for solving
nonlinear PDEs [5, 6, 7, 8]. Among various neural network based PDEs solvers,
the physics-informed neural network (PINN) proposed by Raissi et al. [8] has
drawn considerable attention due to its effectiveness in solving both forward
and inverse nonlinear PDE problems as well as the straightforward
implementation [8, 9]. The core idea of PINN algorithm is to encode the
underlying physical laws into the loss function of the neural network such
that the governing equations are enforced by minimizing the residual loss
function with automatic differentiation [8, 9]. More recently, the PINN and
its variants have been successfully applied to solve various PDEs, including
the Navier-Stokes equation, the Kdv equations [8], heat transfer equations
[10], Fokker-Plank equation [11] and so on. Despite the great progress of
PINNs in solving nonlinear PDEs, the governing equations appeared in previous
works are usually the so-called diffusion-dominated equations, and it is known
that the solutions in these cases are usually smooth enough such that the
traditional PINNs are capable of achieving good predictive accuracy according
to the classical universal approximation theorem [12]. However, as pointed out
by some authors, it is a challenge for current PINNs to predict the solutions
of PDEs with advection-dominant character, in which the solutions often
develop discontinuities in some finite time even if the initial conditions is
smooth enough [13]. To address this issue, some enhanced PINNs are proposed
recently. By using more scattered points around the discontinuous regions, Mao
et al. [14] investigated the possibility of continuous time PINN in solving
the high-speed aerodynamic flows, and found that compared with random or
uniform training points, the predicted solutions with clustered training
points are more efficient in some cases. In addition, another work reported by
Fuks et al. [15] shows that the performance of PINN in solving hyperbolic PDEs
can be enhanced by adding a small diffusion term to the governing PDEs. More
recently, by combining the traditional finite volume methodology with PINNs,
Patel et al. [17] proposed a thermodynamically consistent PINNs for hyperbolic
systems, and the solutions of various PDEs like Euler, Buckley-Leverett
equations are well predicted with this method. Although there are some works
on utilization of PINNs for solving the PDEs with discontinuities, the neural
network appeared in some previous works either needs a priori knowledge of the
solution [14] or the performance of the PINNs largely depends on the neural
network architecture [15, 16].

In this work, we propose a hybrid PINN (hPINN) for nonlinear PDEs, and the
basic idea of this approach is to introduce a discontinuity indicator into the
neural network such that the local solution structures can be classified as
smooth and non-smooth scales. Then the automatic differentiation technique
(one of the most important techniques in deep learning field [18]) is employed
for resolving smooth scales, while the improved weighted essentially non-
oscillatory (WENO) scheme [19], which is a famous shock-capturing scheme, is
adopted to capture discontinuities. The reminder of the present paper is
organized as follows. In Section 2, the details of the proposed approach are
presented. Then, some results for viscous and inviscid Burgers equations are
presented in Section 3, and a brief conclusion is drawn in Section 4.

##  2 Methodology

In this work, we are interested in designing a hybrid PINN model for the
following nonlinear PDE,

${u_{t}}+f{\left(u\right)_{\emph{{x}}}}=\upsilon
g{\left(u\right)_{\emph{{xx}}}}+h\left({\emph{{x}},t}\right),$ (1)

in which $u$ is the scalar variable and usually a function of both time $t$
and space x, $h\left({\emph{{x}},t}\right)$ is the source term, and $\upsilon$
is the diffusion coefficient. In addition, the first derivations in Eq. (1)
are associated with convection while the second derivatives are responsible
for diffusion.

###  2.1 Hybrid physics-informed neural network

Figure 1: Schematic of the hPINN for solving the nonlinear PDEs, in which the
spatial vector ${x}$ is the input, $\sigma$ represents the activation
function, which is set to be the hyperbolic tangent activation function in
this work [8], and all the hyperparameters in the neural network are
initialized by the Xavier initialization method [21].

As presented by Raissi et al. [8], the generally used PINNs can be classified
as the continuous time model and the discrete time model. Compared with the
latter model, the continuous time model has drawn tremendous interest in the
PINN community [14, 15, 17, 16]. However, as pointed out by Raissi et al. [8],
in order to enforce the physical conservation laws , the continuous time model
needs a large amount of collocation points in the entire spatio-temporal
domain, and thus the training of the original continuous time model is usually
prohibitively expensive. Fortunately, this limitation could be well addressed
by employing the discrete time model [8]. In such a case, the present work
intends to solve the forward PDE problems by using the discrete time model,
and the model proposed here can be viewed as the variant of the original
discrete time model.

To have a better understand on the present approach, the original discrete
time PINN is first reviewed. The main idea of this model is to discrete the
the time operator in Eq. (1) by using the classical Runge-Kutta methods with
$q$ stage [8], then we can obtain the following equations

$\begin{split}&{u^{n+{c_{i}}}}\left(x\right)={u^{n}}\left(x\right)-\Delta
t\sum\limits_{j=1}^{q}{{a_{ij}}{\mathscr{N}}}\left[{{u^{n+{c_{j}}}}\left(x\right)}\right],\;i=1,\ldots,q,\\\
&{u^{n+1}}\left(x\right)={u^{n}}\left(x\right)-\Delta
t\sum\limits_{j=1}^{q}{{b_{j}}{\mathscr{N}}}\left[{{u^{n+{c_{j}}}}\left(x\right)}\right],\end{split}$
(2)

in which ${\mathscr{N}}\left[u\right]=f{\left(u\right)_{x}}-\upsilon
g{\left(u\right)_{xx}}-h\left(x\right)$,
${u^{n+{c_{i}}}}\left(x\right)=u\left({{t^{n}}+{c_{i}}\Delta t,x}\right)$, and
the above time discretization scheme may be implicit or explicit, depending on
the values of $\left\\{a_{ij},b_{j},c_{j}\right\\}$. In such a case, the PINN
model usually consists of two sub-networks, in which the first neural network
${{{NN}}(w,b)}$ takes x as input and outputs the values of u at different time
stages, i.e.,
$\left[{{u^{n+{c_{1}}}},\ldots,{u^{n+{c_{q}}}},{u^{n+1}}}\right]$, and they
are further fed into the second network to encode the governing equations as
well as the corresponding boundary conditions. For the forward problems
considered here, since the initial/boundary conditions as well as the
governing equations value are always known, the loss function in the original
discrete time PINN model is given by [8],

$\begin{split}&L={L_{PDE}}+{L_{BC}},\\\
&{L_{PDE}}={\sum\limits_{j=1}^{q}{\sum\limits_{i=1}^{{N_{n}}}{\left|{u_{j}^{n}\left({{x^{n,i}}}\right)-{u^{n,i}}}\right|}}^{2}}+{\left|{u_{q+1}^{n}\left({{x^{n,i}}}\right)-{u^{n,i}}}\right|^{2}},\\\
&{L_{BC}}=\sum\limits_{j=1}^{q}{{{\left|{u_{j}^{n}\left({{x^{b}}}\right)-{u^{b}}}\right|}^{2}}}+{\left|{u_{q+1}^{n}\left({{x^{b}}}\right)-{u^{b}}}\right|^{2}},\end{split}$
(3)

where $N_{n}$ is the number of the collocation points, which can be randomly
or uniformly sampled inside the computational domain,
$\left\\{{{x^{n,i}},{u^{n,i}}}\right\\}_{i=1}^{{N_{n}}}$ is the numerical data
obtained at time-step $t^{n}$, $x^{b}$ represents the collocation points on
the boundary, and $u_{i}^{n}$ is defined as

$\begin{split}&u_{i}^{n}={u^{n+{c_{i}}}}+\Delta
t\sum\limits_{j=1}^{q}{{a_{ij}}{\mathscr{N}}}\left[{{u^{n+{c_{j}}}}\left(x\right)}\right],\;i=1,\ldots,q,\\\
&u_{q+1}^{n}={u^{n+1}}+\Delta
t\sum\limits_{j=1}^{q}{{b_{j}}{\mathscr{N}}}\left[{{u^{n+{c_{j}}}}\left(x\right)}\right].\end{split}$
(4)

As presented by some previous works [15, 16, 17], although the original PINN
works well for various nonlinear PDEs, it may fail to provide reasonable
approximations to the PDEs with dominant convection. To address this issue, an
intuitive way is to employ some discontinuity-capturing methods in the neural
network. In this work, inspired by the weighted essentially nonoscillatory
schemes (WENO) [19, 20], which are a popular class of numerical methods in
capturing discontinuities, an improved WENO scheme [19], dubbed WENO-Z, is
incorporate into the PINN. On the basis of this idea, a simple way is to
discretize the differential operator in the whole computational domain by
using WENO-Z approach. However, this treatment may lead to another two issues:
For one thing, when applying the WENO-Z scheme in the whole domain, the
training cost of the neural network will be increased significantly due to the
local characteristic decomposition and the nonlinear-weights computing [19,
20]. For another, the advantage of the automatic differentiation technique
used in original PINN will be weakened. In this setting, an alternative
solution is to introduce the hybrid concept, for which the efficient automatic
differentiation technique is adopted to compute all the differential operators
of the PDEs in smooth regions while the expensive WENO-Z scheme is just
applied in the vicinity of discontinuities, and it is observed that this idea
is widely employed in scientific computing community [22].

Fig. 1 shows a schematic of the present _hybrid_ PINN (hPINN). Note that the
evolution of the discontinuities depends little on the viscous term (i.e.,
$g(u)_{xx}$), the automatic differentiation technique is adopted to compute
the second derivative in both the smooth and non-smooth regions for its
efficient. For the approximation of the first order derivative ( i.e.,
$f(u)_{x}$) in the non-smooth scale, the fifth-order WENO-Z scheme is used.
For brevity in the presentation, we use one-dimensional scalar nonlinear PDEs
as an example, and assume the scatter points $x_{j}$ are uniform. The cell
size as well as the cells are denoted by $\Delta x={x_{j+1}}-{x_{j}}$ and
${I_{j}}=\left[{{x_{j-1/2}},{x_{j+1/2}}}\right]$ with
${x_{j+1/2}}={x_{j}}+\Delta x/2$, respectively. Then use a conservative finite
difference scheme to the first-order derivative term in Eq. (1), thus obtain

$\frac{{\partial f}}{{\partial
x}}\left|{{}_{x={x_{i}}}}\right.=\frac{{{{\hat{f}}_{i+1/2}}-{{\hat{f}}_{i-1/2}}}}{{\Delta
x}},$ (5)

in which the numerical fluxes ${{\hat{f}}_{i+1/2}}$ and ${{\hat{f}}_{i-1/2}}$
are the positive and negative parts of $f(u)$ at the cell boundaries. To
ensure the numerical stability, the global Lax-Friedrichs flux splitting is
used for calculating the flux $f(u)$

$f\left(u\right)={f^{+}}\left(u\right)+{f^{-}}\left(u\right),$ (6)

where ${f^{\pm}}\left(u\right)=\frac{1}{2}\left({f\left(u\right)\pm\lambda
u}\right)$ and $\lambda=\max\left|{f^{\prime}\left(u\right)}\right|$, and with
this, we have
${{\hat{f}}_{i+1/2}}={{\hat{f}}^{+}}_{i+1/2}+{{\hat{f}}^{-}}_{i+1/2}$. Since
the negative part of the flux ${{\hat{f}}_{i-1/2}}$ is symmetric to the
positive part with respect to $x_{i+\frac{1}{2}}$, here we will only describe
how ${{\hat{f}}^{+}}_{i+1/2}$ is constructed, and for convenience, we will
drop the ”+” sign in the superscript. As presented in [19], the fifth-order
WENO flux ${{\hat{f}}_{i+1/2}}$ is defined as

${{\hat{f}}_{i+1/2}}={\omega_{0}}{{\hat{f}}^{\left(0\right)}}_{i+1/2}+{\omega_{1}}{{\hat{f}}^{\left(1\right)}}_{i+1/2}+{\omega_{2}}{{\hat{f}}^{\left(2\right)}}_{i+1/2},$
(7)

in which ${{\hat{f}}^{\left(k\right)}}_{i+1/2}\;\left({k=0,1,2}\right)$ are
the flux values assigned to three stencils
$\left\\{{{x_{j-2}},{x_{j-1}},{x_{j}}}\right\\}$,
$\left\\{{{x_{j-1}},{x_{j}},{x_{j+1}}}\right\\}$ and
$\left\\{{{x_{j}},{x_{j+1}},{x_{j+2}}}\right\\}$, and they are given by

$\begin{split}&{{\hat{f}}^{\left(0\right)}}_{i+1/2}=\frac{1}{6}\left({2{f_{j-2}}-7{f_{j-1}}+11{f_{j}}}\right),\\\
&{{\hat{f}}^{\left(1\right)}}_{i+1/2}=\frac{1}{6}\left({-{f_{j-1}}+5{f_{j}}+2{f_{j+1}}}\right),\\\
&{{\hat{f}}^{\left(2\right)}}_{i+1/2}=\frac{1}{6}\left({2{f_{j}}+5{f_{j+1}}-{f_{j+2}}}\right).\\\
\end{split}$ (8)

In addition, the nonlinear weights ${\omega_{0}}$, ${\omega_{1}}$ and
${\omega_{2}}$ are calculated as follow [19]

${\omega_{k}}=\frac{{{\alpha_{k}}}}{{{\alpha_{0}}+{\alpha_{1}}+{\alpha_{2}}}},\;{\alpha_{k}}={d_{k}}\left[{1+{{\left({\frac{{{\tau_{5}}}}{{{\beta_{k}}+\varepsilon}}}\right)}^{2}}}\right],k=0,1,2,$
(9)

where $\varepsilon$ is a small positive sensitivity parameter, which is used
to avoid divisions by zero in the weights formulation, and is chosen to be
$10^{-40}$ in this work as suggested in [19]. Coefficient $d_{k}$ is the
optimal weights given by

${d_{0}}=\frac{1}{{10}},{d_{1}}=\frac{6}{{10}},{d_{2}}=\frac{3}{{10}}.$ (10)

$\tau$ a global smoothness indicator used in WENO-Z scheme [19], which is
characterized by

${\tau_{5}}=\left|{{\beta_{0}}-{\beta_{2}}}\right|,$ (11)

In addition, ${{\beta_{k}}}$ is the smoothness indicator, which measures the
smoothness of the solution over a particular stencil, and are defined as

$\begin{split}&{\beta_{0}}=\frac{{13}}{{12}}\left({{f_{j-2}}-2{f_{j-1}}+{f_{j}}}\right)+\frac{1}{4}{\left({{f_{j-2}}-4{f_{j-1}}+3{f_{j}}}\right)^{2}},\\\
&{\beta_{1}}=\frac{{13}}{{12}}\left({{f_{j-1}}-2{f_{j}}+{f_{j+1}}}\right)+\frac{1}{4}{\left({{f_{j+1}}-{f_{j-1}}}\right)^{2}},\\\
&{\beta_{2}}=\frac{{13}}{{12}}\left({{f_{j}}-2{f_{j+1}}+{f_{j+2}}}\right)+\frac{1}{4}{\left({3{f_{j}}-4{f_{j+1}}+{f_{j+2}}}\right)^{2}}.\\\
\end{split}$ (12)

Finally, in order to distinguish the non-smooth scales from the smooth regions
in the neural network, an adaptive discontinuity indicator is needed.
Following the work of Fu [22], the separation between discontinuities and
smooth regions can be achieved by

${\gamma_{k}}=\frac{1}{{{{\left({{\beta_{k}}+\delta}\right)}^{p}}}},\;k=0,1,2,3,$
(13)

where $\delta$ is equal to $10^{-4}$ and $10^{-3}$ for one-dimensional and
multi-dimensional problems, $q$ is a scale separation parameter, the value of
it is set to 6, as suggested in [22], ${\beta_{3}}$ is another smoothness
indicator over the stencil of
$\left\\{{{x_{j+1,}}{x_{j+2,}}{x_{j+3}}}\right\\}$, and it is given by

${\beta_{3}}=\frac{1}{3}\left({{f_{i+1}}\left({22{f_{i+1}}-73{f_{i+2}}+29{f_{i+3}}}\right)+{f_{i+2}}\left({61{f_{i+2}}-49{f_{i+3}}}\right)+10{f_{i+3}}{f_{i+3}}}\right).$
(14)

Then the discontinuity indicator is defined as

${\rm{Flag=}}\left\\{\begin{array}[]{l}0,\;{\chi_{k}}>{C_{T}},\;\forall
k\in\left\\{{0,1,2,3}\right\\}\\\ 1,\;{\rm{otherwise}}\end{array}\right.,$
(15)

in which $\chi_{k}$ the normalized smoothness indicator given by
${\chi_{k}}={{{\gamma_{k}}}\mathord{\left/{\vphantom{{{\gamma_{k}}}{\sum\limits_{k=0}^{3}{{\gamma_{k}}}}}}\right.\kern-1.2pt}{\sum\limits_{k=0}^{3}{{\gamma_{k}}}}}\;\left({k=0,\ldots,3}\right)\;$
and $C_{T}$ is a built-in parameter defined as ${C_{T}}=5\times{10^{-4}}$.

##  3 Computational results and discussion

Figure 2: Comparison of the viscous Burgers equation at $t=0.2$ (a) and
$t=1.0$ (b), in which $\nu=10^{-4}/\pi$, and the reference solution is
computed with WENO-Z scheme using 1000 mesh points.

In this section, the one-dimensional viscous and inviscid Burgers equations
with Dirichlet boundary conditions are chosen to validate the present hPINN,

$\begin{split}&\frac{{\partial u}}{{\partial
t}}+\frac{{\partial\left({\frac{{{u^{2}}}}{2}}\right)}}{{\partial
x}}=\nu\frac{{{\partial^{2}}u}}{{\partial{x^{2}}}},\;x\in\left[{-1,1}\right],\;t\in\left[{0,1}\right],\\\
&u\left({0,x}\right)=-\sin\left({\pi x}\right),\\\
&u\left({t,-1}\right)=u\left({t,1}\right)=0,\end{split}$ (16)

and the predicted results are compared with the original PINN. As for the
neural networks architecture, we employ 5 hidden-layers with 20 neuron in each
layer and the number of residual points are 300 in training process, and the
learning rate is set as $10^{-4}$. In both cases, the network training is
performed up to the loss function is smaller than $10^{-5}$.

Figure 3: Comparison of the inviscid Burgers equation (i.e., $\nu=0.0$) at
$t=0.2$ (a) and $t=1.0$ (b), in which the reference solution is computed with
WENO-Z scheme using 1000 mesh points.

Fig. 2 and Fig. 3 illustrate the distributions of the predict solutions for
both the viscous and inviscid Burgers equations, in which the reference
solutions are obtained by using fifth order WENO-Z scheme with 1000 mesh
points. As depicted in these figures, it can be clearly seen that in regions
where the solution is smooth both the PINN and hPINN work well, when the shock
wave appears in the solution, the original PINN fails to provide reasonable
approximations at the shock position, and this phenomenon is more distinct for
the inviscid Burgers equation (see Fig. 2 and Fig. 3). Different from the
original PINN, the results obtained with present hPINN converges to the
reference solution with a non-oscillatory profile even in the inviscid case,
which further suggests that the performance of present hPINN is better that
previous PINN.

To demonstrate the robustness of the present hPINN, the global relative errors
obtained at various Runge-Kutta stages $q$ and the time-steps $\Delta t$ for
present hPINN are also investigated, and the corresponding results are shown
in Table 1, in which the results obtained by using original PINN are not
incorporate since it fails to predict reasonable results. As shown in this
table, one can find that for a fixed Runge-Kutta stage $q$/time-step $\Delta
t$, the global relative error is always increased/decresed with increasing
time-step $\Delta t$/ Runge-Kutta stage $q$. Moreover, due to the distinct
feature of the neural network, the Runge-Kutta stage used here can be an
arbitrarily larger number, making it possible to solve the nonlinear PDEs with
a very large time step, which is largely different from the traditional
numerical method.

Table 1: Global relative errors obtained with hPINN for different Runge-Kutta stages $q$ and time-step $\Delta t$ at $t=0.6$. | ${\nu}=10^{-4}/\pi$ |  | $\nu=0.0$
---|---|---|---
q | $\Delta t=0.1$ | $\Delta t=0.3$ | $\Delta t=0.6$ |  | $\Delta t=0.1$ | $\Delta t=0.3$ | $\Delta t=0.6$
1 | $2.3438\times 10^{-3}$ | $3.4727\times 10^{-1}$ | $4.3328\times 10^{-1}$ |  | $4.4754\times 10^{-3}$ | $3.8972\times 10^{-2}$ | $5.9224\times 10^{-1}$
4 | $2.8555\times 10^{-3}$ | $9.2065\times 10^{-3}$ | $3.3564\times 10^{-1}$ |  | $2.0912\times 10^{-3}$ | $5.3934\times 10^{-3}$ | $1.2390\times 10^{-1}$
10 | $1.8478\times 10^{-3}$ | $2.6818\times 10^{-2}$ | $4.2223\times 10^{-2}$ |  | $3.0013\times 10^{-3}$ | $5.1584\times 10^{-3}$ | $9.8035\times 10^{-3}$
50 | $2.2558\times 10^{-3}$ | $3.5475\times 10^{-3}$ | $9.6065\times 10^{-2}$ |  | $2.3808\times 10^{-3}$ | $3.7173\times 10^{-3}$ | $1.1018\times 10^{-2}$

##  4 Conclusion

In this paper, by incorporating a discontinuity indicator into the neural
network, a hybrid physics-informed neural network (hPINN) is proposed, in
which the automatic differentiation technique is employ to compute all the
differential operators of the PDEs in smooth regions, while the classical
WENO-Z scheme is adopted to compute the convection term in the vicinity of
discontinuities. Then the viscous and inviscid Burgers are selected to
validated the present model. Based on the results, we observe that the present
hPINN is able to provide reasonable approximation even at the shock position,
and overall performance of the present hPINN is better than the original PINN.
Further, due to the distinct feature of the neural network, the Runge-Kutta
stage used here can be an arbitrarily larger number, making it possible to
solve the nonlinear PDEs with a very large time step, which is very hard or
even impossible in traditional numerical method. Finally, we would like to
point out that although the present work focuses on the one-dimensional
problems, it is easily to extent the present model to the multi-dimensional
problems, which will be considered in ongoing work.

## CRediT authorship contribution statement

Chunyue Lv: Software, Methodology, Investigation. Lei Wang: Conceptualization,
Methodology, Writing - review ${\rm{\&}}$ editing.

## Acknowledgements

This work is financially supported by the National Natural Science Foundation
of China (Grant No. 12002320) and the Fundamental Research Funds for the
Central Universities (Grant No. CUGGC05),.

## References

  * [1] L. Debnath, Nonlinear partial differential equations for scientists and engineers, Birkhauser, Boston, MA, 1997.
  * [2] S. Mazumder, Numerical methods for partial differential equations: finite difference and finite volume methods, first ed., Academic Press, New York, 2015.
  * [3] H.G. Roos, M. Stynes, L. Tobiska, Robust numerical methods for singularly perturbed differential equations: convection-diffusion-reaction and flow problems, Springer, Berlin, 2008.
  * [4] L. Wang, X.G. Yang, H.L. Wang, Z.H. Chai, Z.C. Wei, A modified regularized lattice Boltzmann model for convection-diffusion equation with a source term, Appl. Math. Lett. 112 (2021) 106766.
  * [5] M. Dissanayake, N. Phan-Thien, Neural-network-based approximations for solving partial differential equations, Commun. Numer. Methods Eng. 10 (3) (1994) 195¨C201.
  * [6] J. Sirignano, K. Spiliopoulos, DGM: a deep learning algorithm for solving partial differential equations, J. Comput. Phys. 375 (2018) 1339¨C1364.
  * [7] K. Wu, D. Xiu, Data-driven deep learning of partial differential equations in modal space, J. Comput. Phys. 408 (2020) 109307.
  * [8] M. Raissi, P. Perdikaris, G.E. Karniadakis, Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378 (2019) 686¨C707.
  * [9] G.E. Karniadakis, I.G. Kevrekidis, L. Lu, Physics-informed machine learning, Nat. Rev. Phys. 3 (6) (2021) 422-440.
  * [10] N. Zobeiry, K.D. Humfeld, A physics-informed machine learning approach for solving heat transfer equation in advanced manufacturing and engineering applications, Eng. Appl. Artif. Intell. 101 (2021) 104232.
  * [11] Y. Xu, H. Zhang, Y. Li, Solving Fokker-Planck equation using deep learning. Chaos 30 (1) (2020) 013133.
  * [12] G. Cybenko, Approximation by superpositions of a sigmoidal function, Math. Control. Signals Syst. 2 (4) (1989) 303-314.
  * [13] C.M. Dafermos, C.M. Dafermos, Hyperbolic conservation laws in continuum physics, Springer, Berlin, 2005.
  * [14] Z. Mao, A.D. Jagtap, G.E. Karniadakis, Physics-informed neural networks for high-speed flows, Comput. Meth. Appl. Mech. Eng. 360 (2020) 112789.
  * [15] O. Fuks, H.A. Tchelepi, Limitations of physics informed machine learning for nonlinear two-phase transport in porous media, J. Mach. Learn. Model. Comput. 1 (1) (2020).
  * [16] M.M. Almajid, M.O. Abu-Al-Saud, Prediction of porous media fluid flow using physics informed neural networks, J. Pet. Sci. Eng. 208 (2022) 109205.
  * [17] R.G. Patel, I. Manickam, N.A. Trask, Thermodynamically consistent physics-informed neural networks for hyperbolic systems, J. Comput. Phys. (2021) 110754.
  * [18] A.G. Baydin, B.A. Pearlmutter, A.A. Radul, Automatic differentiation in machine learning: a survey, J. Mach. Learn. Res. 18 (2018).
  * [19] R. Borges, M. Carmona, B. Costa, An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws, J. Comput. Phys. 227 (6) (2008) 3191-3211.
  * [20] G.S. Jiang, C.W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1995) 202¨C228.
  * [21] X. Meng, Z. Li, D. Zhang, PPINN: Parareal physics-informed neural network for time-dependent PDEs, Comput. Meth. Appl. Mech. Eng. 370 (2020) 113250.
  * [22] L. Fu, A hybrid method with TENO based discontinuity indicator for hyperbolic conservation laws, Commun. Comput. Phys. 26 (4) (2019) 973-1007.


Finally, we conclude the section by considering the multiple key-intervention unit case as in Section <ref>.
The Hájek estimator of $\tau^\pi$ is given by

\begin{align*}
\hat{\tau}^\pi_{\text{H\'{a}jek}} =   \frac{\sum_{k=1}^{K}\frac{1}{|\mathcal{S}_k|}\sum_{j \in \mathcal{S}_k}\frac{\mathbbm{1}(\bm{A}^\top_{k\bm{i}^*}\bm{1} = |\bm{i}^*|p^*)}{f_k(\bm{A}^\top_{k\bm{i}^*}\bm{1} = |\bm{i}^*|p^*)}Y^{\text{obs}}_{kj} } {\sum_{k=1}^{K}\frac{1}{|\mathcal{S}_k|}\sum_{j \in \mathcal{S}_k}\frac{\mathbbm{1}(\bm{A}^\top_{k\bm{i}^*}\bm{1} = |\bm{i}^*|p^*)}{f_k(\bm{A}^\top_{k\bm{i}^*}\bm{1} = |\bm{i}^*|p^*)}} = \frac{\hat{\tau}^\pi_{\text{HT}}}{\hat{\lambda}^\pi_{\text{HT}}}.
\end{align*}

The form of the variance of $\hat{\tau}^\pi_{\text{H\'{a}jek}}$ and its estimator is analogous to those in Theorem <ref>, where $\hat{\mu}^\pi_{a,\text{HT}}$ is replaced by  $\hat{\tau}^\pi_{\text{HT}}$ and $\hat{\lambda}^\pi_{a,\text{HT}}$ is replaced by  $\hat{\lambda}^\pi_{\text{HT}}$, and the derivation is analogous to the proof of Theorem <ref>.

§.§ Additional results from the simulation study

8cModel-1     8cModel-2

2*$m_k$     2cBias     2cRMSE     2cCoverage     2cCI length     2cBias     2cRMSE     2cCoverage     2cCI length

HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek

4*$K = 10$      50     0.01     -0.00     0.21     0.13     0.95     0.94     0.86     0.50     0.01     0.00     0.31     0.25     0.95     0.95     1.23     0.98

100     -0.01     -0.01     0.18     0.12     0.95     0.93     0.70     0.47     -0.02     -0.02     0.27     0.25     0.94     0.94     1.05     0.96

250     0.00     0.00     0.13     0.10     0.95     0.95     0.51     0.42     0.00     0.00     0.22     0.21     0.94     0.95     0.86     0.81

500     -0.00     -0.00     0.13     0.11     0.94     0.93     0.47     0.40     0.00     0.00     0.20     0.19     0.95     0.96     0.79     0.76

4*$K = 50$       50     0.00     -0.00     0.09     0.06     0.95     0.94     0.37     0.23     -0.00     -0.00     0.13     0.11     0.96     0.95     0.52     0.44

100     0.00     0.00     0.07     0.05     0.95     0.95     0.29     0.21     0.00     0.00     0.12     0.11     0.94     0.94     0.47     0.43

250     0.00     0.00     0.06     0.05     0.96     0.95     0.23     0.19     0.00     0.00     0.10     0.09     0.95     0.95     0.40     0.37

500     -0.00     -0.00     0.06     0.05     0.94     0.95     0.21     0.19     -0.01     -0.01     0.10     0.09     0.94     0.95     0.38     0.37

Bias, RMSE, coverage and length of $95\%$ confidence intervals for $\mu^\pi_1$ under stochastic intervention $\pi^{(1)}_k(\cdot)$.

8cModel-1     8cModel-2

2*$m_k$     2cBias     2cRMSE     2cCoverage     2cCI length     2cBias     2cRMSE     2cCoverage     2cCI length

HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek

4*$K = 10$      50     0.03     0.00     0.49     0.26     0.96     0.97     2.04     1.15     0.03     0.00     0.56     0.35     0.96     0.98     2.39     1.65

100     -0.02     -0.01     0.41     0.25     0.95     0.97     1.64     1.07     -0.04     -0.03     0.45     0.34     0.97     0.99     1.98     1.64

250     0.00     0.00     0.28     0.21     0.96     0.98     1.13     0.96     0.00     0.00     0.34     0.29     0.97     0.99     1.50     1.38

500     -0.00     -0.00     0.27     0.21     0.94     0.97     1.01     0.93     -0.00     -0.00     0.31     0.27     0.97     0.98     1.35     1.31

4*$K = 50$      50     0.00     -0.00     0.22     0.12     0.96     0.97     0.89     0.53     -0.00     -0.01     0.24     0.16     0.98     0.99     1.03     0.76

100     0.00     0.00     0.16     0.11     0.95     0.97     0.67     0.47     0.00     0.00     0.20     0.15     0.96     0.97     0.84     0.71

250     0.00     0.00     0.13     0.10     0.97     0.97     0.52     0.44     0.00     0.01     0.15     0.13     0.98     0.98     0.69     0.62

500     -0.00     -0.00     0.12     0.10     0.94     0.97     0.45     0.44     -0.01     -0.01     0.14     0.13     0.97     0.99     0.63     0.63

Bias, RMSE, coverage and length of $95\%$ confidence intervals for $\text{DE}^\pi$ under stochastic intervention $\pi^{(1)}_k(\cdot)$.

8cModel-1     8cModel-2

2*$m_k$     2cBias     2cRMSE     2cCoverage     2cCI length     2cBias     2cRMSE     2cCoverage     2cCI length

HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek

4*$K = 10$      50     0.09     -0.00     1.99     0.29     0.94     0.96     7.67     1.05     0.09     -0.01     2.03     0.62     0.95     0.98     7.87     2.28

100     0.00     0.01     1.96     0.20     0.96     0.94     7.63     0.91     -0.02     -0.01     2.01     0.62     0.96     0.99     7.84     2.31

250     0.03     -0.01     1.97     0.23     0.95     0.96     7.78     0.90     0.02     -0.04     2.04     0.69     0.96     0.99     8.07     2.47

500     -0.01     -0.00     2.06     0.18     0.95     0.94     7.83     0.98     -0.00     0.00     2.07     0.66     0.96     0.99     7.88     2.43

4*$K = 50$      50     0.03     0.01     0.86     0.12     0.95     0.95     3.39     0.46     0.03     0.01     0.88     0.29     0.94     0.95     3.46     1.09

100     -0.01     -0.00     0.86     0.10     0.95     0.96     3.41     0.40     -0.02     -0.01     0.88     0.28     0.95     0.96     3.48     1.09

250     -0.00     0.00     0.87     0.08     0.95     0.94     3.38     0.32     -0.01     -0.00     0.90     0.22     0.95     0.95     3.50     0.86

500     0.00     0.00     0.85     0.07     0.96     0.95     3.43     0.27     0.00     0.01     0.89     0.21     0.96     0.94     3.57     0.78

Bias, RMSE, coverage and length of $95\%$ confidence intervals for $\mu^\pi_1$ under stochastic intervention $\pi^{(2)}_k(\cdot)$.

8cModel-1     8cModel-2

2*$m_k$     2cBias     2cRMSE     2cCoverage     2cCI length     2cBias     2cRMSE     2cCoverage     2cCI length

HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek     HT     Hájek

4*$K = 10$      50     0.02     -0.01     1.36     0.41     0.96     0.96     5.33     1.63     0.02     -0.03     1.47     0.70     0.96     0.98     5.79     2.78

100     -0.00     0.01     1.21     0.36     0.96     0.96     4.74     1.82     -0.03     -0.01     1.41     0.72     0.96     0.99     5.52     2.93

250     -0.03     -0.00     1.11     0.31     0.95     0.95     4.42     1.77     -0.04     -0.03     1.27     0.71     0.96     0.99     5.01     2.87

500     -0.01     0.00     1.12     0.28     0.94     0.95     4.26     2.29     -0.01     0.00     1.28     0.67     0.96     0.99     4.90     3.15

4*$K = 50$      50     -0.02     0.02     0.61     0.18     0.95     0.96     2.39     0.72     -0.02     0.02     0.67     0.33     0.95     0.95     2.66     1.29

100     0.01     -0.00     0.54     0.16     0.94     0.95     2.12     0.64     -0.00     -0.01     0.58     0.31     0.96     0.97     2.32     1.24

250     0.01     0.00     0.50     0.13     0.95     0.97     1.97     0.54     0.01     -0.00     0.54     0.25     0.94     0.97     2.11     1.03

500     -0.01     0.00     0.48     0.11     0.96     0.96     1.93     0.48     -0.01     0.01     0.51     0.24     0.96     0.96     2.08     0.95

Bias, RMSE, coverage and length of $95\%$ confidence intervals for $\text{DE}^\pi$ under stochastic intervention $\pi^{(2)}_k(\cdot)$.