# Ion acceleration at two collisionless shocks in a multicomponent plasma

Rajesh Kumar  Graduate School of Science, Osaka University, Japan  Youichi
Sakawa<EMAIL_ADDRESS>Institute of Laser Engineering, Osaka
University, Japan  Takayoshi Sano  Institute of Laser Engineering, Osaka
University, Japan  Leonard N. K. D$\ddot{\rm o}$hl  York Plasma Institute,
Department of Physics, University of York, United Kingdom  Nigel Woolsey  York
Plasma Institute, Department of Physics, University of York, United Kingdom
Alessio Morace  Institute of Laser Engineering, Osaka University, Japan

###### Abstract

Intense laser-plasma interactions are an essential tool for the laboratory
study of ion acceleration at a collisionless shock. With two-dimensional
particle-in-cell calculations of a multicomponent plasma we observe two
electrostatic collisionless shocks at two distinct longitudinal positions when
driven with a linearly-polarized laser at normalized laser vector potential
$a_{0}$ that exceeds 10. Moreover, these shocks, associated with protons and
carbon ions, show a power-law dependence on $a_{0}$ and accelerate ions to
different velocities in an expanding upstream with higher flux than in a
single-component hydrogen or carbon plasma. This results from an electrostatic
ion two-stream instability caused by differences in the charge-to-mass ratio
of different ions. Particle acceleration in collisionless shocks in
multicomponent plasma are ubiquitous in space and astrophysics, and these
calculations identify the possibility for studying these complex processes in
the laboratory.

##  I Introduction

Collisionless shocks under ambient magnetic field are ubiquitous in space and
astrophysical plasmas, and are believed to be sources for high-energy
particles or cosmic-rays Sagdeev and Shapiro (1973); Bell (1978); Blandford
and Ostriker (1978); Wu (1984); Ball and Melrose (2001); Hoshino (2001);
Sakawa  _et al._ (2016). Multiple collisionless shocks occur in plasmas
associated with planetary systems Sauer (1996); Shimazu (2001); Mazelle  _et
al._ (2004), where multicomponent plasmas occur as planetary material mixes
with the solar wind. In the magnetospheres of planets, such as Mars and Venus,
see Bertucci  _et al._ (2011), multicomponent plasmas occur and ions of
differing charge-to-mass ratio likely play a role. Jarvinen  _et al._ (2018)
discuss the role of oxygen in an induced Martian magnetosphere, where oxygen
is likely introduced by the past solar wind bombardment of water on the
unmagnetized surface of Mars. Multiple-reflection of solar-wind protons at the
Martian bow-shock was recorded across a shock by Mars Express and described by
Yamauchi  _et al._ (2012). These observations and the Voyager missions, see
for example Gurnett  _et al._ (2013), show multiple collisionless shocks are
associated with planetary and stellar systems. Borisov and Fraenz (2016)
illustrates this for Mars and Venus where the formation of a second
collisionless shock, in a region of magnetic pile-up between the bow shock and
ionosphere Bertucci  _et al._ (2005), results from the presence of planetary
oxygen ions and solar wind protons.

Collisionless shocks occur in much more extreme astrophysical systems Warren
_et al._ (2005); Caprioli  _et al._ (2017); Metzger (2020) such as supernova
remnants where a reverse shock, an inward-propagating collisionless shock,
heats stellar ejecta material containing a mixture of protons and heavy ions
Warren  _et al._ (2005). Warren  _et al._ (2005) observe localized regions
where strong line emission of Fe and Si ions occur in the reverse-shock heated
ejecta. Yamaguchi  _et al._ (2014) illustrate collisionless electron heating
at the front of the reverse shock caused by a cross-shock potential created by
charge deflection. Understanding of collisionless shocks and the associated
particle acceleration processes in multicomponent plasmas is of general
importance in space, astrophysics, and plasma physics.

While multiple collisionless shocks are expected in such systems, it is not
possible to observe them because of the limited resolution of the remote
sensing. It is possible that future spatially-resolved measurements using
multi-point spacecraft clusters might observe double-shock structures. Cohen
_et al._ (2019) and Broll  _et al._ (2018) demonstrate in situ spatially
resolved proton reflections Cohen  _et al._ (2019) and multi-ion (solar wind
protons and He2+ contamination) reflections Broll  _et al._ (2018) from a
shock in the Earth’s magnetosphere with the magnetospheric multiscale (MMS)
cluster. Laboratory experiments are a unique way of obtaining spatially
resolved measurements of collisionless shocks. They can provide tests of
understanding of particle acceleration in multiple collisionless shocks.
Numerical simulations by Schaeffer  _et al._ (2020) demonstrate the formation
of two collisionless shocks as a laser-ablated plasma acts as a piston pushes
on a magnetized multicomponent CH plasma. Laboratory studies show how ion
separation in unmagnetized multicomponent plasma is a common occurrence Byvank
_et al._ (2020); Rambo and Denavit (1994); Bellei  _et al._ (2014); Sio  _et
al._ (2019); Rinderknecht  _et al._ (2018). As examples, Byvank  _et al._
(2020) use merging plasma jets at oblique angles to observe ion and shock-
front separation when using jets that contain a mixture of He and Ar.
Rinderknecht  _et al._ (2018) observe ion velocity separation in a laser-
driven collisional shock generated in a multicomponent plasma, and ion-species
separation is predicted in inertial confinement fusion experiments as a strong
shock enters the fuel containing multiple-ion species Bellei  _et al._ (2014);
Sio  _et al._ (2019).

Continuing advances in high-intensity laser technology Danson  _et al._ (2019)
drives the development of compact, high-flux sources of energetic ionsDaido
_et al._ (2012); Macchi  _et al._ (2013). These sources may prove useful for
many applications Bulanov  _et al._ (2014); Li  _et al._ (2006); Roth  _et
al._ (2001). Among the many ion acceleration mechanisms being pursued Snavely
_et al._ (2000); Wilks  _et al._ (2001); Wagner  _et al._ (2016); Scott  _et
al._ (2018); Esirkepov  _et al._ (2004); Henig  _et al._ (2009a); Macchi  _et
al._ (2009); Kim  _et al._ (2016); Higginson  _et al._ (2018); YIN  _et al._
(2006); Henig  _et al._ (2009b); Stark  _et al._ (2019), collisionless shock
acceleration (CSA) of ions Denavit (1992); Silva  _et al._ (2004); Fiuza  _et
al._ (2012); Haberberger  _et al._ (2011); Tresca  _et al._ (2015); Zhang  _et
al._ (2015, 2017); Antici  _et al._ (2017); Chen  _et al._ (2017); Pak  _et
al._ (2018); Polz  _et al._ (2019); Ota  _et al._ (2019); Kumar  _et al._
(2019) is of particular relevance to space and astrophysical shocks. With CSA,
ions located ahead of an unmagnetized electrostatic (and collisionless) shock
Sakawa  _et al._ (2016) are reflected by the electrostatic potential of the
shock to twice the shock velocity Fiuza  _et al._ (2012).

Unmagnetized electrostatic collisionless shocks Forslund and Shonk (1970);
Forslund  _et al._ (1971) are rare among space Balogh, A and Treumann (2013)
and astrophysical systems, since shocks occur in collisionless magnetized
plasma. However, there are common and important collisionless processes
involved in both type of shocks Silva  _et al._ (2004). For example, particle
acceleration occurs in collisionless shocks Drury (1983); JONES (1991),
reflected particles excite two-stream instabilities Ohira and Takahara (2008);
Treumann (2009), reflected ions cause shock dissipation and reformation
Treumann (2009); Balogh, A and Treumann (2013); Madanian  _et al._ (2020),
effects of cross-shock electrostatic potential Bale and Mozer (2007); Cohen
_et al._ (2019), and so on. Bale et al. Bale  _et al._ (2005) describe shock
dissipation due to ion reflection in terms of the Cluster satellite mission.
Therefore, understanding of collisionless shocks and the associated particle
acceleration processes in multicomponent plasmas is of general importance in
space and astrophysical shocks. The study of collisionless shocks and particle
interaction is possible with laser-plasma systems. In this work the
colllisionless shock is mediated by an electrostatic interaction.

In laser-plasma experiments hydrogen and carbon are ever-present on the
surfaces of solid targets and inevitably result in multicomponent plasmas. A
number of studies Zhang  _et al._ (2017); Antici  _et al._ (2017); Pak  _et
al._ (2018); Ota  _et al._ (2019) specifically use multicomponent thin-foil
targets such as plastic (CH) or Mylar ($\rm C_{10}H_{8}O_{4}$), and in Kumar
_et al._ (2019) we reported on how target composition influences CSA by
comparing C2H3Cl, $\rm CH$, $\rm He^{3}H$, and H. Inclusion of a high atomic-
number element like Cl, results in partial ionization, to Cl15+, enabling the
study of a material with $\langle Z\rangle/\langle A\rangle<0.5$. An
electrostatic ion two-stream instability (EITI) excited in the multicomponent
plasma is central to the ion acceleration process with CSA accelerating
protons Kumar  _et al._ (2019) and heavier ions to the same velocity.

In comparison, the radiation pressure acceleration studies by Zhang  _et al._
(2009), which use circularly-polarized laser pulses and a three-layer
“sandwich” target containing protons and heavier-ions, show the emergence of
two shock fronts. One shock is associated with protons and the other with
heavier ions, the different species of ion are accelerated in different fields
to different velocities.

In this paper we examine, using the two-dimensional (2D) particle-in-cell
(PIC) simulation code EPOCH Arber  _et al._ (2015), the physical conditions
for the appearance of collisionless shocks and ion acceleration in a
multicomponent plasma formed from C2H3Cl and $\rm CH$ targets. We use a
linearly p-polarized laser pulse, and for a normalized vector potential
$a_{0}\geq 10$ to show the existence of two collisionless shock fronts. These
shocks are associated with the proton and C6+ ion populations. The shock front
accompanying the proton population propagates faster than the shock
accompanying the C6+ ions. As a result, CSA of protons and C6+ ions occurs at
different shocks and longitudinal locations in the plasma, producing ion
populations at different velocities.

##  II Particle-in-cell simulation

We study four values: $a_{0}=3.35$, $10$, $20$, and $33$, where $a_{0}=3.35$
corresponds to $1.4\times 10^{19}~{}\mathrm{W/cm^{2}}$ for the wavelength of
$1~{}\mathrm{\mu m}$. The simulated laser pulse uses a Gaussian temporal
profile with $1.5~{}\mathrm{ps}$ full-width-at-half-maximum. Figure 1 shows
the normalized initial electron density profile used in PIC simulations for
$a_{0}=3.35$. The simulated targets use a longitudinal ($x$-direction in Figs.
2 to 5) density profile consisting of an exponentially increasing
$5~{}\mathrm{\mu m}$ scale-length laser-irradiated front region,
$5~{}\mathrm{\mu m}$ uniform central region, and an exponentially decreasing
profile with $30~{}\mathrm{\mu m}$ scale-length rear region as the back of the
target. Details of the simulations including the target density profiles at
$a_{0}=3.35$ are given in Kumar  _et al._ (2019). When $a_{0}$ is varied, the
maximum electron density is increased to match the relativistic critical
density $a_{0}n_{cr}$, where $n_{cr}=1.12\times 10^{21}~{}\mathrm{cm^{-3}}$ is
the critical plasma density to the laser at $1~{}\mathrm{\mu m}$. The charge
states $Z_{i}$ of protons, C-ions, and Cl-ions are $1$, $6$, and $15$,
respectively. The corresponding ion density for each material is calculated
from the quasineutral plasma condition.

Figure 1: The normalized initial electron density profile used in PIC
simulations for $a_{0}=3.35$. The laser is from the left-hand side of the
simulation box. The density profile consists of an exponentially increasing
$5~{}\mathrm{\mu m}$ scale-length laser-irradiated front region, followed by
$5~{}\mathrm{\mu m}$ uniform central region, and an exponentially decreasing
rear-side profile with $30~{}\mathrm{\mu m}$ scale-length. To avoid boundary
effects, the simulations use 40 $\mu$m and 100 $\mu$m vacuum regions at the
front and rear of the target, respectively.

At $a_{0}=3.35$, as the relativistic electrons move through the plasma, the
inertia of the more massive ions sets up an electrostatic field, $E_{x}$. The
exponentially decreasing density profile on the rear side of the target
results in an electrostatic field or target-normal-sheath-acceleration field,
$E_{\rm TNSA}$ Kumar  _et al._ (2019). This TNSA field occurs in the upstream
region and results in the upstream ions moving at velocity $v_{0}$ in the
longitudinal direction Grismayer and Mora (2006).

Figure 2: Phase-space of (a) protons and (b) C6+ ions for $\rm C_{2}H_{3}Cl$
plasma at $a_{0}$ = 3.35 and at t = 4.0 ps. The horizontal lines represent the
lower threshold $v_{L}^{i}$ for ion reflection and the shock velocity in
proton density $V_{sh}^{\rm P}$. The color scale shows the number of ions in a
log scale.

To accelerate the ions via the CSA mechanism, the potential energy at the
collisionless shock must exceed the kinetic energy of the upstream expanding
ions. In other words, the electrostatic potential $\phi$ at the shock front
satisfies the following condition Tidman and Krall (1971):
$Z_{i}$$e\phi\geq\frac{1}{2}$$A_{i}$$m_{p}(V_{sh}^{i}-v_{0}^{i})^{2}$, where
$e$ is the electric charge, $A_{i}$ is the ion mass number, $m_{p}$ is the
proton mass, $V^{i}_{sh}$ is the shock velocity, and the superscript $i$
represents the different ion species. The lower ion-velocity threshold
($v_{L}^{i}$) for ion reflection and CSA is
$v_{L}^{i}=V_{sh}^{i}-\sqrt{2(Z_{i}/A_{i})}e\phi/m_{p}$ Kumar  _et al._
(2019). Therefore, CSA occurs for

$v_{L}^{i}\leq v_{0}^{i}\leq V_{sh}^{i}.$ (1)

Equation (1) represents the lower $v_{L}^{i}$ and upper $V_{sh}^{i}$ bounds in
$v_{0}^{i}$ for ion reflection. All ions with velocities $v_{0}^{i}$ between
$v_{L}^{i}$ and $V_{sh}^{i}$ are reflected at the collisionless shock and
leave with velocity 2$V_{sh}^{i}-v_{0}^{i}$ and the maximum velocity is
$V_{sh}^{i}+\sqrt{2(Z_{i}/A_{i})e\phi/m_{p}}$ = 2$V_{sh}^{i}-v_{L}^{i}$. For
protons $Z_{i}=A_{i}=1$, and the lower threshold is $v_{L}^{\rm P}=V_{sh}^{\rm
P}-\sqrt{2e\phi/m_{p}}$.

##  III Results

###  III.1 Double-shock formation

Figure 2 shows the phase-space of protons and C6+ ions at $a_{0}$ = 3.35. A
significant population of protons satisfy Eq. (1) and as result are
accelerated at the collisionless shock [Fig. 2(a)]. In comparison, relatively
few C6+ ions are reflected by the same collisionless shock, as this requires
$V_{sh}^{\rm C}$ = $V_{sh}^{\rm P}$ [Fig. 2(b)]. The lower threshold velocity
for carbon ions, $v_{L}^{\rm C}$ is slightly larger than $v_{L}^{\rm P}$ as
the charge-to-mass ratio, $Z_{i}/A_{i}$, is a factor of two smaller for C6+.
Furthermore, because of the smaller $Z_{i}/A_{i}$, the expansion velocity
$v_{0}^{\rm C}$ driven by $E_{\rm TNSA}$ in the upstream region is lower than
equivalent process for protons. This causes $v_{0}^{\rm C}$ to drop below
$v_{L}^{\rm C}$. This is illustrated in Fig. 2(b) which highlights how few C6+
ions are accelerated. Indeed, some of the energetic C6+ ions in Fig. 2(b)
likely originate early in time from the laser interaction at the front surface
of the plasma. We conclude that a negligible number of C6+ ions are
accelerated via the CSA mechanism for $a_{0}=3.35$.

Simulations at low-intensity ($a_{0}=3.35$) generate a single-shock. At higher
intensity, a key finding is the appearance of two distinct collisionless
shocks. Figure 3 shows results at $a_{0}=10$, and Fig. 3(a) illustrates the
longitudinal electrostatic field $E_{x}$, averaged over the y-axis, and
potential $\phi$ at $t=3.0~{}\mathrm{ps}$. Large amplitude changes in $E_{x}$
and $\phi$ are present at two different longitudinal positions. Large changes
in the normalized proton and C6+ ion densities are indicated, respectively, by
the dotted ($x\approx 112~{}\rm\mu m$) and solid ($x\approx 126~{}\rm\mu m$)
vertical lines in Fig. 3(b). Figure 3(c) shows how the normalized ion
populations have evolved 1 ps later at $t=4.0~{}\rm ps$. It is clear that the
position of the jump in proton and C6+ ion densities are different.

Figure 3: (a) The electrostatic field $E_{x}$ (left axis, blue line) and
potential $\phi$ (right axis, red line) at $t=3.0~{}\rm ps$. The normalized
proton $n_{H}/n_{cr}$ (left axis, blue line) and carbon $n_{\rm C}/n_{cr}$
(right axis, red line) densities at (b) $t=3.0~{}\rm ps$ and (c) $t=4.0~{}\rm
ps$ in $\rm C_{2}H_{3}Cl$ plasma for $a_{0}$ = 10. (a) is shown across a
narrow longitudinal range compared to (b) and (c).

Multiple shock structures are seen in the phase-space and velocity spectra in
the first and second columns, respectively, of Fig. 4. The three sets of data
are for single-component H, single-component C, and multicomponent $\rm
C_{2}H_{3}Cl$ plasmas at $t=4.0~{}\rm ps$. The positions of the shock fronts
highlighted in Figs. 4(c) and 4(e) are at the same longitudinal locations as
the jumps in $n_{\rm H}/n_{cr}$ and $n_{\rm C}/n_{cr}$ identified in Fig.
3(c). A large number of protons and some of the C6+ ions have velocities
greater than $v_{L}^{i}$ and so CSA increases the velocity of these ions to
$2V_{sh}^{i}-v_{L}^{i}$. In the $\rm C_{2}H_{3}Cl$ plasma, collisionless
shocks associated with the protons and separately with the C6+ ions accelerate
the protons and C6+ ions to different velocities. Figures 4(c) and 4(e)
indicate that the multicomponent $\rm C_{2}H_{3}Cl$ plasma develops, in the
expanding upstream, a broad velocity distribution within the proton and C6+
ion populations. This is driven by an electrostatic ion two-stream instability
(EITI) that arises from the velocity differences between the proton population
with $Z_{i}/A_{i}=1$, and the heavier C6+ ions with $Z_{i}/A_{i}=0.5$. We
refer to this as heavy-ion EITI or HI-EITI Kumar  _et al._ (2019). We find
that the HI-EITI decelerates some upstream protons, while it accelerates some
C6+ ions with velocities below $v_{L}^{\rm C}$ to velocities that exceed this
lower threshold, and thereby increases the population of C6+ ions available
for CSA Kumar  _et al._ (2019). Furthermore, the CSA reflected-ion population,
which moves at high velocity, causes an additional EITI with the slower moving
expanding plasma that forms the upstream. We refer to this as reflected-ion
EITI or RI-EITI Kumar  _et al._ (2019). Overall, RI-EITI accelerates the
slower upstream expanding ions towards higher velocity and promotes some ions,
both protons and C6+ ions, with velocities below $v_{L}^{i}$ to velocities
above the lower threshold. This further increases the ion population available
for CSA Kumar  _et al._ (2019). In a multicomponent $\rm C_{2}H_{3}Cl$ plasma,
both RI-EITI and HI-EITI play essential roles in enabling the acceleration of
C6+ ions.

Figure 4: The first and second columns indicate the ion phase-space and the
velocity spectra respectively for $a_{0}=10$ and at $t=4.0~{}\rm ps$. The
velocity spectra are taken in the upstream region immediately in front of the
shock across $\Delta x=3~{}\rm\mu m$. (a), (b) Results for protons from a
single-component H plasma. (c), (d) Results for protons from a $\rm
C_{2}H_{3}Cl$ plasma. (e), (f) Results for C6+ ions from a $\rm C_{2}H_{3}Cl$
plasma; and (g), (h) C6+ ions in single-component C plasma. The vertical lines
on phase-space in panels (a), (c), (e), and (g) identify the position of the
shock front. In panels (b), (d), (f) and (h), moving left to right, the dotted
lines indicate the positions of the lower threshold velocity ($v_{L}^{i}$),
shock velocity ($V_{sh}^{i}$), and the maximum velocity of the reflected ions
($2V_{sh}^{i}-v_{L}^{i}$). The color scale shows the number of ions in a log
scale.

Figure 5: The phase-space for (a) protons and (b) C6+ ions from a $\rm CH$
plasma for $a_{0}=10$ and at $t=4.0~{}\rm ps$. The vertical lines identify the
positions of the shock front associated with the protons (solid line) and the
C6+ ions (dashed line). The horizontal lines indicate the positions of the
lower threshold velocity ($v_{L}^{\rm P}$) and shock velocity ($V_{sh}^{\rm
P}$) of protons, and the maximum velocity of the reflected protons
($2V_{sh}^{\rm P}-v_{L}^{\rm P}$).

The color scale shows the number of ions in a log scale.

###  III.2 Double-step shock acceleration

In the previous Section, the multiple-shock (proton shock and C6+-ion shock)
formation is described in a multicomponent $\rm C_{2}H_{3}Cl$ plasma at
$a_{0}=10$, and protons and C6+ ions are reflected and accelerated by each
shock once. In this Section, we illustrate C6+-ion acceleration is a double-
step process with reflections at each shock in a multicomponent $\rm CH$
plasma at $a_{0}=10$.

Figures 5 show the phase-space for protons [Fig. 5(a)] and C6+ ions [Fig.
5(b)] in a $\rm CH$ plasma at $a_{0}=10$ and $t=4.0~{}\rm ps$. We see that in
this $\rm CH$ plasma, the high-mass C6+ ions are reflected and accelerated
twice; first at the C6+ ion-shock ($x\approx 143~{}\rm\mu m$) and second at
the proton-shock ($x\approx 177~{}\rm\mu m$) to the velocity $V/c=0.23$. This
is a clear observation of double-step multiple-shock acceleration of high-mass
C6+ ions in a multicomponent plasma. This double-step shock acceleration of
C6+ ions is clearly seen in Fig. 5 but not in a $\rm C_{2}H_{3}Cl$ plasma
[Fig. 4]. This is caused by a slightly faster proton-shock velocity of
$V_{sh}^{\rm P}/c=0.22$ in a $\rm CH$ plasma compared with $V_{sh}^{\rm
P}/c=0.20$ in a $\rm C_{2}H_{3}Cl$ plasma. As a result, in a $\rm CH$ plasma
$V_{sh}^{\rm P}$ is larger than the velocity of the pre-accelerated C6+ ions,
which are reflected and accelerated by the C6+-ion shock and likely originated
from the laser interaction at the front surface of the plasma early in time.
This results in the second acceleration of C6+ ions by the proton-shock. In
the case of a $\rm C_{2}H_{3}Cl$ plasma, $V_{sh}^{\rm P}$ is nearly equal to
the velocity of the pre-accelerated C6+ ions, and the second acceleration of
C6+ ions is not observed.

Furthermore, the respective deceleration and acceleration of expanding proton
and C6+ ion populations, as a result of HI-EITI, are more apparent in a $\rm
CH$ plasma compared with a $\rm C_{2}H_{3}Cl$ plasma.

###  III.3 The $a_{0}$ dependence of plasma parameters

Simulations show the formation of two collisionless shocks at $a_{0}\geq 10$,
and CSA of a significant number of C6+ ions in multicomponent plasmas. This is
qualitatively different from simulations at $a_{0}=3.35$ which show, see in
Fig. 2(b), a single shock. To understand the importance of increasing $a_{0}$,
we extend our numerical investigation of CSA to $a_{0}=20$ and $33$ in a $\rm
C_{2}H_{3}Cl$ plasma. These simulations confirm the existence of two
collisionless shocks and indicate that the Mach number depends on $a_{0}$.

In Fig. 6(a) we compare, at $t=4.0~{}\rm ps$, the upstream electron energy
distributions for different $a_{0}$ and fit these with two-dimensional-
relativistic (2D-relativistic) Maxwellian functions. The distributions at
$a_{0}=10$, $20$, and $33$ are described by a two-temperature fit representing
a bulk population and an energetic tail, while at $a_{0}$ = 3.35 the
distribution is described by a single temperature. The bulk and tail
Maxwellian components are shown for $a_{0}=33$. Figure 6(b) shows an $a_{0}$
power-law dependence for temperatures associated with the bulk and high-energy
parts of the electron distributions. The fitted electron temperatures do not
depend on the target material, as the laser intensity and the electron
densities are not material dependent but determined by $a_{0}$ Kumar  _et al._
(2019).

Figure 6: (a) The electron energy distribution taken at $t=4.0~{}\rm ps$ in
the upstream region of the shock front for different laser intensities
corresponding to $a_{0}=3.35$ ($\times$), $10$ (+), $20$ (-), and $33$ ($|$)
for $\rm C_{2}H_{3}Cl$. A sum of two (bulk and tail) 2D relativistic
Maxwellian is used to fit to the electron energy distribution shown by the
solid lines for $a_{0}=10$, $20$, and $33$. The bulk (dotted line) and tail
(dashed line) components for $a_{0}=33$ are shown. (b) The electron
temperatures as a function of $a_{0}$.

The shock velocity $V_{sh}^{i}$ and the mean velocity $v_{m}^{i}$ of the
expanding ions for all values of $a_{0}$ are higher in single-component H and
C plasmas, compared with a multicomponent $\rm C_{2}H_{3}Cl$ plasma.
Furthermore, the difference between $V_{sh}^{i}$ and $v_{m}^{i}$,
$v_{df}^{i}=V_{sh}^{i}-v_{m}^{i}$, increases with $a_{0}$ as a power-law
except at the highest intensity, where $a_{0}$ = 33, which results from
significant-levels ion reflection depleting or dissipating the collisionless
shock Liseykina  _et al._ (2015).

Figure 7: The $a_{0}$ dependence of (a) shock velocities $V_{sh}^{i}$, (b)
mean velocities of the expanding ions $v_{m}^{i}$, (c) ion-acoustic velocities
$c_{s}^{i}$, (d) difference between the shock velocity and mean velocity of
the expanding ions $v_{df}^{i}=V_{sh}^{i}-v_{m}^{i}$, and (e) the
corresponding Mach number $M^{i}=v_{df}^{i}/c_{s}^{i}$ at t = 4.0 ps for
protons in a single-component H plasma ($\bigcirc$), C6+ ions in a single-
component C plasma ($\triangle$), and protons (∙) and C6+ ions
($\blacktriangle$) in a multicomponent $\rm C_{2}H_{3}Cl$ plasma.

Figures 7(a) and 7(b) represent the shock velocity ($V^{i}_{sh}$) and the mean
velocity ($v^{i}_{m}$), respectively, of the expanding protons and C6+ ions as
a function of $a_{0}$ in a single-component H plasma, single-component C
plasma, and $\rm C_{2}H_{3}Cl$ plasmas. The proton and C6+-ion $V_{sh}^{i}$
and $v_{m}^{i}$ are always larger for the single-component H plasma and
single-component C plasma compared with the multicomponent $\rm C_{2}H_{3}Cl$
plasma, and follow the trend $V^{\rm P}_{\rm H}>V^{\rm P}_{\rm
C_{2}H_{3}Cl}>V^{\rm C}_{\rm C}>V^{\rm C}_{\rm C_{2}H_{3}Cl}$ for all laser
intensities. Here superscripts P and C denote protons and C6+ ions,
respectively, with different plasmas indicated by subscripts. The ordering of
velocities results from differences in the average charge-to-mass ratio
$\langle Z\rangle/\langle A\rangle$, that is as $V^{i}_{sh}$ and $v^{i}_{m}$
are predominantly determined by the ion-acoustic velocity ($c^{i}_{s}$) and
velocity of ions due to $E_{\rm TNSA}$, respectively. Differences in the hole-
boring velocity, which depends on $\sqrt{\langle Z\rangle/\langle A\rangle}$,
explains why $V^{\rm P}_{\rm H}>V^{\rm P}_{\rm C_{2}H_{3}Cl}$ and $V^{\rm
C}_{\rm C}>V^{\rm C}_{\rm C_{2}H_{3}Cl}$ Kumar  _et al._ (2019). As a result,
the shock velocity in a single-component H plasma (with $\langle
Z\rangle/\langle A\rangle$ = 1) is larger than that in $\rm C_{2}H_{3}Cl$
($\langle Z\rangle/\langle A\rangle$ = 0.48), and shock velocity for C6+ ion
in a single-component C plasma ($\langle Z\rangle/\langle A\rangle$ = 0.50) is
larger than that in $\rm C_{2}H_{3}Cl$.

Ion-acoustic waves are excited in proton and C6+ ion populations and using the
bulk electron temperatures $T_{e}$ to derive an ion-acoustic velocity,
$c_{s}^{i}=\sqrt{(Z_{i}/A_{i})T_{e}/m_{p}}$, we find that the associated Mach
numbers, $M^{i}=v_{df}^{i}/c_{s}^{i}$, scale as a power law in $a_{0}$. The
ion-acoustic velocities for protons ($c_{s}^{\rm P}$) and C6+ ions
($c_{s}^{\rm C}$) are calculated using the bulk temperature of the plasma. The
upstream bulk temperatures in a single-component H plasma, single-component C
plasma, and multicomponent $\rm C_{2}H_{3}Cl$ plasma are the same, as a result
the $c_{s}^{i}$ depends on the $\sqrt{\langle Z\rangle/\langle A\rangle}$.
This is shown in Fig. 7(c). The $c_{s}^{i}$, indicated by the solid lines,
scale with $a_{0}$ as a power-law. The difference between the shock velocity
and mean velocity of the expanding ions, i.e.,
$v_{df}^{i}=V_{sh}^{i}-v_{m}^{i}$, is shown in Fig. 7(d) and increases as a
power-law with $a_{0}$ except at $a_{0}$ = 33.

The ratio between the $v_{df}^{i}$ and $c_{s}^{i}$ yields the Mach number
$M^{i}=v_{df}^{i}/c_{s}^{i}$, this is shown in Fig. 7(e). In comparison with
$M^{\rm P}$, the Mach number for protons, $M^{\rm C}$, the Mach number for C6+
ions, has a strong scaling with $a_{0}$. Notice in multicomponent $\rm
C_{2}H_{3}Cl$, $M^{\rm C}<1$ for $a_{0}$ = 3.35, and no shock is associated
with the C6+ ions. Furthermore, $M^{\rm P}$ in a single-component H plasma
decreases with $a_{0}$, this occurs as $v_{m}^{\rm P}$ scales faster with
$a_{0}$ than $V_{sh}^{\rm P}$, causing a slow scaling of $v_{df}^{\rm P}$
compared with $c_{s}^{\rm P}$ as $a_{0}$ increases.

Figure 8: (a) The temporal evolution of shock positions $X_{sh}$ and shock
velocities $V_{sh}$ for protons in a single-component H plasma (shown in red)
and a multicomponent $\rm C_{2}H_{3}Cl$ plasma (shown in blue) for $a_{0}$ =
3.35. $X_{sh}$ data are shown as open circles (a single-component H
plasma:$\bigcirc$) and open triangles (a multicomponent $\rm C_{2}H_{3}Cl$
plasma:$\triangle$). The derivative of $X_{sh}$ with respect to time gives
$V_{sh}$. $X_{sh}$ (dotted lines) and $V_{sh}$ (solid lines) rise
exponentially with time. (b) The temporal evolution of $X_{sh}$ and $V_{sh}$
for protons (shown in red) and C6+ ions (shown in blue) in a multicomponent
$\rm C_{2}H_{3}Cl$ plasma for $a_{0}$ = 33. $X_{sh}$ data are shown as open
circles (protons:$\bigcirc$) and open triangles (C6+ ions:$\triangle$). The
time dependencies of $X_{sh}$ and $V_{sh}$ are best represented by a third-
order (dotted lines) and second-order (solid lines) order polynomials,
respectively.

###  III.4 Shock dissipation

For $a_{0}$ = 33 shock dissipation, driven by ion reflection, becomes more
pronounced. This reduces the shock velocity Liseykina  _et al._ (2015).
Evidence for this is seen in Fig. 7(d) of $v_{df}^{i}$ and in Fig. 7(e) of
$M^{i}$ which illustrate a power-law trend for $a_{0}=3.35$, $10$, and $20$ up
to the end of the simulations at $t=4.0$ ps. Simulations at $a_{0}=33$ show
significant shock dissipation from $t\approx 2.5$ ps. For $a_{0}=3.35$, $10$,
and $20$, shock velocities increase exponentially with time until $t=4.0~{}\rm
ps$, in contrast, for $a_{0}=33$, the shock velocity increases to $t<2.5~{}\rm
ps$ then dissipates, which results in low Mach numbers at $a_{0}=33$ for
single- and multicomponent plasmas Liseykina  _et al._ (2015).

The temporal variation of shock positions ($X_{sh}$) in a single-component H
plasma and a multicomponent $\rm C_{2}H_{3}Cl$ plasma for $a_{0}$ = 3.35 are
shown in Fig. 8(a). The derivative $dX_{sh}/dt$ gives the shock velocity
$V_{sh}$. Since there is an exponential drop in the density at the rear-side
of the target, $V_{sh}$ increases exponentially as a function of time for all
target materials. For $a_{0}$ = 3.35 [see Fig. 8(a)], 10, and 20, $X_{sh}$ and
$V_{sh}$ rise exponentially with time. Comparing this to Fig. 8(b), we see
that at $a_{0}$ = 33 the temporal evolution of $X_{sh}$ and $V_{sh}$ for
protons and C6+ ions in a multicomponent $\rm C_{2}H_{3}Cl$ plasma is slower.
Indeed, the time dependencies of $X_{sh}$ and $V_{sh}$ are best represented by
third- and second-order polynomials, respectively. This slower temporal
evolution results from enhanced ion reflection at the shock which increases
shock dissipation Liseykina  _et al._ (2015).

Figure 9: The $a_{0}$ dependence of (a) energy $E$ per nucleon and (b) the
number $dN/dE$ at $E$ of reflected protons and C6+ ions at the peak of the
energy distribution at t = 4.0 ps for protons in single-component H
($\bigcirc$), C6+ ions in single-component C ($\triangle$), and protons (∙)
and C6+ ions ($\blacktriangle$) in $\rm C_{2}H_{3}Cl$ plasmas.

###  III.5 The $a_{0}$ dependence of ion acceleration

Figures 9(a) and 9(b) show how the energy $E$ and the number of reflected ions
at the peak of the energy distribution $dN/dE$ depend on $a_{0}$. In the $\rm
C_{2}H_{3}Cl$ plasma there are no C6+-associated shocks at $a_{0}=3.35$ as
$M^{\rm C}<1$. The energies of the reflected ions are always larger in single-
component H or C plasma when compared with multicomponent $\rm C_{2}H_{3}Cl$
plasma [Fig. 9(a)]. This is a feature of smaller $V_{sh}^{i}$ and amplitude of
$\phi$ in the multicomponent $\rm C_{2}H_{3}Cl$ compared with the single-
component H or C plasma as shown in Fig. 4.

Figures 10(a) and 10(b) show the spatial profile of electrostatic potentials
$\phi$ in a multicomponent $\rm C_{2}H_{3}Cl$ plasma and a single-component H
plasma, respectively, at $t=$ 2.5 (blue curve) and 4.0 ps (red curve) for
$a_{0}=3.35$. The vertical lines indicate the position of the shock fronts.
These highlight that $\phi$ is smaller in a multicomponent $\rm C_{2}H_{3}Cl$
plasma compared to a single-component H plasma. The smaller $\phi$ is a result
of a lower $\langle Z\rangle/\langle A\rangle$ plasma. In a single-component H
plasma, the gradient $d\phi/dx$ is large. This potential jump is associated
with the shock and is necessary for ion acceleration. It is produced by a
charge separation between electrons and ions. This feature is not observed in
the multicomponent $\rm C_{2}H_{3}Cl$ plasma because the charge separation
between electrons and ions is smeared out over a larger volume by the heavier
C and Cl ions, as a result the amplitude and gradient associated with $\phi$
are smaller.

Our PIC results indicate that $V_{sh}$ is smaller at lower $\langle
Z\rangle/\langle A\rangle$, this is explained by recognising that the hole-
boring velocity Robinson  _et al._ (2009)

$V_{HB}=c\sqrt{a_{0}^{2}\frac{\langle Z\rangle}{\langle
A\rangle}\frac{m_{e}}{m_{p}}\frac{n_{\rm cr}}{n_{e}}},$ (2)

determines the velocity of the piston driving the collisionless shock. Given
that $a_{0}$ and $n_{e}$ are same for all target materials, $V_{HB}$ has
relative dependence on only $\sqrt{\langle Z\rangle/\langle A\rangle}$,
maximizing $V_{sh}$ when $\langle Z\rangle/\langle A\rangle$ is largest, i.e.,
for a single-component H plasma.

Figure 10: The spatial profile of electrostatic potentials in (a) a
multicomponent $\rm C_{2}H_{3}Cl$ plasma and (b) a single-component hydrogen
plasma at $t=2.5$ (blue curves) and 4.0 ps (red curves) for $a_{0}=3.35$. The
vertical lines indicate the position of the shock fronts.

For the multicomponent $\rm C_{2}H_{3}Cl$ plasma at $a_{0}>3.35$ the flux of
the reflected protons and C6+ ions is higher [see Fig. 9(b)]. It is important
to note that more protons are accelerated in multicomponent plasma as
$v_{L}^{\rm P}$ is lower in comparison with the single-component H plasma. For
C6+ ion acceleration, HI-EITI, which is present only in a multicomponent
plasma, broadens the velocity distribution of the expanding C6+ ions towards
higher velocity. This results in more C6+ ions being available for CSA in
comparison to single-component C plasmas.

These results confirm our earlier observation Kumar  _et al._ (2019) that only
proton collisionless shocks were observed in multicomponent plasmas at $a_{0}$
= 3.35. In this work Mach numbers $M$ = 1.6 - 1.7 were calculated with a
critical Mach number needed for the proton reflection and CSA. These values
were derived using ion-acoustic velocities based on a $\langle
Z\rangle/\langle A\rangle$, where $\langle Z\rangle$ and $\langle A\rangle$
are the respective averages of $\langle Z_{i}\rangle$ and $\langle
A_{i}\rangle$ across all ion species in a plasma. Here, we use ion-specific
ion-acoustic velocities to describe the two collisionless shocks. For protons
$c_{s}^{\rm P}$ determines the Mach number $M^{\rm P}$ of a proton
collisionless shock, and ions satisfying the reflection condition given by Eq.
(1) are accelerated even when the Mach number is less than $M_{\rm cr}$
defined in Ref. Kumar  _et al._ (2019).

##  IV Summary

Two-dimensional PIC simulations are used to investigate the evolution of
electrostatic collisionless shocks and CSA of protons and heavy ions in
multicomponent plasmas. The interaction of a high-intensity p-polarized laser
with $\rm C_{2}H_{3}Cl$ and $\rm CH$ plasmas leads to the formation of the two
shock fronts in proton and C6+-ion populations. Both shocks have different
amplitudes of the shock potential and propagate with different velocities. The
electron temperature, shock velocities, and Mach numbers for shocks associated
with proton ($M^{\rm P}$) and C6+ ions ($M^{\rm C}$) scale as a power-laws
with the normalized laser intensity $a_{0}$. In the multicomponent $\rm
C_{2}H_{3}Cl$ plasma, $M^{\rm C}$ scales faster with $a_{0}$ compared to
$M^{\rm P}$. At $a_{0}=3.35$, as $M^{\rm C}<1$, a C6+ ion-shock does not form.
On increasing $a_{0}$, shock formation with CSA of protons and C6+ ions occurs
at different location and velocities. Double-step shock acceleration is
investigated in a $\rm CH$ plasma, in which the pre-accelerated C6+ ions are
further accelerated at the proton-shock. A broadening upwards of the C6+ ion
velocity distribution, as a result of a HI-EITI, is important and increases
the number of C6+ ions accelerated. For $a_{0}$ = 33 shock dissipation, driven
by ion reflection, becomes more pronounced. This results in the reduction of
the shock velocity. Moreover, modern ultra-intense, picosecond duration lasers
enable the laboratory study of the formation and modification of collisionless
shocks as ions are accelerated in multicomponent plasmas. These topics are
important to space physics, astrophysics, and plasma physics.

##  V Acknowledgement

This research was partially supported by Japan Society for the Promotion of
Science (JSPS) KAKENHI Grant No. JP15H02154, JP17H06202, JP19H00668, JSPS
Core-to-Core Program B. Asia-Africa Science Platforms Grant No.
JPJSCCB20190003, EPSRC grant EP/L01663X/1 and EP/P026796/1.

## References

  * Sagdeev and Shapiro (1973) R. Sagdeev and V. Shapiro, JETP Lett. (USSR) (Engl. Transl.), v. 17, no. 7, pp. 279-282 (1973).
  * Bell (1978) A. R. Bell, Monthly Notices of the Royal Astronomical Society 182, 147 (1978).
  * Blandford and Ostriker (1978) R. D. Blandford and J. P. Ostriker, The Astrophysical Journal 221, L29 (1978).
  * Wu (1984) C. S. Wu, Journal of Geophysical Research 89, 8857 (1984).
  * Ball and Melrose (2001) L. Ball and D. B. Melrose, Publications of the Astronomical Society of Australia 18, 361 (2001).
  * Hoshino (2001) M. Hoshino, Progress of Theoretical Physics Supplement 143, 149 (2001).
  * Sakawa  _et al._ (2016) Y. Sakawa, T. Morita, Y. Kuramitsu, and H. Takabe, Advances in Physics: X 1, 425 (2016).
  * Sauer (1996) K. Sauer, Geophysical Research Letters 23, 3643 (1996).
  * Shimazu (2001) H. Shimazu, JOURNAL OF GEOPHYSICAL RESEARCH 106, 8333 (2001).
  * Mazelle  _et al._ (2004) C. Mazelle, D. Winterhalter, and K. Sauer, Space Science Reviews 111, 115 (2004).
  * Bertucci  _et al._ (2011) C. Bertucci, F. Duru, N. Edberg, M. Fraenz, C. Martinecz, K. Szego, and O. Vaisberg,  _Space Science Reviews_, Vol. 162 (2011) pp. 113–171.
  * Jarvinen  _et al._ (2018) R. Jarvinen, D. A. Brain, R. Modolo, A. Fedorov, and M. Holmström, Journal of Geophysical Research: Space Physics 123, 1678 (2018).
  * Yamauchi  _et al._ (2012) M. Yamauchi, Y. Futaana, A. Fedorov, R. A. Frahm, E. Dubinin, R. Lundin, J. A. Sauvaud, J. D. Winningham, S. Barabash, and M. Holmström, Earth, Planets and Space 64, 61 (2012).
  * Gurnett  _et al._ (2013) D. A. Gurnett, W. S. Kurth, L. F. Burlaga, and N. F. Ness, Science 341, 1489 (2013).
  * Borisov and Fraenz (2016) N. Borisov and M. Fraenz, Physics of Plasmas 23, 122109 (2016).
  * Bertucci  _et al._ (2005) C. Bertucci, C. Mazelle, M. H. Acuña, C. T. Russell, and J. A. Slavin, Journal of Geophysical Research: Space Physics 110, A01209 (2005).
  * Warren  _et al._ (2005) J. S. Warren, J. P. Hughes, C. Badenes, P. Ghavamian, C. F. McKee, D. Moffett, P. P. Plucinsky, C. Rakowski, E. Reynoso, and P. Slane, The Astrophysical Journal 634, 376 (2005).
  * Caprioli  _et al._ (2017) D. Caprioli, D. T. Yi, and A. Spitkovsky, Physical Review Letters 119, 171101 (2017), arXiv:1704.08252 .
  * Metzger (2020) B. D. Metzger, Living Reviews in Relativity 23, 1 (2020), arXiv:1910.01617 .
  * Yamaguchi  _et al._ (2014) H. Yamaguchi, K. A. Eriksen, C. Badenes, J. P. Hughes, N. S. Brickhouse, A. R. Foster, D. J. Patnaude, R. Petre, P. O. Slane, and R. K. Smith, The Astrophysical Journal 780, 136 (2014).
  * Cohen  _et al._ (2019) I. J. Cohen, S. J. Schwartz, K. A. Goodrich, N. Ahmadi, R. E. Ergun, S. A. Fuselier, M. I. Desai, E. R. Christian, D. J. McComas, G. P. Zank, J. R. Shuster, S. K. Vines, B. H. Mauk, R. B. Decker, B. J. Anderson, J. H. Westlake, O. Le Contel, H. Breuillard, B. L. Giles, R. B. Torbert, and J. L. Burch, Journal of Geophysical Research: Space Physics 124, 3961 (2019).
  * Broll  _et al._ (2018) J. M. Broll, S. A. Fuselier, K. J. Trattner, S. J. Schwartz, J. L. Burch, B. L. Giles, and B. J. Anderson, Geophysical Research Letters 45, 49 (2018).
  * Schaeffer  _et al._ (2020) D. B. Schaeffer, W. Fox, J. Matteucci, K. V. Lezhnin, A. Bhattacharjee, and K. Germaschewski, Physics of Plasmas 27, 042901 (2020).
  * Byvank  _et al._ (2020) T. Byvank, S. J. Langendorf, C. Thoma, and S. C. Hsu, Physics of Plasmas 27, 042302 (2020).
  * Rambo and Denavit (1994) P. W. Rambo and J. Denavit, Physics of Plasmas 1, 4050 (1994).
  * Bellei  _et al._ (2014) C. Bellei, H. Rinderknecht, A. Zylstra, M. Rosenberg, H. Sio, C. K. Li, R. Petrasso, S. C. Wilks, and P. A. Amendt, Physics of Plasmas 21, 056310 (2014).
  * Sio  _et al._ (2019) H. Sio, J. A. Frenje, A. Le, S. Atzeni, T. J. T. Kwan, M. Gatu Johnson, G. Kagan, C. Stoeckl, C. K. Li, C. E. Parker, C. J. Forrest, V. Glebov, N. V. Kabadi, A. Bose, H. G. Rinderknecht, P. Amendt, D. T. Casey, R. Mancini, W. T. Taitano, B. Keenan, A. N. Simakov, L. Chacón, S. P. Regan, T. C. Sangster, E. M. Campbell, F. H. Seguin, and R. D. Petrasso, Physical Review Letters 122, 035001 (2019).
  * Rinderknecht  _et al._ (2018) H. G. Rinderknecht, H. S. Park, J. S. Ross, P. A. Amendt, S. C. Wilks, J. Katz, N. M. Hoffman, G. Kagan, E. L. Vold, B. D. Keenan, A. N. Simakov, and L. Chacón, Physics of Plasmas 25, 056312 (2018).
  * Danson  _et al._ (2019) C. N. Danson, C. Haefner, J. Bromage, and T. Butcher, High Power Laser Science and Engineering 7, e54 (2019).
  * Daido  _et al._ (2012) H. Daido, M. Nishiuchi, and A. S. Pirozhkov, Reports on progress in physics. Physical Society (Great Britain) 75, 056401 (2012).
  * Macchi  _et al._ (2013) A. Macchi, M. Borghesi, and M. Passoni, Reviews of Modern Physics 85, 751 (2013).
  * Bulanov  _et al._ (2014) S. V. Bulanov, J. J. Wilkens, T. Z. Esirkepov, G. Korn, G. Kraft, S. D. Kraft, M. Molls, and V. Khoroshkov, Physics-Uspekhi 57, 1149 (2014).
  * Li  _et al._ (2006) C. K. Li, F. H. Séguin, J. A. Frenje, J. R. Rygg, R. D. Petrasso, R. P. J. Town, P. A. Amendt, S. P. Hatchett, O. L. Landen, A. J. Mackinnon, P. K. Patel, V. A. Smalyuk, T. C. Sangster, and J. P. Knauer, Physical Review Letters 97, 135003 (2006).
  * Roth  _et al._ (2001) M. Roth, T. E. Cowan, M. H. Key, S. P. Hatchett, C. Brown, W. Fountain, J. Johnson, D. M. Pennington, R. A. Snavely, S. C. Wilks, K. Yasuike, H. Ruhl, P. Pegoraro, S. V. Bulanov, E. M. Campbell, M. D. Perry, and H. Powell, Physical review letters 86, 436 (2001).
  * Snavely  _et al._ (2000) R. A. Snavely, M. H. Key, S. P. Hatchett, T. E. Cowan, M. Roth, T. W. Phillips, M. A. Stoyer, E. A. Henry, T. C. Sangster, M. S. Singh, S. C. Wilks, A. MacKinnon, A. Offenberger, D. M. Pennington, K. Yasuike, A. B. Langdon, B. F. Lasinski, J. Johnson, M. D. Perry, and E. M. Campbell, Physical Review Letters 85, 2945 (2000).
  * Wilks  _et al._ (2001) S. C. Wilks, A. B. Langdon, T. E. Cowan, M. Roth, M. Singh, S. Hatchett, M. H. Key, D. Pennington, A. MacKinnon, and R. A. Snavely, Physics of Plasmas 8, 542 (2001).
  * Wagner  _et al._ (2016) F. Wagner, O. Deppert, C. Brabetz, P. Fiala, A. Kleinschmidt, P. Poth, V. A. Schanz, A. Tebartz, B. Zielbauer, M. Roth, T. Stöhlker, and V. Bagnoud, Physical Review Letters 116, 205002 (2016).
  * Scott  _et al._ (2018) G. G. Scott, D. C. Carroll, S. Astbury, R. J. Clarke, C. Hernandez-Gomez, M. King, A. Alejo, I. Y. Arteaga, R. J. Dance, A. Higginson, S. Hook, G. Liao, H. Liu, S. R. Mirfayzi, D. R. Rusby, M. P. Selwood, C. Spindloe, M. K. Tolley, F. Wagner, E. Zemaityte, M. Borghesi, S. Kar, Y. Li, M. Roth, P. McKenna, and D. Neely, Physical Review Letters 120, 204801 (2018).
  * Esirkepov  _et al._ (2004) T. Esirkepov, M. Borghesi, S. V. Bulanov, G. Mourou, and T. Tajima, Physical Review Letters 92, 175003 (2004).
  * Henig  _et al._ (2009a) A. Henig, S. Steinke, M. Schnürer, T. Sokollik, R. Hörlein, D. Kiefer, D. Jung, J. Schreiber, B. M. Hegelich, X. Q. Yan, J. Meyer-ter-Vehn, T. Tajima, P. V. Nickles, W. Sandner, and D. Habs, Physical Review Letters 103, 245003 (2009a), arXiv:0908.4057 .
  * Macchi  _et al._ (2009) A. Macchi, S. Veghini, and F. Pegoraro, Physical Review Letters 103, 085003 (2009).
  * Kim  _et al._ (2016) I. J. Kim, K. H. Pae, I. W. Choi, C.-L. Lee, H. T. Kim, H. Singhal, J. H. Sung, S. K. Lee, H. W. Lee, P. V. Nickles, T. M. Jeong, C. M. Kim, and C. H. Nam, Physics of Plasmas 23, 070701 (2016).
  * Higginson  _et al._ (2018) A. Higginson, R. J. Gray, M. King, R. J. Dance, S. D. R. Williamson, N. M. H. Butler, R. Wilson, R. Capdessus, C. Armstrong, J. S. Green, S. J. Hawkes, P. Martin, W. Q. Wei, S. R. Mirfayzi, X. H. Yuan, S. Kar, M. Borghesi, R. J. Clarke, D. Neely, and P. McKenna, Nature Communications 9, 724 (2018).
  * YIN  _et al._ (2006) L. YIN, B. J. ALBRIGHT, B. M. HEGELICH, and J. C. FERNÁNDEZ, Laser and Particle Beams 24, 291 (2006).
  * Henig  _et al._ (2009b) A. Henig, D. Kiefer, K. Markey, D. C. Gautier, K. A. Flippo, S. Letzring, R. P. Johnson, T. Shimada, L. Yin, B. J. Albright, K. J. Bowers, J. C. Fernández, S. G. Rykovanov, H.-C. Wu, M. Zepf, D. Jung, V. K. Liechtenstein, J. Schreiber, D. Habs, and B. M. Hegelich, Physical Review Letters 103, 045002 (2009b).
  * Stark  _et al._ (2019) D. J. Stark, L. Yin, and B. J. Albright, Physics of Plasmas 26 123101 (2019), 10.1063/1.5121430.
  * Denavit (1992) J. Denavit, Physical Review Letters 69, 3052 (1992).
  * Silva  _et al._ (2004) L. O. Silva, M. Marti, J. R. Davies, R. A. Fonseca, C. Ren, F. Tsung, and W. B. Mori, Physical Review Letters 92, 015002 (2004).
  * Fiuza  _et al._ (2012) F. Fiuza, A. Stockem, E. Boella, R. A. Fonseca, L. O. Silva, D. Haberberger, S. Tochitsky, C. Gong, W. B. Mori, and C. Joshi, Physical Review Letters 109, 215001 (2012).
  * Haberberger  _et al._ (2011) D. Haberberger, S. Tochitsky, F. Fiuza, C. Gong, R. A. Fonseca, L. O. Silva, W. B. Mori, and C. Joshi, Nature Physics 8, 95 (2011).
  * Tresca  _et al._ (2015) O. Tresca, N. P. Dover, N. Cook, C. Maharjan, M. N. Polyanskiy, Z. Najmudin, P. Shkolnikov, and I. Pogorelsky, Physical Review Letters 115, 094802 (2015).
  * Zhang  _et al._ (2015) H. Zhang, B. F. Shen, W. P. Wang, Y. Xu, Y. Q. Liu, X. Y. Liang, Y. X. Leng, R. X. Li, X. Q. Yan, J. E. Chen, and Z. Z. Xu, Physics of Plasmas 22, 013113 (2015).
  * Zhang  _et al._ (2017) H. Zhang, B. F. Shen, W. P. Wang, S. H. Zhai, S. S. Li, X. M. Lu, J. F. Li, R. J. Xu, X. L. Wang, X. Y. Liang, Y. X. Leng, R. X. Li, and Z. Z. Xu, Physical Review Letters 119, 164801 (2017).
  * Antici  _et al._ (2017) P. Antici, E. Boella, S. N. Chen, D. S. Andrews, M. Barberio, J. Böker, F. Cardelli, J. L. Feugeas, M. Glesser, P. Nicolaï, L. Romagnani, M. Scisciò, M. Starodubtsev, O. Willi, J. C. Kieffer, V. Tikhonchuk, H. Pépin, L. O. Silva, E. D. Humières, and J. Fuchs, Scientific Reports 7, 16463 (2017), arXiv:1708.02539 .
  * Chen  _et al._ (2017) S. N. Chen, M. Vranic, T. Gangolf, E. Boella, P. Antici, M. Bailly-Grandvaux, P. Loiseau, H. Pépin, G. Revet, J. J. Santos, A. M. Schroer, M. Starodubtsev, O. Willi, L. O. Silva, E. D’humières, and J. Fuchs, Scientific Reports 7, 1 (2017).
  * Pak  _et al._ (2018) A. Pak, S. Kerr, N. Lemos, A. Link, P. Patel, F. Albert, L. Divol, B. B. Pollock, D. Haberberger, D. Froula, M. Gauthier, S. H. Glenzer, A. Longman, L. Manzoor, R. Fedosejevs, S. Tochitsky, C. Joshi, and F. Fiuza, Physical Review Accelerators and Beams 21, 103401 (2018), arXiv:1810.08190 .
  * Polz  _et al._ (2019) J. Polz, A. P. Robinson, A. Kalinin, G. A. Becker, R. A. Fraga, M. Hellwing, M. Hornung, S. Keppler, A. Kessler, D. Klöpfel, H. Liebetrau, F. Schorcht, J. Hein, M. Zepf, R. E. Grisenti, and M. C. Kaluza, Scientific Reports 9, 16534 (2019).
  * Ota  _et al._ (2019) M. Ota, A. Morace, R. Kumar, S. Kambayashi, S. Egashira, M. Kanasaki, Y. Fukuda, and Y. Sakawa, High Energy Density Physics 33, 100697 (2019).
  * Kumar  _et al._ (2019) R. Kumar, Y. Sakawa, L. N. K. Döhl, N. Woolsey, and A. Morace, Physical Review Accelerators and Beams 22, 043401 (2019).
  * Forslund and Shonk (1970) D. Forslund and C. Shonk, Physical Review Letters 25, 281 (1970).
  * Forslund  _et al._ (1971) D. W. Forslund, R. L. Morse, and C. W. Nielson, Physical Review Letters 27, 1424 (1971).
  * Balogh, A and Treumann (2013) R. A. Balogh, A and Treumann,  _Physics of Collisionless Shocks: Space Plasma Shock Waves_ , Vol. 12 (2013).
  * Drury (1983) L. O. Drury, Reports on Progress in Physics 46 973(1983).
  * JONES (1991) F. C. JONES, Space Science Reviews 58, 259 (1991).
  * Ohira and Takahara (2008) Y. Ohira and F. Takahara, The Astrophysical Journal 688, 320 (2008), arXiv:0808.3195 .
  * Treumann (2009) R. A. Treumann, The Astronomy and Astrophysics Review 17, 409 (2009).
  * Madanian  _et al._ (2020) H. Madanian, S. J. Schwartz, J. S. Halekas, and L. B. Wilson, Geophysical Research Letters 47, e2020GL088309 (2020).
  * Bale and Mozer (2007) S. D. Bale and F. S. Mozer, Physical Review Letters 98, 205001 (2007).
  * Bale  _et al._ (2005) S. D. Bale, P. J. Kellogg, F. S. Mozer, T. S. Horbury, and H. Reme, Physical Review Letters 94, 215002 (2005).
  * Zhang  _et al._ (2009) X. Zhang, B. Shen, L. Ji, F. Wang, Z. Jin, X. Li, M. Wen, and J. R. Cary, Physical Review Special Topics - Accelerators and Beams 12, 021301 (2009).
  * Arber  _et al._ (2015) T. D. Arber, K. Bennett, C. S. Brady, A. Lawrence-Douglas, M. G. Ramsay, N. J. Sircombe, P. Gillies, R. G. Evans, H. Schmitz, A. R. Bell, and C. P. Ridgers, Plasma Physics and Controlled Fusion 57 , 113001 (2015), 10.1088/0741-3335/57/11/113001.
  * Grismayer and Mora (2006) T. Grismayer and P. Mora, Physics of Plasmas 13, 032103 (2006).
  * Tidman and Krall (1971) D. A. Tidman and N. A. Krall,  _Shock Waves in Collisionless Plasmas_ (Wiley-Interscience, New York, 1971).
  * Liseykina  _et al._ (2015) T. V. Liseykina, G. I. Dudnikova, V. A. Vshivkov, and M. A. Malkov, Journal of Plasma Physics 81, 495810507 (2015).
  * Robinson  _et al._ (2009) A. P. Robinson, P. Gibbon, M. Zepf, S. Kar, R. G. Evans, and C. Bellei, Plasma Physics and Controlled Fusion 51 1024004 (2009), 10.1088/0741-3335/51/2/024004.


[a]Weidong Jin

# VERITAS follow-up observation of the blazar TXS 0506+056

on behalf of the VERITAS Collaboration  RileyAnne Sharpe

###### Abstract

The gamma-ray blazar TXS 0506+056 was found with an enhanced gamma-ray
emission state in spatial and temporal coincidence with the IceCube high
energy neutrino event IC170922A. This is the most significant association by
far between a high-energy neutrino event and a blazar in a flaring state.
Studying the time evolution and spectral behavior of the blazar emission may
help in identifying the sources of the diffuse neutrino flux observed by
IceCube and the origin of energetic cosmic rays. TXS 0506+056 was detected by
the VERITAS gamma-ray observatory with a significance of 5.8 standard
deviations above 110 GeV in a 35 hour data set collected between September 23,
2017 and February 6, 2018. Here we will present results from recent VERITAS
observations and an associated multiwavelength campaign, collected between
October 10, 2018 to March 1, 2021. A relatively quiet very high energy gamma-
ray emission state was observed during this time period, and flux upper limits
are used to constrain the potential variability of this blazar.

##  1 Introduction

On September 22, 2017, the IceCube Collaboration reported the detection of a
high-energy muon neutrino (290 TeV) event, labeled IceCube-170922A, of
potential astrophysical origin. This event was sent out as an alert [1] and
prompted a multiwavelength campaign from the radio to very-high-energy (VHE, E
$>$ 100 GeV) gamma-ray band. The location of the neutrino was consistent with
that of the blazar TXS 0506+056 (z = 0.3365 $\pm$ 0.0010) which at the time
was observed by _Fermi_ -LAT to be flaring in gamma rays. The correlation of
the neutrino detection with the gamma-ray flare of TXS 0506+056 is
statistically significant at the 3$\sigma$ level [2]. Follow-up observations
led to the first detection of this blazar in the VHE emission region by the
MAGIC telescopes [3]. An initial observation taken by VERITAS did not yield a
detection, but a follow-up observation taken between September 23, 2017 and
February 6, 2018 (MJD 58019-58155) yielded a detection of the source above 110
GeV with a statistical significance of 5.8$\sigma$ based on a 35-hr data set
[4].

Although no significant excess of cumulative neutrino emission was found by
searching gamma-ray blazars detected by 2nd _Fermi_ -LAT AGN catalog (2LAC),
and the maximum contribution from blazars in 2LAC to the all-sky astrophysical
neutrino flux is constrained to be no more than $\sim$27% [5], this did not
rule out that individual blazars could be potential neutrino counterparts.
Since 2017, TXS 0506+056 remains the most significant correlation between a
blazar and a high-energy neutrino. More recently, other blazars such as 4FGL
J0658.6+0636 (more details in this conference in [6]) has also been identified
in the 90% uncertainty region of well-reconstructed IceCube neutrino events.
Thus TXS 0506+056 plays a key role in establishing connections between high-
energy neutrinos and astrophysical sources. Studying the time evolution and
spectral behavior of this blazar may help in identifying the sources of the
diffuse neutrino flux observed by IceCube and the origin of energetic cosmic
rays.

##  2 VERITAS observation of blazar TXS 0506+056 from 2018 to 2021

VERITAS [7] consists of an array of four 12-m imaging atmospheric Cherenkov
telescopes located at the Fred Lawrence Whipple Observatory (FLWO) in southern
Arizona, USA (31∘ 40’ N, 110∘ 57’ W, 1.3 km a.s.l.). Each telescope equipped
with a camera containing 499 photomultiplier tubes covering a field of view of
3.5∘. VERITAS has maximum sensitivity in the range of 100 GeV to 30 TeV. The
angular resolution of VERITAS is $\sim 0.1^{\circ}$ at 1 TeV (for 68%
containment) and the energy resolution is 15-25% at the same energy. In its
current configuration, VERITAS can detect a source ($5\sigma$ significance)
with a flux of 1% of the steady flux from the Crab Nebula in less than 25
hours [8]. VERITAS data were analyzed using standard analysis tools [9, 10,
11].

VERITAS has implemented a wide-ranging neutrino follow-up program to search
for VHE gamma rays associated with IceCube neutrino events [12], including
follow-up of real-time IceCube neutrino alerts (more details at this
conference in [13]) and long-term monitoring of the candidate neutrino blazar
TXS 0506+056. VERITAS collected 61 hours of quality-selected observations of
TXS 0506+056 from from Oct 10, 2018 (MJD 58401) to March 7, 2021 (MJD 59280)
with an average zenith angle of 28.8∘. Observations were performed using the
standard “wobble” observation mode [14] with a 0.5∘ offset in each of four
cardinal directions. The analysis yielded a detection of the source above 190
GeV with a statistical significance of 3.4$\sigma$. The integral flux above an
energy threshold of 190 GeV is $(1.34\pm 0.40)\times 10^{-12}$ cm-2 s-1, which
corresponds to $(0.52\pm 0.16)\%$ of the Crab Nebula flux (C.U.) under the
same energy threshold [15, 16]. Fig. 1 shows a significance sky map with TXS
0506+056 in the center. A light curve of the integral flux above an energy
threshold of 190 GeV is shown in Fig. 2. Each point represents the average
flux level of the source with one month per bin. The vertical gray lines
represent individual observations. A 95$\%$ confidence level upper limit is
calculated when significance is less than 2$\sigma$. For the upper limit
calculation, this analysis used the best fitting spectral index ($\Gamma$=4.8)
obtained from VERITAS’s follow-up observation of TXS 0506+056 from year 2017
to 2018 [4]. A relatively quiet VHE gamma-ray emission state was observed.

Figure 1: VERITAS significance sky map with TXS 0506+056 in the center.

##  3 Multiwavelength observations

The multiwavelength light curves are composed of high energy gamma-ray data
from the _Fermi_ -LAT light curve repository
(LCR)111https://fermi.gsfc.nasa.gov/ssc/data/access/lat/LightCurveRepository/,
X-ray data from the Neil Gehrels _Swift_ Observatory X-ray telescope
(XRT)222https://www.swift.ac.uk/user_objects/ and Nuclear Spectroscopic
Telescope Array (NuSTAR), as well as optical data from the All-Sky Automated
Survey for Supernovae (ASAS-SN) Sky Patrol333https://asas-sn.osu.edu/ (Fig.
2). The LCR is a public library of light curves for sources in the 10 year
Fermi LAT point source (4FGL) catalog [17]. The data are presented on time
scales of weeks and photon energy between 100 MeV and 100 GeV. The observation
made by ASAS-SN Sky Patrol are available in two optical band: the _V_
(original two stations) or _g_ (three new stations) band filters and three
dithered 90s exposures [18, 19].

NuSTAR and _Swift_ observations were collected during a monitoring campaign of
TXS 0506+056 to support the VERITAS observations in search for evidence of
hadronic emission. The campaign, organized during NuSTAR GO Cycle 5 between
MJD 58550 and 58950, consisted of 4 $\times$ 20 ks NuSTAR exposures separated
by $\sim 60$ days to sample similar time scales as the candidate neutrino
“flare” of 2014-2015 and was supplemented by 10 $\times$ 3 ks _Swift_
exposures interleaved with the NuSTAR epochs. The epochs were selected to
coincide with the middle of the VERITAS dark observation period so that quasi-
simultaneous observations were possible. The NuSTAR and _Swift_ -XRT
observations were reduced using the nupipeline and xrtpipeline routines,
respectively, available in HEASoft 6.28. Background-subtracted count rate
light curves were derived from the reduced data in two energy ranges: 0.3 - 10
keV for _Swift_ -XRT and 5-80 keV for NuSTAR. The light curves shown in Fig. 2
show no indication of strong flaring in X-rays with the exception of an
elevated state observed by _Swift_ on MJD 58902 with no apparent counterpart
in _Fermi_. The flux level that corresponds to this rate ($\sim 4\times
10^{-12}$ erg cm-2 s-1) did not meet the pre-defined criteria to trigger an
additional NuSTAR exposure during the monitoring program and no coincident
VERITAS observations were possible at the time as it happened outside the
VERITAS dark observing period. The detailed analysis of the X-ray data is
ongoing.

Figure 2: Multiwavelength light curves of blazar TXS 0506+056 from Oct 10,
2018 to March 7, 2021. The data set is composed of VHE gamma-ray data from
VERITAS, high-energy gamma-ray data from _Fermi_ -LAT LCR, X-ray data from
_Swift_ and NuSTAR, and optical data (not corrected for extinction) from ASAS-
SN Sky Patrol.

##  4 Summary and outlook

VERITAS’s previous follow-up observations of TXS 0506+056 during 2017/2018
season gave a flux of $(1.79\pm 0.56)\times 10^{-12}$ cm-2 s-1, corresponding
to $(0.70\pm 0.22)\%$ of C.U. (the analysis was redone with energy threshold
corrected to 190 GeV and updated instrument response functions) [4]. TXS
0506+056 is found to be in a quiet state from 2018 to 2021 and shows a
consistent flux level $((0.52\pm 0.16)\%$ of C.U.) compared to previous
VERITAS follow-up observations. MAGIC detected enhanced VHE gamma-ray
emissions on Dec $1^{st}$ (MJD 58453) and $3^{rd}$ 2018 (MJD 58455) and its
flux was reported as comparable with the flare detected by MAGIC in October
2017 [20]. This enhanced emission was not detected by VERITAS because VERITAS
didn’t have observation on the first date and the data quality is bad on the
second date due to poor weather condition. In addition to that, no extended
VHE gamma-ray emissions were detected by VERITAS, consistent with a low state
of the source during this period. Clear variability is observed in optical,
X-ray and high-energy gamma-ray bands. Quantitative analysis will be presented
in an upcoming publication.

## Acknowledgements

This research is supported by grants from the U.S. Department of Energy Office
of Science, the U.S. National Science Foundation and the Smithsonian
Institution, by NSERC in Canada, and by the Helmholtz Association in Germany.
This research used resources provided by the Open Science Grid, which is
supported by the National Science Foundation and the U.S. Department of
Energy’s Office of Science, and resources of the National Energy Research
Scientific Computing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility operated under Contract No. DE-AC02-05CH11231. We
acknowledge the excellent work of the technical support staff at the Fred
Lawrence Whipple Observatory and at the collaborating institutions in the
construction and operation of the instrument.

This research was partially supported by NASA grant NuSTAR GO-5277. This
research has made use of data obtained with NuSTAR, a project led by Caltech,
funded by NASA and managed by NASA/JPL, and has utilized the NUSTARDAS
software package, jointly developed by the ASDC (Italy) and Caltech (USA).

## References

  * [1] C. Kopper and E. Blaufuss, _IceCube-170922A - IceCube observation of a high-energy neutrino candidate event._ , _GRB Coordinates Network_ 21916 (2017) 1.
  * [2] The IceCube Collaboration et al., _Multimessenger observations of a flaring blazar coincident with high-energy neutrino icecube-170922a_ , _Science_ 361 (2018) .
  * [3] S. Ansoldi, L.A. Antonelli, C. Arcaro, D. Baack, A. Babić, B. Banerjee et al., _The blazar txs 0506+ 056 associated with a high-energy neutrino: Insights into extragalactic jets and cosmic-ray acceleration_ , _The Astrophysical Journal Letters_ 863 (2018) L10.
  * [4] A.U. Abeysekara, A. Archer, W. Benbow, R. Bird, A. Brill, R. Brose et al., _Veritas observations of the bl lac object txs 0506+056_ , _The Astrophysical Journal_ 861 (2018) L20.
  * [5] M.G. Aartsen, K. Abraham, M. Ackermann, J. Adams, J.A. Aguilar, M. Ahlers et al., _The contribution of fermi-2lac blazars to diffuse tev–pev neutrino flux_ , _The Astrophysical Journal_ 835 (2017) 45.
  * [6] R. de Menezes et al., _Multi-Messenger observations of the $\gamma$-ray blazar 4FGL J0658.6+0636 consistent with an IceCube high-energy neutrino_, in _37th International Cosmic Ray Conference (ICRC2021)_ , International Cosmic Ray Conference, 2021.
  * [7] J. Holder et al., _The first VERITAS telescope_ , _Astroparticle Physics_ 25 (2006) 391 [astro-ph/0604119].
  * [8] N. Park et al., _Performance of the VERITAS experiment_ , _Proc. of 34th ICRC, Vol 34, 771_ 34 (2015) 771 [1508.07070].
  * [9] G. Maier and J. Holder, _Eventdisplay: An Analysis and Reconstruction Package for Ground-based Gamma-ray Astronomy_ , in _35th International Cosmic Ray Conference (ICRC2017)_ , vol. 301 of _International Cosmic Ray Conference_ , p. 747, Jan., 2017 [1708.04048].
  * [10] P. Cogan, _VEGAS, the VERITAS Gamma-ray Analysis Suite_ , in _International Cosmic Ray Conference_ , vol. 3 of _International Cosmic Ray Conference_ , pp. 1385–1388, Jan., 2008 [0709.4233].
  * [11] M.K. Daniel, _The VERITAS standard data analysis_ , in _International Cosmic Ray Conference_ , vol. 3 of _International Cosmic Ray Conference_ , pp. 1325–1328, Jan., 2008 [0709.4006].
  * [12] M. Santander, _Recent results from the VERITAS multi-messenger program_ , in _36th International Cosmic Ray Conference (ICRC2019)_ , vol. 36 of _International Cosmic Ray Conference_ , p. 782, July, 2019, https://arxiv.org/abs/1909.05228 [1909.05228].
  * [13] K. Satalecka et al., _Searching for VHE gamma-ray emission associated with IceCube neutrino alerts using FACT, H.E.S.S., MAGIC, and VERITAS_ , in _37th International Cosmic Ray Conference (ICRC2021)_ , International Cosmic Ray Conference, 2021.
  * [14] V.P. Fomin, A.A. Stepanian, R.C. Lamb, D.A. Lewis, M. Punch and T.C. Weekes, _New methods of atmospheric Cherenkov imaging for gamma-ray astronomy. I. The false source method_ , _Astroparticle Physics_ 2 (1994) 137.
  * [15] A.M. Hillas, C.W. Akerlof, S.D. Biller, J.H. Buckley, D.A. Carter-Lewis, M. Catanese et al., _The Spectrum of Teravolt Gamma Rays from the Crab Nebula_ , _The Astrophysical Journal_ 503 (1998) 744.
  * [16] K. Meagher and VERITAS Collaboration, _Six years of VERITAS observations of the Crab Nebula_ , in _34th International Cosmic Ray Conference (ICRC2015)_ , vol. 34 of _International Cosmic Ray Conference_ , p. 792, July, 2015 [1508.06442].
  * [17] S. Abdollahi et al., _Fermi Large Area Telescope Fourth Source Catalog_ , _The Astrophysical Journal Supplement Series_ 247 (2020) 33 [1902.10045].
  * [18] B.J. Shappee, J.L. Prieto, D. Grupe, C.S. Kochanek, K.Z. Stanek, G. De Rosa et al., _The man behind the curtain: X-rays drive the uv through nir variability in the 2013 active galactic nucleus outburst in ngc 2617_ , _The Astrophysical Journal_ 788 (2014) 48.
  * [19] C.S. Kochanek, B.J. Shappee, K.Z. Stanek, T.W.S. Holoien, T.A. Thompson, J.L. Prieto et al., _The All-Sky Automated Survey for Supernovae (ASAS-SN) Light Curve Server v1.0_ , _Publications of the Astronomical Society of the Pacific_ 129 (2017) 104502 [1706.07060].
  * [20] K. Satalecka, E. Bernardini, W. Bhattacharyya, M. Cerruti, V. Fallah Ramazani, L. Foffano et al., _MAGIC and MWL monitoring of the blazar TXS 0506+056 in the 2018/2019 season_ , in _Proceedings of 36th International Cosmic Ray Conference — PoS(ICRC2019)_ , vol. 358, p. 783, 2019, DOI.

## Full Authors List: VERITAS Collaboration

C. B. Adams1, A. Archer2, W. Benbow3, A. Brill1, J. H. Buckley4, M. Capasso5,
J. L. Christiansen6, A. J. Chromey7, M. Errando4, A. Falcone8, K. A. Farrell9,
Q. Feng5, G. M. Foote10, L. Fortson11, A. Furniss12, A. Gent13, G. H.
Gillanders14, C. Giuri15, O. Gueta15, D. Hanna16, O. Hervet17, J. Holder10, B.
Hona18, T. B. Humensky1, W. Jin19, P. Kaaret20, M. Kertzman2, T. K. Kleiner15,
S. Kumar16, M. J. Lang14, M. Lundy16, G. Maier15, C. E McGrath9, P.
Moriarty14, R. Mukherjee5, D. Nieto21, M. Nievas-Rosillo15, S. O’Brien16, R.
A. Ong22, A. N. Otte13, S. R. Patel15, S. Patel20, K. Pfrang15, M. Pohl23,15,
R. R. Prado15, E. Pueschel15, J. Quinn9, K. Ragan16, P. T. Reynolds24, D.
Ribeiro1, E. Roache3, J. L. Ryan22, I. Sadeh15, M. Santander19, G. H.
Sembroski25, R. Shang22, D. Tak15, V. V. Vassiliev22, A. Weinstein7, D. A.
Williams17, and T. J. Williamson10

1Physics Department, Columbia University, New York, NY 10027, USA 2Department
of Physics and Astronomy, DePauw University, Greencastle, IN 46135-0037, USA
3Center for Astrophysics $|$ Harvard & Smithsonian, Cambridge, MA 02138, USA
4Department of Physics, Washington University, St. Louis, MO 63130, USA
5Department of Physics and Astronomy, Barnard College, Columbia University, NY
10027, USA 6Physics Department, California Polytechnic State University, San
Luis Obispo, CA 94307, USA 7Department of Physics and Astronomy, Iowa State
University, Ames, IA 50011, USA 8Department of Astronomy and Astrophysics, 525
Davey Lab, Pennsylvania State University, University Park, PA 16802, USA
9School of Physics, University College Dublin, Belfield, Dublin 4, Ireland
10Department of Physics and Astronomy and the Bartol Research Institute,
University of Delaware, Newark, DE 19716, USA 11School of Physics and
Astronomy, University of Minnesota, Minneapolis, MN 55455, USA 12Department of
Physics, California State University - East Bay, Hayward, CA 94542, USA
13School of Physics and Center for Relativistic Astrophysics, Georgia
Institute of Technology, 837 State Street NW, Atlanta, GA 30332-0430 14School
of Physics, National University of Ireland Galway, University Road, Galway,
Ireland 15DESY, Platanenallee 6, 15738 Zeuthen, Germany 16Physics Department,
McGill University, Montreal, QC H3A 2T8, Canada 17Santa Cruz Institute for
Particle Physics and Department of Physics, University of California, Santa
Cruz, CA 95064, USA 18Department of Physics and Astronomy, University of Utah,
Salt Lake City, UT 84112, USA 19Department of Physics and Astronomy,
University of Alabama, Tuscaloosa, AL 35487, USA 20Department of Physics and
Astronomy, University of Iowa, Van Allen Hall, Iowa City, IA 52242, USA
21Institute of Particle and Cosmos Physics, Universidad Complutense de Madrid,
28040 Madrid, Spain 22Department of Physics and Astronomy, University of
California, Los Angeles, CA 90095, USA 23Institute of Physics and Astronomy,
University of Potsdam, 14476 Potsdam-Golm, Germany 24Department of Physical
Sciences, Munster Technological University, Bishopstown, Cork, T12 P928,
Ireland 25Department of Physics and Astronomy, Purdue University, West
Lafayette, IN 47907, USA


# Pulsar timing array results sheds light on Hubble tension during the end of
inflation

M. Bousder1, E. Salmani1,2, A.Riadsolh${}^{3},$<EMAIL_ADDRESS>H. Ez-Zahraouy1,2, A. El Fatimy4,5 and M. El Belkacemi3
1Laboratory of Condensend Matter and Interdisplinary Sciences, Department of
physics,
Faculty of Sciences, Mohammed V University in Rabat, Morocco
2CNRST Labeled Research Unit (URL-CNRST), Morocco
3Laboratory of Conception and Systems (Electronics, Signals and Informatics)
Faculty of Sciences, Mohammed V University in Rabat, Morocco
4Central European Institute od Technology,
CEITEC BUT, Purkyňova 656/123, 61200 Brno, Czech Republic
5Departement of Physics, Université Mohammed VI Polytechnique, Ben Guerir
43150, Morocco

###### Abstract

Recently, pulsar timing array (PTA) collaborations, including NANOGrav, have
reported evidence of a stochastic gravitational wave background within the nHz
frequency range. It can be interpreted by gravitational waves from preheating
era. In this context, we demonstrate that the emission of this stochastic
gravitational wave background can be attributed to fluctuations occurring at
the end of inflation, thus giving rise to the Hubble tension issue. At the
onset of inflation, the value of the frequency of the gravitational wave
signal stood at $f=0.08nHz$, but it rapidly transitioned to $f=1nHz$ precisely
at the end of inflation. However, just before the end of inflation, a phase
characterized by curvature perturbation is known to occur, causing a swift
increase in the frequency.

Keywords: Pulsar timing array, gravitational wave, inflation, Hubble tension.

##  1 Introduction

The stochastic gravitational wave (GW) signals at nHZ frequencies in the
15-year data set [1, 2] and 12.5-year data [3, 4], has recently garnered
significant attention. This data has been published by the North American
Nanohertz Observatory for Gravitational Waves (NANOGrav) collaboration. The GW
signal may arise following the power-law distribution of abundance:
$\Omega_{GW}\propto f^{\zeta}$. NANOGrav’s observations indicate that the
exponent $\zeta$ lies within the range of $\left(-1.5,0.5\right)$ at a
frequency of $f=5.5nHz$. In response to the NANOGrav results, a series of
recent gravitational wave models, [5, 6, 7], have been proposed. Often the
authors linked this signal with cosmic inflation [8, 9] and from axion
inflation [10]. This signal has the potential to be construed as an outcome of
the stochastic gravitational wave background emanating from primordial black
holes (PBHs) formed during inflation, [11, 12, 13, 14, 15, 16]. Conversely,
alternative studies propose that cosmic strings might serve as the origin of
these waves, as indicated in [17].
The authors in [18] explore gravitational wave production during the end of
inflation (preheating). This process generates significant inhomogeneities,
leading to a stochastic background of gravitational waves within the comoving
Hubble horizon at the end of inflation. Importantly, their results emphasize
that the current amplitude of these gravitational waves is independent of the
inflationary energy scale. The analysis focuses on a specific model with an
inflationary energy scale of approximately $10^{9}GeV$. This discovery
highlights the potential for a new observational perspective on inflationary
physics and encourages further investigations into stochastic gravitational
wave backgrounds across the Hz to GHz frequency range, characterized by an
amplitude of $\Omega_{GW}h^{2}\sim 10^{-11}$. Moreover, their computational
approach has broader applications for understanding gravitational waves
originating from various inhomogeneous processes in the early universe. The
paper [19] considers a scenario in which the end of inflation is treated as a
sudden event, and it provides mathematical formulations for essential
parameters such as the spectrum, spectral tilt, and non-Gaussianity. These
formulas are then applied to analyze a minimal extension of the initial hybrid
inflation model, shedding light on the intricate dynamics of the early
universe.
The Measurements of the current rate of expansion of the Universe conducted on
a local scale, utilizing observations of Type Ia Supernovae (SNeIa) through
the Hubble Space Telescope (HST), exhibit values that tend to cluster around
$H_{0}\simeq\left(73.24\pm 1.74\right)km.s^{-1}.Mpc^{-1}$ [20]. In contrast,
calculations of the Hubble constant derived from measurements of the Cosmic
Microwave Background (CMB) indicate a value of approximately
$H_{0}\simeq\left(67.27\pm 0.66\right)km.s^{-1}.Mpc^{-1}$ [21]. This
introduces a deviation of 4$\sigma$ or more from the local measurements. One
of the more effective resolutions to the Hubble tension [22, 23, 24, 25, 26]
involves incorporating an early dark energy component [27] or introducing
primordial magnetic fields [28]. Additionally, numerous studies have
approached this issue by imposing constraints.
The purpose of the present paper is to introduce an interpretation of the
Hubble tension from NANOGrav signal at the end of inflation.
This paper is organized as follows. In Sec. 2, we review the stochastic
gravitational wave background from the tensor perturbation and the
gravitational wave spectrum. In Sec. 3, we analyze the Hubble tension in
inflationary epoch. In Sec. 4, the e-folding number, slow roll parameter and
spectral index are studied. In Sec. 5, we study the energy density spectrum in
the end of inflation. We conclude our findings in Sec. 6.

##  2 Amplitude of the GW signal

The action of the gravitational wave sources represented by the tensor
perturbation $h_{ij}$ $(i,j=1,2,3)$ is

$S=\int d\tau d^{3}x\sqrt{-\delta}\left[\frac{-\delta^{\mu\nu}}{64\pi
G}\partial_{\mu}h_{ij}\partial_{\nu}h^{ij}+\frac{1}{2}\Pi_{ij}h^{ij}\right],$
(2.1)

where $G$ is the Newtonian gravitational constant, $\tau$ is the conformal
time and $\Pi_{ij}$ corresponds to an anisotropic stress, encompassing the
impact of a macroscopic magnetic field. Here $\delta_{\mu\nu}\ $in the
Minkowski metric of flat spacetime [29]
$g_{\mu\nu}=\delta_{\mu\nu}+h_{\mu\nu}$ $(\mu,\nu=0,1,2,3).$
The Friedmann-Robertson-Walker (FRW) spacetime [30]:

$ds^{2}=\left(ic^{2}dt\right)^{2}-a^{2}(t)\left(\delta_{ij}+h_{ij}\right)dx^{i}dx^{j},$
(2.2)

where $a\left(t\right)$ is the scale factor with $ad\tau=dt$. The Hubble
parameter is defined as $H=\dot{a}/a$ where $\dot{a}$ signifies the derivative
of $a$ with respect to time. we write the conditions on the tensor $h_{ij}$ as
$\partial^{i}h_{ij}=h_{\text{ \ }i}^{i}=0$.
We can perform a Fourier transformation of the gravitational wave field by
[31, 32, 33]:

$h_{ij}\left(\tau,\vec{x}\right)=\int\frac{d^{3}\vec{k}}{\left(2\pi\right)^{3/2}}\sum_{s=+,\times}\left[h_{k}\left(\tau\right)\epsilon_{ij}C\left(\vec{k}\right)e^{i\vec{k}\vec{x}}+cc\right],$
(2.3)

where $k=\left|\vec{k}\right|=2\pi f$ is the GW wavenumber and $cc$ stands for
the complex conjugate term. Here, the $\epsilon_{ij}$ satisfy the conditions
$\delta^{ij}\epsilon_{ij}=0$. The energy density of the gravitational waves is
expressed as follows [30]:

$\rho_{GW}=\frac{1}{64\pi
G}\left\langle\left(\partial_{t}h_{ij}\right)^{2}+\frac{1}{a^{2}}\left(\nabla
h_{ij}\right)^{2}\right\rangle,$ (2.4)

where the bracket $\left\langle\bullet\right\rangle$ describes the spatial
average. To describe the spectral amplitude of GWs, we employ the
dimensionless parameter $\Omega_{GW}=\rho_{GW}/\rho_{crit}$ where
$\rho_{crit}=\frac{3H_{0}^{2}}{8\pi G}$ . This parameter signifies the energy
density of GWs within each logarithmic frequency interval:

$\Omega_{GW}\left(f\right)=\frac{1}{\rho_{crit}}\frac{d\rho_{GW}}{d\ln f},$
(2.5)

where $H_{0}$ is the current Hubble constant and $f$ is the frequency of the
GW signal. An approximation for the gravitational wave spectrum can be made by
employing the power-law model:

$\Omega_{GW}\left(f\right)=\Omega_{GW\ast}\left(\frac{f}{f_{\ast}}\right)^{n_{s}\left(f\right)}\text{,}$
(2.6)

where $\Omega_{GW\ast}=\Omega_{GW}\left(f=f_{\ast}\right)$. The spectral index
is

$n_{s}\left(f\right)=\left\\{\begin{array}[]{c}n_{GW1}\text{ \ \ \ for \ \
}f<f_{\ast}\\\ n_{GW2}\text{ \ \ \ \ for \ }f>f_{\ast}\end{array}\right..$
(2.7)

We notice that for $f=f_{\ast}$, we get $n_{GW1}=n_{GW2}$. If $n_{s}=-1$, we
obtain $f\Omega_{GW}=f_{\ast}\Omega_{GW\ast}$.

##  3 Hubble tension from GW signal

Hubble tension from GW signal [34, 13]. The current amplitude of the
gravitational wave signal observed today is given by [35, 36]

$\Omega_{GW}\left(f\right)=\frac{\Omega_{R,0}}{24}\mathcal{P}_{h}\left(f\right),$
(3.1)

where $\Omega_{R,0}\approx 8.6\times 10^{-5}$ is the radiation energy density
today and $\mathcal{P}_{h}\left(f\right)$ is the power spectrum of tensor
fluctuations, which varies with frequency, at the moment of exiting the cosmic
horizon. We can establish a connection between the e-folding number $N$ and
the GW frequency of the signal $f$, such that [35]:

$N=N_{CMB}+\ln\frac{k_{CMB}}{0.002Mpc^{-1}}-44.9-\ln\frac{f}{10^{2}Hz}.$ (3.2)

We notice that $N=0$ at the end of inflation. While, in the CMB scale we have
$k_{CMB}=0.002Mpc^{-1}$and $N_{CMB}\sim 50-60$, in this case we get

$N=N_{CMB}-44.9-\ln\frac{f}{10^{2}Hz}.$ (3.3)

Or equivalently

$N=\ln\frac{10^{2}e^{N_{CMB}-44.9}Hz}{f}.$ (3.4)

Thus, the e-folding number can be expressed as

$N=\ln\frac{f_{\ast}}{f},$ (3.5)

where $f_{\ast}=10^{2}e^{\left(N_{CMB}-44.9\right)}Hz$ is the typical
frequency. Based on this relation, it’s noteworthy that the state represented
by $f_{\ast}=f$ describes the end of inflation. Conversely, in order to align
the value of $N$ with observational data, we need to take into account the
following condition:

$f<f_{\ast}.$ (3.6)

This condition is comparable with $n_{s}=n_{GW1}$ (2.7). On the other hand,
the equations (2.6) and (3.5) take this form,

$\Omega_{GW}\left(f\right)=\Omega_{GW\ast}e^{n_{s}\left(f\right)N\text{ }}.$
(3.7)

At the end of inflation, we find
$\Omega_{GW}\left(f\right)\approx\Omega_{GW\ast}$. The Planck CMB data implies
that $n_{s}\approx 0.965\pm 0.004$ [21]. It is easy to see from $N_{CMB}\sim
50$ that

$\Omega_{GW}\left(f\right)=\Omega_{GW\ast}e^{48.25\pm 0.2\text{ }}.$ (3.8)

For $\Omega_{GW}\approx 9.3_{-4.0}^{+5.8}\times 10^{-9}$ [1, 2] we have the
critical density of GWs as follows $\Omega_{GW\ast}\sim 10^{-29}-10^{-30}$.
The number of e-folds can be given by [39, 40]:

$N=\frac{1}{\epsilon_{H}}\int_{H}^{H_{end}}d\ln
H^{\prime}=\frac{1}{\epsilon_{H}}\ln\frac{H_{end}}{H}.$ (3.9)

Here, $\epsilon_{H}=-\frac{\dot{H}}{H^{2}}\ll 1$ is the slow roll parameter.
This relationship corresponds with observations data, because if we take
$H=H_{end}$ we find the value of $N$ at the end of inflation ($N=0$). We can
additionally calculate the $H$ concerning $f$ using the following approach:

$H\left(f\right)=H_{end}\left(\frac{f}{f_{\ast}}\right)^{\epsilon_{H}}.$
(3.10)

In the scenario where $f_{\ast}=f$, we arrive at $H\left(f\right)=H_{end}$.
Conversely, in the case of (3.6), we achieve

$H\left(f\right)<H_{end}.$ (3.11)

We recall that the assessment of the present rate of the Universe’s expansion
yields $H_{0}\simeq\left(73.24\pm 1.74\right)km.s^{-1}.Mpc^{-1}$ [20].
Conversely, measurements of the Hubble constant based on the CMB result is
$H_{0}\simeq\left(67.27\pm 0.66\right)km.s^{-1}.Mpc^{-1}$ [21]. The Hubble
tension’s behavior reveals a diminishing value of $H$ as time elapses. This
shows that according to (3.11) the Hubble parameter $H\left(f\right)$
characterizes the inflationary epoch. This is underscored by the fact that
$H_{end}$ signifies the end of inflation. Consequently, the frequencies $f$
detected by the NANOGrav collaboration are indicative of the inflationary era.

##  4 Inflationary scenario from e-folding number

Thereafter we want to seek the relationship between the current amplitude of
the gravitational wave signal and frequency. We recall that $\dot{N}=-H$ [41],
i.e. $\dot{N}=-\frac{\dot{f}}{f}$ (3.5). Consequently, this leads to

$H=\frac{\dot{f}}{f}.$ (4.1)

This result implies that the Hubble parameter represents the logarithmic
variation of GW frequency. If we compare (4.1) with the definition of the
Hubble parameter $H=\frac{\dot{a}}{a}$, we find $f=af_{0}$ with $f_{0}$ is a
constant with frequency unit. This proportionality between the scale factor
and frequency can establish a link between the Hubble tension problem and the
NANOGrav 15-yr result. The slow roll parameter
$\epsilon_{H}=-\frac{\dot{H}}{H^{2}}$ can be written as

$\epsilon_{H}=1-\frac{\ddot{f}f}{\dot{f}^{2}}.$ (4.2)

Introducing these parameters

$\sigma=\frac{\ddot{f}}{\dot{f}^{2}}\text{, \
}\lambda=\frac{\dddot{f}}{\ddot{f}}\frac{1}{\dot{f}}.$ (4.3)

From this model, the slow-roll parameter
$\eta=\epsilon_{H}-\frac{\dot{\epsilon}_{H}}{2H\epsilon_{H}}$ [41] is given by

$\eta=1-\sigma f-\frac{f}{2\left(1-\sigma
f\right)}\left(2\sigma^{2}f-\lambda\sigma f-\sigma\right).$ (4.4)

In the inflationary scenario, the $\epsilon_{H}$ grows as well
$d\epsilon_{H}/dN=2\epsilon_{H}\left(\eta-\epsilon_{H}\right)$. The spectral
index of perturbations is equal to [41]:

$n_{s}=1-6\epsilon_{H}+2\eta.$ (4.5)

Thus, the spectral index (see Appendix A) can be written as

$n_{s}=-3+4\sigma f-\frac{f}{1-\sigma f}\left(2\sigma^{2}f-\lambda\sigma
f-\sigma\right).$ (4.6)

If $\lambda=2\sigma$, we can estimate the spectral index by the ratio
$n_{s}\approx-3+4\sigma f+\frac{\sigma f}{1-\sigma f}$. If $f=0$ we obtain
$n_{s}=-3$. From (3.7) and (4.6) we get

$\Omega_{GW}\left(f\right)=\Omega_{GW\ast}\exp N(-3+4\sigma f-\frac{\sigma
f}{1-\sigma f}\left(2\sigma f-\lambda f-1\right)),$ (4.7)

where $\Omega_{GW\ast}$ is the amplitude at the typical frequency $f_{\ast}$.
The amplitude at $f=0$ is given by
$\Omega_{GW}\left(f\right)=\Omega_{GW\ast}e^{-6N\text{ }}$.

##  5 Stochastic GW signals at preheating

Thereafter, we want to study the spectral index from the relationship (3.10)
instead of (4.1):
$H\left(f\right)=H_{end}\left(\frac{f}{f_{\ast}}\right)^{\epsilon_{H}}.$ We
calculate the derivative of this parameter with respect to time (see Appendix
A), and we determine

$\frac{\dot{H}\left(f\right)}{H\left(f\right)}=\dot{\epsilon}_{H}\ln\left(\frac{f}{f_{\ast}}\right)+\epsilon_{H}\frac{\dot{f}}{f}.$
(5.1)

Using $\epsilon_{H}=-\frac{\dot{H}}{H^{2}}$, we can relate the slow roll
parameter $\epsilon_{H}$ to GW frequency as

$-\frac{H\epsilon_{H}}{\dot{\epsilon}_{H}}=\ln\left(\frac{f}{f_{\ast}}\right)+\frac{\epsilon_{H}\dot{f}}{\dot{\epsilon}_{H}f},$
(5.2)

and hence it can be written as

$\epsilon_{H}-\frac{\dot{\epsilon}_{H}}{2H\epsilon_{H}}=\epsilon_{H}+\frac{1}{2}\left[\ln\left(\frac{f}{f_{\ast}}\right)+\frac{\epsilon_{H}\dot{f}}{\dot{\epsilon}_{H}f}\right]^{-1}.$
(5.3)

Therefore, the slow-roll parameter
$\eta=\epsilon_{H}-\frac{\dot{\epsilon}_{H}}{2H\epsilon_{H}}$ [41] is given by
$\eta=\epsilon_{H}+\frac{1}{2}\left[\frac{\epsilon_{H}\dot{f}}{\dot{\epsilon}_{H}f}+\ln\left(\frac{f}{f_{\ast}}\right)\right]^{-1}$.
Thus, the spectral index $n_{s}=1-6\epsilon_{H}+2\eta$ can be written as

$n_{s}=1-4\epsilon_{H}+\left[\frac{\epsilon_{H}\dot{f}}{\dot{\epsilon}_{H}f}+\ln\left(\frac{f}{f_{\ast}}\right)\right]^{-1}.$
(5.4)

If $\sigma f=1$, one obtain $n_{s}=1+\ln^{-1}\left(\frac{f}{f_{\ast}}\right)$,
i.e.
$\Omega_{GW}\left(f\right)=\Omega_{GW\ast}e^{\left(1+\ln^{-1}\left(\frac{f}{f_{\ast}}\right)\right)N\text{
}}$. Next, we consider $\sigma=\frac{\ddot{f}}{\dot{f}^{2}}$ to be a constant
and we obtain $\epsilon_{H}=1-\sigma f$, this yields
$\frac{\dot{\epsilon}_{H}}{\epsilon_{H}}=\frac{-\sigma\dot{f}}{1-\sigma f}$.
The $n_{s}$ with respect to $f$ is given by

$n_{s}=-3+4\sigma f+\left[-\frac{1-\sigma f}{\sigma
f}+\ln\left(\frac{f}{f_{\ast}}\right)\right]^{-1},$ (5.5)

and the energy density spectrum can be written as

$\Omega_{GW}\left(f\right)=\Omega_{GW\ast}\exp N\left(-3+4\sigma
f+\left[-\frac{1-\sigma f}{\sigma
f}+\ln\left(\frac{f}{f_{\ast}}\right)\right]^{-1}\right).$ (5.6)

In Figure (1) we plot the $n_{s}$ as a function of the frequency $f$. In this
illustration, it’s evident that the curves representing the function
$n_{s}(f)$ intersect with the range of observed values. While, in Figures
(2)-(3) we plot the energy density spectrum $\Omega_{GW}h^{2}$ as it varies
with frequency $f$.

Figure 1:  Evolution of model parameters $n_{s}$ (5.5) for $\sigma=5\times
10^{7}$ and $nHz\leq f_{\ast}\leq 2.5nHz.$ Figure 2: Evolution of energy
density spectrum $\Omega_{GW}h^{2}$ (5.6) for $\sigma=5\times 10^{7}$ and
$N=1$ with $nHz\leq f_{\ast}\leq 2.5nHz.$ Figure 3: Region of the
gravitational wave spectrum for $\sigma=10^{7}$, $\Omega_{GW\ast}\sim
10^{-14}$, $N=2$ and $f_{\ast}=0.5nHz.$

In Figure (3), we observe two distinct regions delineating the GW frequency
characteristics of the early universe. The interval $0nHz\leq f\leq 20nHz$
corresponds to the outcomes associated with the CMB. Notably, the interval
$f\sim 25nHz-65nHz$ characterizes the dataset generated by the NANOGrav 15-yr
result at the (near) end of inflation $N\sim 2$. These results directly match
the observed data [1, 2, 36]. This constraint on $f$ has significant
implications for the energy density spectrum of gravitational waves,
$\Omega_{GW}h^{2}$, at PTA frequencies. A vanishing GW frequency, which would
imply a nearly scale-invariant GW spectrum, corresponds to the first instance
of the universe. The study suggests that there exists an intermediate range of
frequencies between those of CMB and PTA, which lacks energy density spectrum.
components.
We consider that $\sigma f\ll 1$ and $-3+4\sigma f-\frac{\sigma f}{1-\sigma
f}=0$, thus, (5.5) can be written as $n_{s}\approx-\left(\frac{\sigma
f}{1-\sigma f}\right)^{2}\ln\left(\frac{f}{f_{\ast}}\right)$. The energy
density spectrum is estimated as

$\Omega_{GW}\left(f\right)h^{2}=\Omega_{GW\ast}h^{2}\left(\frac{f}{f_{\ast}}\right)^{-\left(\frac{\sigma
f}{1-\sigma f}\right)^{2}N\text{ }}.$ (5.7)

If we use $k=2\pi f$ with $k_{\ast}=2\pi f_{\ast}$ is the characteristic
scale, we find that
$\Omega_{GW}\left(f\right)=\Omega_{GW\ast}\exp\left(-\ln\left(\frac{k}{k_{\ast}}\right)\left(\frac{\sigma
f}{1-\sigma f}\right)^{2}N\right)$, which exactly corresponds to the
primordial power spectrum of curvature perturbations
$\mathcal{P}_{\mathcal{R}}\left(f\right)$ [37, 38], also it checked equation
(3.1). The relation (5.7) is in good agreement with (2.6) and the NANOGrav’s
observations, suggesting that $\Omega_{GW}\propto f^{\zeta}$ with
$\zeta\in\left(-1.5,0.5\right)$ at a frequency of $f=5.5nHz$. Combining this
with (5.7) we get

$\zeta=-\left(\frac{\sigma f}{1-\sigma f}\right)^{2}N<0.$ (5.8)

In the scenario where $\sigma\sim 10^{7}$, as illustrated in Figures (2)-(3),
we determine that $\zeta=-0.143876N<0$. This conclusion subsequently implies
the following condition: $0\leq N\leq 10.425$. Furthermore, this range ensures
that $\zeta$ remains greater than or equal to $-1.5$. In the tables below (1),
we analyze numerical representations of the relationship (5.8).

$N$ | $\zeta$ | $\sigma\left[s\right]$
---|---|---
$60$ | $-1.5$ | $2.482\times 10^{7}$
$60$ | $-0.5$ | $1.520\times 10^{7}$
$60$ | $-10^{-6}$ | $2.346\times 10^{4}$
$50$ | $-1.5$ | $2.684\times 10^{7}$
$50$ | $-0.5$ | $1.652\times 10^{7}$
$2$ | $-0.5$ | $6.060\times 10^{7}$
Table 1: The stochastic gravitational wave background based on the value of
the Hubble constant for $f=5.5nHz$.

Based on these results, it is worth noting that the values of $\zeta$
precisely satisfy the data conditions indicating that the relationship (3.10)
aligns well with the observational data. The range $0\leq N\leq 10.425$
signifies the epoch marking the end of inflation. During this period, there is
the emission of stochastic gravitational waves, which are then detected by the
NANOGrav collaboration. From the relationship (5.8) we deduce
$N\left(f\right)=\left(\frac{1-\sigma f}{\sigma
f}\right)^{2}\left|\zeta\right|$ and
$f=\frac{1}{\sigma\left(1+\sqrt{\frac{N}{\left|\zeta\right|}}\right)}$. From
this, we plot the curve of $N$ as a function of $f$ in figure (4).

Figure 4: $N$ curve as a function of $f$, such as $\sigma=1.5\times 10^{7}s$
and $\zeta=-0.5.$

In Figure (4), we observe that when considering the CMB with $30\leq N\leq 60$
[42], the frequency of the GW signal $f$ remains consistently close to the
value $f=0.08nHz$. However, as we approach the end of inflation, the range of
$N\leq 1$ values encompasses two distinct frequency values. Finally, at the
precise end of inflation ($N=0$), the frequency settles at a singular value,
which is equal to $f=1nHz$. The sharp increase in frequency observed towards
the end of the inflationary period indicates a significant cosmic variation.
This change is due to the curvature perturbation [19].
The state of the universe at the the end of inflation is known as the
preheating [18]. Typically, following the end of inflation, our Universe will
transition into a radiation-dominated era, during which perturbations will
once again enter the observable horizon.

##  6 Conclusion

In summary, the study has explored the stochastic gravitational wave produced
during the end of inflation. In this context, it has been established that the
Hubble parameter is linked to the logarithmic variation of GW frequency. By
analyzing the frequencies detected by the NANOGrav collaboration in relation
to the end of inflation marked by $H_{end}$, we have unveiled a significant
connection. This connection suggests a logarithmic relationship between the
Hubble parameter and gravitational wave frequency, represented as $f\propto
a$, where $a$ represents the scale factor. This proportional link between
scale factor and frequency could shed light on the Hubble tension problem and
its connection to the NANOGrav 15-yr results. Furthermore, we present plots of
the spectral index $n_{s}$ as a function of frequency $f$, showing how
$n_{s}\left(f\right)$ intersects with observed values. We have also explored
the energy density spectrum $\Omega_{GW}h^{2}$ in relation to frequency,
revealing two distinct regions in Figure (3), corresponding to GW frequencies
associated with the early universe and those generated by NANOGrav 15-yr
results at the end of inflation. These results align with observed data,
emphasizing the importance of frequency constraints for the energy density
spectrum of gravitational waves at PTA frequencies. Notably, an absence of
energy density spectrum components is identified in an intermediate frequency
range, bridging the CMB and PTA domains. In Figure (4), we have observed the
behavior of GW signal frequency concerning the CMB for different inflationary
parameters. Remarkably, as we approach the end of inflation ($N\leq 1$), two
distinct frequency values emerge. At the precise end of inflation ($N=0$), the
frequency stabilizes at a singular value, $f=1nHz$. This study proposes that
primordial universe scenarios, which rely on the emission of the stochastic
gravitational waves at the end of inflation, can offer a new period of the
emission of this stochastic gravitational wave background. This could result
in a greater likelihood of observing distinct signals in PTA observations.

##  7 Appendix A

1- Inflationary scenario from e-folding number
From $\dot{N}=-H$ we have $\dot{N}=-\frac{\dot{f}}{f}=-H$, which leads to
$H=\frac{\dot{f}}{f}$. From this we get

$\dot{H}=\frac{\ddot{f}f-\dot{f}^{2}}{f^{2}}.$ (7.1)

We employ the constraint
$\epsilon_{H}=-\frac{\dot{H}}{H^{2}}=\frac{\dot{f}^{2}-\ddot{f}f}{\dot{f}^{2}}$derived
from $\epsilon_{H}=-\frac{\dot{H}}{H^{2}}$ and we find

$\epsilon_{H}=1-\frac{\ddot{f}f}{\dot{f}^{2}}=1-\sigma f,$ (7.2)

where $\sigma=\frac{\ddot{f}}{\dot{f}^{2}}$ and
$\lambda=\frac{\dddot{f}}{\ddot{f}}\frac{1}{\dot{f}}$. Using
$2H\epsilon_{H}=2\frac{\dot{f}^{2}-\ddot{f}f}{\dot{f}f}$ and
$\frac{1}{2H\epsilon_{H}}=\frac{\dot{f}f}{2\left(\dot{f}^{2}-\ddot{f}f\right)}$
we can calculate the slow-roll parameter
$\eta=\epsilon_{H}-\frac{\dot{\epsilon}_{H}}{2H\epsilon_{H}}$. We have

$\dot{\epsilon}_{H}=2\left(\frac{\ddot{f}}{\dot{f}}\right)^{2}\frac{f}{\dot{f}}-\frac{\dddot{f}}{\dot{f}^{2}}f-\frac{\ddot{f}}{\dot{f}},$
(7.3)

which leads to

$\eta=1-\frac{\ddot{f}f}{\dot{f}^{2}}-\frac{f}{2\left(1-\frac{\ddot{f}f}{\dot{f}^{2}}\right)}\left(2\left(\frac{\ddot{f}}{\dot{f}}\right)^{2}\frac{f}{\dot{f}^{2}}-\frac{\dddot{f}}{\dot{f}}\frac{f}{\dot{f}^{2}}-\frac{\ddot{f}}{\dot{f}^{2}}\right).$
(7.4)

So we can relate the slow-roll parameter to frequency of the GW signal $f$ as

$\eta=1-\sigma f-\frac{f}{2\left(1-\sigma
f\right)}\left(2\sigma^{2}f-\lambda\sigma f-\sigma\right).$ (7.5)

We can additionally calculate the spectral index using the following approach:
the spectral index $n_{s}=1-6\epsilon_{H}+2\eta$ is

$n_{s}=1-6+6\sigma f+2-2\sigma f-\frac{f}{\left(1-\sigma
f\right)}\left(2\sigma^{2}f-\lambda\sigma f-\sigma\right).$ (7.6)

So, we find equation (4.5):

$n_{s}=-3+4\sigma f-\frac{f}{\left(1-\sigma
f\right)}\left(2\sigma^{2}f-\lambda\sigma f-\sigma\right).$ (7.7)

2- Inflationary scenario from Hubble parameter
From the relationship (3.10)
$H\left(f\right)=H_{end}\left(\frac{f}{f_{\ast}}\right)^{\epsilon_{H}}$ we
have

$\dot{H}\left(f\right)=\frac{d}{dt}\left(\epsilon_{H}\ln\left(\frac{f}{f_{\ast}}\right)\right)H\left(f\right),$
(7.8)

where

$\frac{d}{dt}\left(\epsilon_{H}\ln\left(\frac{f}{f_{\ast}}\right)\right)=\dot{\epsilon}_{H}\ln\left(\frac{f}{f_{\ast}}\right)+\epsilon_{H}\frac{\dot{f}}{f},$

which leads to

$\frac{\dot{H}\left(f\right)}{H\left(f\right)}=\left(\dot{\epsilon}_{H}\ln\left(\frac{f}{f_{\ast}}\right)+\epsilon_{H}\frac{\dot{f}}{f}\right).$
(7.9)

Using $\epsilon_{H}=-\frac{\dot{H}}{H^{2}}$, so we find equation (5.1)
$-H\epsilon_{H}=\dot{\epsilon}_{H}\ln\left(\frac{f}{f_{\ast}}\right)+\epsilon_{H}\frac{\dot{f}}{f}.$

## References

  * [1] Afzal, A., et al. and NANOGrav Collaboration, 2023, The Astrophysical Journal Letters, 951(1), L11.
  * [2] G. Agazie et al. and NANOGrav Collaboration, 2023, The Astrophysical Journal Letters, 951(1), L8.
  * [3] Arzoumanian, Z., et al. and NANOGrav Collaboration, 2020, The Astrophysical Journal Letters, 905(2), L34.
  * [4] Arzoumanian, Z., et al. and NANOGrav Collaboration, 2021, Physical review letters, 127(25), 251302.
  * [5] Antoniadis, J., et al., 2022, Monthly Notices of the Royal Astronomical Society, 510(4), 4873-4887.
  * [6] Middleton, H., et al., 2021, Monthly Notices of the Royal Astronomical Society: Letters, 502(1), L99-L103.
  * [7] Zhang, Z., Cai, C., Su, Y. H., Wang, S., Yu, Z. H., and Zhang, H. H., 2023, arXiv preprint arXiv:2307.11495.
  * [8] Vagnozzi, S., 2023, Journal of High Energy Astrophysics, 39, 81-98.
  * [9] An, H., Su, B., Tai, H., Wang, L. T., and Yang, C., 2023, arXiv preprint arXiv:2308.00070.
  * [10] Niu, X., and Rahat, M. H., 2023, arXiv preprint arXiv:2307.01192.
  * [11] Mansoori, S. A. H., Felegray, F., Talebian, A., and Sami, M., 2023, Journal of Cosmology and Astroparticle Physics, 2023(08), 067.
  * [12] Cheung, K., Ouseph, C. J., and Tseng, P. Y., 2023, arXiv preprint arXiv:2307.08046.
  * [13] Bousder, M., Riadsolh, A., Fatimy, A. E., Belkacemi, M. E., and Ez-Zahraouy, H., 2023, arXiv preprint arXiv:2307.10940.
  * [14] Bousder, M. (2022). A new constant behind the rotational velocity of galaxies. Journal of Cosmology and Astroparticle Physics, 2022(01), 015.
  * [15] Bousder, M., Salmani, E., & Ez-Zahraouy, H. (2023). Entropy as logarithmic term of the central charge and modified Friedmann equation in AdS/CFT correspondence. Journal of High Energy Astrophysics, 38, 49-57.
  * [16] Bousder, M., & Bennai, M. (2021). Particle-antiparticle in 4D charged Einstein-Gauss-Bonnet black hole. Physics Letters B, 817, 136343.
  * [17] Ellis, J., Lewicki, M., Lin, C., and Vaskonen, V., 2023, arXiv preprint arXiv:2306.17147.
  * [18] Easther, R., Giblin Jr, J. T., and Lim, E. A., 2007, Physical Review Letters, 99(22), 221301.
  * [19] Lyth, D. H., 2005, Journal of Cosmology and Astroparticle Physics, 2005(11), 006.
  * [20] Riess, A. G., et al., 2016, The Astrophysical Journal, 826(1), 56\.
  * [21] Aghanim, N., et al., 2020, Astronomy & Astrophysics, 641, A6.
  * [22] Di Valentino, et al., 2021, Classical and Quantum Gravity, 38(15), 153001.
  * [23] Vagnozzi, S., 2023, Universe, 9(9), 393.
  * [24] Schiavone, T., Montani, G., and Bombacigno, F., 2023, Monthly Notices of the Royal Astronomical Society: Letters, 522(1), L72-L77.
  * [25] Yang, Y., Lu, X., Qian, L., & Cao, S., 2023, Monthly Notices of the Royal Astronomical Society, 519(4), 4938-4950.
  * [26] Monjo, R., & Campoamor-Stursberg, R., 2023, Classical and Quantum Gravity.
  * [27] Poulin, V., Smith, T. L., Karwal, T., and Kamionkowski, M., 2019, Physical review letters, 122(22), 221301.
  * [28] Jedamzik, K., and Pogosian, L.,2020, Physical Review Letters, 125(18), 181302.
  * [29] Christensen, N.,2018, Reports on Progress in Physics, 82(1), 016903..
  * [30] Kuroyanagi, S., Chiba, T., and Takahashi, T., 2018, Journal of Cosmology and Astroparticle Physics, 2018(11), 038..
  * [31] Zhao, W., Zhang, Y., You, X. P., and Zhu, Z. H.,2013, Physical Review D, 87(12), 124012.
  * [32] Liu, X. J., Zhao, W., Zhang, Y., and Zhu, Z. H., 2016, Physical Review D, 93(2), 024031.
  * [33] Vagnozzi, S.,2021, Monthly Notices of the Royal Astronomical Society: Letters, 502(1), L11-L15.
  * [34] Bian, L., Ge, S., Li, C., Shu, J., and Zong, J., 2022, arXiv preprint arXiv:2212.07871.
  * [35] Domcke, V., Pieroni, M., and Binétruy, P. ,2016, Journal of Cosmology and Astroparticle Physics, 2016(06), 031.
  * [36] Niu, X., & Rahat, M. H., 2023, arXiv preprint arXiv:2307.01192.
  * [37] Liu, L., Chen, Z. C., and Huang, Q. G., 2023, arXiv preprint arXiv:2307.14911.
  * [38] Chang, Z., Kuang, Y. T., Wu, D., and Zhou, J. Z., 2023, arXiv preprint arXiv:2309.06676.
  * [39] Bousder, M., Riadsolh, A., El Belkacemi, M. and Ez-Zahraouy, H., 2023, Annals of Physics, 458, 169441.
  * [40] Bousder, M., Salmani, E., El Fatimy, A., & Ez-Zahraouy, H., 2023, Annals of Physics, 452, 169282.
  * [41] Maartens, R., Wands, D., Bassett, B. A., and Heard, I. P., 2000, Physical Review D, 62(4), 041301.
  * [42] Giarè, W., Forconi, M., Di Valentino, E., & Melchiorri, A., 2023, Monthly Notices of the Royal Astronomical Society, 520(2), 1757-1773.