# Very high-energy gamma rays from GRB 180720B and GRB 190829A with external
Compton emission

Monica Barnard,1 Soebur Razzaque,1,2,3 and Jagdish C. Joshi 4,1
1Centre for Astro-Particle Physics (CAPP), Department of Physics, University
of Johannesburg, PO Box 524, Auckland Park 2006, South Africa
2Department of Physics, The George Washington University, Washington, DC
20052, USA
3 National Institute for Theoretical and Computational Sciences (NITheCS),
South Africa
4 Aryabhatta Research Institute of Observational Sciences (ARIES), Manora
Peak, Nainital 263001, India
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Gamma-ray bursts (GRBs) comprise of short, bright, energetic flashes of
emission from extragalactic sources followed by a longer afterglow phase of
decreased brightness. Recent discoveries of very-high-energy (VHE, $\gtrsim
100$ GeV) afterglow emission from GRB 180720B and GRB 190829A by H.E.S.S. have
raised questions regarding the emission mechanism responsible. We interpret
these observed late-time emission to be the result of inverse Compton emission
of ultra-relativistic electrons in the GRB blastwave in an external radiation
field, i.e., external Compton (EC), considering both the wind and interstellar
medium scenarios. We present predictions of multiwavelength light curves and
energy spectra, ranging from optical to VHE, and include the synchrotron and
synchrotron self-Compton (SSC) radiation mechanisms as well. We corrected the
EC and SSC model for the $\gamma$-ray attenuation by absorption of photons
through their interaction with the extragalactic background light (EBL). We
compared our results to multiwavelength data and found that EC gives a
satisfactory fit for a given set of fixed model parameters for GRB 180720B,
whereas SSC result in a better fit for GRB 190829A. For both GRBs a wind
environment is preferred over constant density inter-stellar medium, and the
Cosmic Microwave Background as the external radiation field. However, with
more data and an effective optimisation tool we can find a more robust fit of
the model, implying better constraints on the GRB environment and the particle
energy requirements for the emission observed at late times. This has
consequences for future observations of GRBs at these extreme energies.

###### keywords:

(transients:) gamma-ray bursts – (stars:) gamma-ray burst: individual: 180720B
– (stars:) gamma-ray burst: individual: 190829A

††pubyear: 2023††pagerange: Very high-energy gamma rays from GRB 180720B and
GRB 190829A with external Compton emission–Very high-energy gamma rays from
GRB 180720B and GRB 190829A with external Compton emission

##  1 Introduction

Gamma-ray burst (GRB) afterglow radiation is characterised by their diverse
emission across broad wavelengths, making them fascinating multiwavelength
emitters (see, for reviews, Piran, 2004; Mészáros, 2006; Kumar & Zhang, 2015).
Their multiwavelength features can be explained using the external-shock
models by including detail processes, i.e., reverse shock, structured jets,
arbitrary circumburst medium density profile, etc (Wang et al., 2015).
Synchrotron radiation in the external shocks is very useful to explain the
radiation up to $\sim 10$ GeV (see, e.g., Barniol Duran & Kumar, 2011; Gehrels
& Razzaque, 2013; Gill & Granot, 2022). However, processes such as the
inverse-Compton scattering (see, e.g., Chiang & Dermer, 1999; Sari & Esin,
2001; Zhang & Mészáros, 2001; Fan et al., 2008), proton-synchrotron radiation
(see, e.g., Razzaque et al., 2010; Ghisellini et al., 2020; Isravel et al.,
2022; Zhang et al., 2023) and photopion interactions (see, e.g., Böttcher &
Dermer, 1998; Asano et al., 2009; Sahu & Fortín, 2020) can extend the spectrum
of gamma-rays beyond 10 GeV, as expected in theoretical models. However, the
confirmation of VHE signal in afterglows came recently by the ground-based
Cherenkov telescopes; i.e., Major Atmospheric Gamma Imaging Cherenkov
Telescopes (MAGIC) (MAGIC Collaboration et al., 2019b) and High Energy
Stereoscopic System (H.E.S.S.) (Abdalla et al., 2021), respectively. The VHE
signal from GRB 180720B (Abdalla et al., 2019), GRB 190829A (Abdalla et al.,
2021), and GRB 190114C (MAGIC Collaboration et al., 2019b, a) have provided us
with an opportunity to investigate the physical properties of these GRBs.

The VHE signal has been modelled using synchrotron-self-Compton (SSC)
mechanism for GRB 190114C by MAGIC Collaboration et al. (2019b); Wang et al.
(2019); Joshi & Razzaque (2021); Derishev & Piran (2021); Yamasaki & Piran
(2022); Klinger et al. (2023). However, hadronic $p\gamma$ interaction models
are also used to study this burst (Sahu & Fortín, 2020). The VHE observations
are also used to constrain the intergalactic magnetic field values (Wang et
al., 2020; Da Vela et al., 2023) and the energetics and radiative efficiency
of GRB fireball (Li & Wang, 2023). SSC and $p\gamma$ processes have also been
used to study the broadband radiation from GRB 180720B (Fraija et al., 2019;
Sahu & Fortín, 2020) and GRB 190829A (Abdalla et al., 2021; Fraija et al.,
2021; Zhang et al., 2021a; Sahu et al., 2022). Further, Huang et al. (2022a)
used the VHE afterglow emission to study the maximum energy of electrons
accelerated in those blastwaves. Recently, an unusually bright burst GRB
221009A, also termed as the brightest of all time (B.O.A.T.), has been
discovered (Dichiara et al., 2022; Galanti et al., 2022; Huang et al., 2022b;
Veres et al., 2022) and its afterglow has been detected from radio to VHE
energies (Williams et al., 2023; LHAASO Collaboration et al., 2023; Lesage et
al., 2023). This burst was also followed by the IceCube Neutrino Observatory
but not detected (Abbasi et al., 2023). Further, based on the multiwavelength
observations, detection of $\sim 18$ TeV photons (Huang et al., 2022b) has
been used to explore new physics (Baktash et al., 2022; Finke & Razzaque,
2023; Nakagawa et al., 2023; Li & Ma, 2023; Cheung, 2022; Smirnov & Trautner,
2022), GRB jet composition studies (Das & Razzaque, 2023; Alves Batista, 2022;
Liu et al., 2023; Laskar et al., 2023; Zhang et al., 2023), and jet structure
(Sato et al., 2023; Ren et al., 2023) and jitter radiation as a possible
mechanism for its production (Mao & Wang, 2023) etc.

In this paper, we will focus on the origin of the VHE signal from GRB 180720B
(Abdalla et al., 2019) and GRB 190829A (Abdalla et al., 2021) and study the
role played by the external-Compton (EC) radiation mechanism and fit
multiwavelength data. However, we will also include the contributions due to
the synchrotron and SSC mechanism. In the EC scenarios, relativistic electrons
in the blastwave upscatter photons from an external radiation field.
Previously, the EC model has been applied for the VHE radiation from GRB
190829A by Zhang et al. (2021b), where they used _Swift_ -XRT detected flare
photons at $\sim 10^{3}$ s, as target photons. We considered the low energy
seed photons to be from Cosmic Microwave Background (CMB) origin. In Section 2
we briefly describe the blastwave evolution, particle transport, the different
radiation mechanisms considered and the implementation thereof, as well as the
absorption of sub-TeV photons in the blastwave. In Section 3 we apply our
model to the before mentioned GRBs and compare our results to the multi-
wavelength observations, with the discussion following in Section 4.

##  2 Multi-wavelength emission model

###  2.1 GRB blastwave evolution

The GRB blastwave with isotropic-equivalent initial kinetic energy of $E_{k}$
(erg) and initial bulk Lorentz factor $\Gamma$ interacts with circumburst
medium of density $n(R)$ and generates external shocks (Blandford & McKee,
1976). Here $R$ is the distance from the central engine. The blastwave
accumulates matter from the surrounding environment and decelerates. Here we
have only used the adiabatic expansion scenario in a constant-density
interstellar medium (ISM) and in a wind environment. We have also used
formulas and their numerical values from Joshi & Razzaque (2021). The
blastwave deceleration time is then

$\displaystyle t_{\rm dec}=\left[\frac{17E_{k}(1+z)^{3}}{64\pi
nm_{p}c^{5}\Gamma^{8}}\right]^{1/3}\,;~{}~{}~{}~{}{\rm ISM}$ (1)

in ISM, where $n(R)=n_{0}~{}{\rm cm^{-3}}$, and

$\displaystyle t_{\rm dec}=\frac{9E_{k}(1+z)}{16\pi
Am_{p}c^{3}\Gamma^{4}}\,;~{}~{}~{}~{}{\rm Wind}$ (2)

in wind, where $n(R)=AR^{-2}~{}{\rm cm^{-3}}$. For a stellar wind of velocity
$v_{w}=10^{8}v_{8}~{}{\rm cm~{}s^{-1}}$ with mass loss rate
$\dot{M}_{w}=10^{-5}\dot{M}_{-5}M_{\odot}~{}{\rm yr^{-1}}$, we have
$A=\dot{M}_{w}/(4\pi v_{w}m_{p})=3.02\times 10^{35}A_{\star}~{}{\rm cm^{-1}}$,
and $A_{\star}\equiv\dot{M}_{-5}/v_{8}$. The time evolution of the Lorentz
factor of the shocked gas in the blastwave can be written in the ISM and wind
scenario as (see, e.g., Razzaque, 2013)

$\displaystyle\Gamma_{\rm
g}(t)=\frac{\Gamma}{\zeta}\begin{cases}\left(\frac{t_{\rm
dec}}{t}\right)^{3/8}\,;&{\rm ISM}\\\ \left(\frac{t_{\rm
dec}}{t}\right)^{1/4}\,;&{\rm Wind}\end{cases}$ (3)

where the parameter $\zeta$ depends on the details of the blastwave evolution,
in particular on the similarity parameter (Blandford & McKee, 1976). For our
calculations, we have used $\zeta=\sqrt{2}$. The corresponding radius of the
blastwave evolves as

$\displaystyle R(t)=\begin{cases}\frac{16\Gamma^{2}_{\rm
g}(t)ct}{(1+z)}\,;&{\rm ISM}\\\ \frac{8\Gamma^{2}_{\rm g}(t)ct}{(1+z)}\,;&{\rm
Wind}\end{cases}$ (4)

###  2.2 Particle acceleration and spectrum

A population of non-thermal electrons in the shocked gas is assumed to be
accelerated in the magnetic field in the shocked region. The energy density of
this magnetic field $u^{\prime}_{B}=B^{\prime 2}/8\pi$ is assumed to be a
fraction $\epsilon_{B}$ of the energy density of the shocked gas in the
blastwave frame, $u^{\prime}_{\rm g}=4n(R)m_{p}c^{2}\Gamma_{\rm g}^{2}$, for a
strong shock. Thus the magnetic field is

$\displaystyle B^{\prime}=\sqrt{32\pi\epsilon_{B}n(R)m_{p}c^{2}}\Gamma_{\rm
g}.$ (5)

The electrons also cool due to synchrotron and SSC radiation. The resulting
electron Lorentz factor spectrum in the slow-cooling regime, isotropic in the
blastwave frame, relevant for late afterglow emission, is given by

$\displaystyle
n^{\prime}(\gamma^{\prime})=\frac{k}{4\pi}\times\begin{cases}(\gamma^{\prime}/\gamma_{c}^{\prime})^{-p};~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\gamma^{\prime}_{m}\leq\gamma^{\prime}\leq\gamma^{\prime}_{c}\\\
(\gamma^{\prime}/\gamma_{c}^{\prime})^{-p-1};~{}~{}~{}~{}~{}~{}~{}\gamma^{\prime}_{c}<\gamma^{\prime}\leq\gamma^{\prime}_{s}\end{cases},$
(6)

with $\gamma^{\prime}$ the Lorentz factor in the blastwave frame, and $p$ the
electron index. Here

$\displaystyle\gamma^{\prime}_{m}(t)=\left[\frac{m_{p}}{m_{e}}\epsilon_{e}\frac{p-2}{p-1}\Gamma_{\rm
g}\right];\,p>2,$ (7)

is the minimum Lorentz factor, assuming a fraction $\epsilon_{e}$ of
$u^{\prime}_{g}$ is acquired by the shock-accelerated electrons. The the
cooling Lorentz factor is given by

$\displaystyle\gamma^{\prime}_{c}(t)=\left[\frac{6\pi
m_{e}c^{2}(1+z)}{\sigma_{T}cB{{}^{\prime}}^{2}(t)t\Gamma_{\rm
g}(1+Y+\Gamma_{\rm g}^{2}u_{\rm ex}/u_{B}^{\prime})}\right]$ (8)

with $u_{\rm ex}$ the energy density of the external radiation field, and is
found by equating the cooling time scale and dynamic time scale, and

$\displaystyle\gamma^{\prime}_{s}(t)=\left[\frac{6\pi
e}{\phi\sigma_{T}B^{\prime}(t)(1+Y+\Gamma_{\rm g}^{2}u_{\rm
ex}/u_{B}^{\prime})}\right]^{1/2}$ (9)

is the maximum or saturation Lorentz factor, found by equating the cooling
time scale to the acceleration time scale. In the above equations, Eq. 7 is
only valid when
$\gamma^{\prime}_{s}(t)^{(p-2)}\gg\gamma^{\prime}_{m}(t)^{(p-2)}$, $m_{p}$ and
$m_{e}$ is the proton and electron rest masses respectively, $\sigma_{T}$ the
Thomson cross section, $c$ the speed of light, and $\phi\geq 1$ is the number
of gyro-radius required for accelerating electrons to $\gamma^{\prime}_{s}$.
The Compton $Y$-parameter in the slow-cooling, Thomson regime is given by

$\displaystyle
Y=\sqrt{\epsilon_{e}/\epsilon_{B}}\times\begin{cases}(t/t_{0}^{ssc})^{(2-p)/[2(4-p)]};~{}~{}~{}{\rm
Adiabatic-ISM}\\\ (t/t_{0}^{ssc})^{(2-p)/(4-p)};~{}~{}~{}~{}~{}~{}{\rm
Adiabatic-Wind}\\\ \end{cases},$ (10)

with $t_{0}$ the transition time from the fast- to slow-cooling spectra, and
$Y\approx\sqrt{\epsilon_{e}/\epsilon_{B}}\approx 1$ (implying
$\epsilon_{e}\approx\epsilon_{B}$; see Equation (13) in Joshi & Razzaque
(2021)). Note that the second part of the equation is rather small for $p\sim
2.0-2.3$ (Sari & Esin, 2001). $Y\approx 1$ value is imposed when we
investigate the EC emission and need to reduce the SSC flux.

The normalisation constant $k$ in Equation (6) can be found from the integrals
$\int
d\Omega\int_{\gamma^{\prime}_{m}}^{\gamma^{\prime}_{s}}n^{\prime}(\gamma^{\prime})d\gamma=4n(R)\Gamma_{g}$,
where $n(R)$ is the particle density in the surrounding environment in the lab
frame, as

$\displaystyle k=\frac{4n(R)\Gamma_{\rm
g}}{\int_{\gamma^{\prime}_{m}}^{\gamma^{\prime}_{c}}(\gamma^{\prime}/\gamma_{c}^{\prime})^{-p}d\gamma^{\prime}+\int_{\gamma^{\prime}_{c}}^{\gamma^{\prime}_{s}}(\gamma^{\prime}/\gamma_{c}^{\prime})^{-p-1}d\gamma^{\prime}}.$
(11)

Using the transformation $\gamma/\gamma^{\prime}=\Gamma_{\rm g}$ and
$n(\gamma,\mu)=n^{\prime}(\gamma^{\prime})(\gamma/\gamma^{\prime})^{2}$ (with
$\gamma$ the Lorentz factor in the lab frame, $\mu=\cos{\theta}$, and $\theta$
the angle of the blastwave with respect to the observer’s line of sight), one
obtains the lab frame electron distribution as

$\displaystyle N(\gamma,\mu)$ $\displaystyle=$ $\displaystyle
n^{\prime}(\gamma^{\prime})\Gamma_{\rm g}^{3}V$ (12) $\displaystyle=$
$\displaystyle\frac{k}{4\pi}\Gamma_{g}^{3}V\times\begin{cases}(\gamma^{\prime}/\gamma_{c}^{\prime})^{-p};~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\gamma^{\prime}_{m}\leq\gamma^{\prime}\leq\gamma^{\prime}_{c}\\\
(\gamma^{\prime}/\gamma_{c}^{\prime})^{-p-1};~{}~{}~{}~{}~{}~{}~{}\gamma^{\prime}_{c}<\gamma^{\prime}\leq\gamma^{\prime}_{s}\,.\end{cases}$

Here $V=(4/3)\pi R^{3}$ is the volume of the blastwave in the lab frame.

###  2.3 External Compton radiation

Electrons can lose energy by Compton scattering photons of an external
radiation field. This external Compton (EC) process is relevant when the
energy density of the external radiation field in the jet frame exceeds that
of the energy density in the magnetic field, i.e., $u_{\rm
ex}\Gamma_{g}^{2}\gtrsim
u^{\prime}_{B}=4n(R)m_{p}c^{2}\Gamma_{g}^{2}\epsilon_{B}$. This gives a
condition on the number density of particles in the environment as
$n(R)\lesssim u_{\rm ex}/(4m_{p}c^{2}\epsilon_{B})$ for EC to be important.
For CMB with energy density $4.17\times 10^{-13}(1+z)^{4}$ erg/cm3, the
requirement is then

$n(R)\lesssim 7\times 10^{-5}(1+z)^{4}\epsilon_{B}^{-1}~{}{\rm cm}^{-3}\,.$
(13)

For a wind-type environment this condition leads to a time scale of

$t_{\rm EC,wind}\gtrsim 5.7\times
10^{5}\,\epsilon_{B,-2}A_{\star,-2}^{2}E_{54}^{-1}\,{\rm s},$ (14)

when VHE emission by EC is important. Here
$\epsilon_{B,-2}=\epsilon_{B}/0.01$, $A_{\star,-2}=A_{\star}/0.01$ and
$E_{54}=E_{k}/10^{54}$ erg.

To calculate the EC spectrum we followed the same numerical calculations as
given in Georganopoulos et al. (2001), to study the blastwave evolution
(Section 2.1) and the particle transport (Section 2.2) of GRB afterglow
emission. This model assumes an isotropic distribution of soft mono-energetic
photons upscattered by non-thermal electrons described by a power-law density
distribution. We adapted this model by describing the electron density as a
broken power-law (see Equation (6)), assuming that non-thermal ultra-
relativistic electrons in the blastwave upscatter low-energy photons from an
external radiation field via IC.

In the Klein-Nishina (KN) and Thomson limits we obtain the specific luminosity
$L$ per energy $\epsilon$ interval per solid angle $\Omega$ in the same manner
as Georganopoulos et al. (2001), using the electron distribution in Equation
(12), scattering rate of the photons
$dN_{p}/dtd\epsilon=3\sigma_{T}cf(x)/4\epsilon_{0}\gamma^{2}$, observed photon
energy $\epsilon$ (in units of $m_{e}c^{2}$), the photon number density
$n_{p}={u_{\rm ex}}/\epsilon_{0}m_{e}c^{2}$, the photon energy density $u_{\rm
ex}$ and the target photon energy $\epsilon_{0}$ (of the external radiation
field), the expression yields

$\begin{split}\epsilon\frac{dL}{d\epsilon
d\Omega}&=D^{3+p}\frac{3kV\sigma_{\rm T}c{u_{\rm
ex}}}{16\pi}\left(\frac{\epsilon}{\epsilon_{0}}\right)^{2}\\\
&\times\begin{cases}\int_{\gamma_{\rm m}}^{\gamma_{\rm
c}}\gamma^{-(2+p)}\gamma_{c}^{p}f(x)d\gamma;~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\gamma_{m}\leq\gamma\leq\gamma_{c}\\\
\int_{\gamma_{\rm c}}^{\gamma_{\rm
s}}\gamma^{-(3+p)}\gamma_{c}^{(1+p)}f(x)d\gamma;~{}~{}~{}~{}~{}~{}\gamma_{c}<\gamma\leq\gamma_{s},\end{cases}\end{split}$
(15)

in erg-s-1, and assuming $\Gamma_{\rm g}\approx D$ with $D$ being the Doppler
factor, and $k$ the normalisation factor in Equation (11). In the Klein-
Nishina case, for energies $\epsilon_{\rm min,KN}\leq\epsilon\leq\epsilon_{\rm
max,KN}$, the lower limit of integration in Equation (15) is found by setting
$x=1$, and is given by

$\displaystyle\gamma_{\rm min}$ $\displaystyle=$
$\displaystyle\frac{\epsilon\epsilon_{0}+\sqrt{\epsilon^{2}\epsilon_{0}^{2}+\epsilon\epsilon_{0}}}{2\epsilon_{0}},$
(16)

when $\gamma_{\rm min}>\gamma_{m}$ we used $\gamma_{\rm min}$, and vice versa.

Additionally, we considered the effect of the extragalactic background light
(EBL) $\gamma$-ray attenuation on our EC model. Thus, we fit an attenuated
model of the form

$\displaystyle\frac{dN}{dE}=\left(\frac{dN}{dE}\right)_{\rm
EC}e^{-\tau(E,z)}=\left(\frac{L}{4\pi D_{\rm L}^{2}}\right)_{\rm
EC}e^{-\tau(E,z)},$ (17)

where $(dN/dE)_{\rm EC}$ is the intrinsic EC energy flux (calculated by
dividing Equation (15) by $4\pi D_{\rm L}^{2}$, with $D_{\rm L}$ the
luminosity distance of the source), the exponential term corresponding to the
attenuation, $\tau$ is the energy-dependent optical depth for a source at
redshift $z$.

###  2.4 Synchrotron and Synchrotron self-Compton radiation

In addition to the external Compton described above electrons will also emit
radiation due to the synchrotron and SSC processes. To estimate these
radiation spectra we have followed the formalism as discussed in Sari et al.
(1998); Sari & Esin (2001); Joshi & Razzaque (2021). For the GRBs discussed
here, the late emission occurs in the slow cooling regime for which the
synchrotron spectrum can be calculated using

$\displaystyle F_{\nu,\rm{slow}}$
$\displaystyle=f_{\nu,\rm{max}}\begin{cases}(\frac{\nu}{\nu_{m}})^{1/3};~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\nu\leq\nu_{m}\\\
(\frac{\nu}{\nu_{m}})^{-(p-1)/2};~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\nu_{m}<\nu<\nu_{c}\\\
(\frac{\nu_{c}}{\nu_{m}})^{-(p-1)/2}(\frac{\nu}{\nu_{c}})^{-p/2};~{}~{}~{}~{}~{}~{}\nu\geq\nu_{c}\,.\end{cases}$
(18)

Here $\nu$ is the synchrotron radiation frequency and $\nu_{m}$ is the
synchrotron frequency corresponding to the minimum electron Lorentz factor
$\gamma^{\prime}_{m}$, $\nu_{c}$ corresponds to the cooling Lorentz factor
$\gamma^{\prime}_{c}$ and $\nu_{s}$ corresponds to the maximum value of the
electron Lorentz factor $\gamma^{\prime}_{s}$. The value of the maximum
synchrotron flux density $f_{\nu,\rm{max}}$ is taken from Razzaque (2013).
During the SSC scattering the synchrotron photons are upscattered by the non-
thermal electrons and in the Thomson regime the upscattered frequencies are
$\nu_{m}^{\rm ssc}\simeq 2\gamma_{m}^{\prime 2}\nu_{a},\nu_{c}^{\rm ssc}\simeq
2\gamma_{c}^{\prime 2}\nu_{c}$. The SSC spectra have been calculated using the
smooth approximation, where for each electron energy the possible interactions
with the synchrotron photons are integrated, as discussed by Sari & Esin
(2001). However, KN effects will start above Lorentz factor
$\gamma^{\prime}_{\rm KN}\simeq m_{e}c^{2}\Gamma_{g}/h\nu_{c}(1+z)$ (Joshi &
Razzaque, 2021) and corresponding maximum photon energy in the Thomson regime
for the wind case would be $E_{\gamma,\rm cut}^{\rm ssc,Th}=5.5~{}{\rm
TeV}A_{\star,-2}^{3/2}\epsilon_{B,-1}^{3/2}t_{2}^{-1}(1+Y)^{2}$. Above this
energy the SSC spectrum is steeper due to the KN effects as discussed in Nakar
et al. (2009) and the same is applied to our SSC spectrum.

##  3 Multiwavelength modelling of GRBs with late time VHE Emission

We model optical to VHE gamma-ray data from the afterglows of GRB 180720B and
GRB 190829A detected by the H.E.S.S. in this section. We have considered
synchrotron radiation (SR), SSC and EC from relativistic electrons accelerated
in the blastwaves of these GRBs. Attenuation of VHE photons in the EBL is
included in our modelling. We have considered both the ISM and wind
environments but the wind environment gives a better fit to data. We also
considered infrared, optical and CMB as the external radiation field, but the
CMB gives a better fit to the VHE data. Below, we have discussed the
observational properties and simultaneous interpretation of their light curves
and SEDs.

###  3.1 GRB 180720B

Figure 1: Multiwavelength light curves including optical (purple; r-band), the
_Swift_ -BAT spectra (15 keV-150 keV) extrapolated to the XRT band (0.3–10
keV) for a combined light curve (grey; http://www.swift.ac.uk/xrt_curves), and
_H.E.S.S._ (100-440 GeV; Abdalla et al. 2019) observations, for long GRB
180720B in the wind case for the CMB target photons ($\epsilon_{0}=8.9\times
10^{-4}$ eV). We indicate $t_{\rm dec,w}$ (black dashed line), and show the
synchrotron (purple and grey lines), SSC (orange lines) and EC (darkred
lines). The dashed (Finke et al., 2010) and dotted (Inoue et al., 2013) lines
take into account the EBL attenuation at energies $>100$ GeV. See relevant
physical parameters listed in Table 1. Figure 2: Multiwavelength SED of GRB
180720B corresponding to the light curves in Figure 1. The grey box shows the
spectrum and uncertainty of the X-ray data (an extrapolation to the time
window of H.E.S.S. observation), and the purple point is the optical data
(Fraija et al., 2019). The cyan and green lines are for the joint spectral
fits of the synchrotron with SSC and EC respectively. The SSC and EC
components are corrected for EBL attenuation using the models of Finke et al.
(2010) (dashed lines) and Inoue et al. (2013) (dotted lines).

GRB 180720B, with a redshift $z$=0.653 and a $T_{90}=48.9\pm 0.4$ s, is one of
the brightest VHE sources and was detected by H.E.S.S. at energies up to 440
GeV in the late afterglow at $T_{0}+10$ hours (Abdalla et al., 2019). During
its prompt phase, the isotropic energy released in gamma-rays is $(6\pm
0.1)\times 10^{53}$ erg (Abdalla et al., 2019). The gamma-ray observations by
Fermi-LAT also infer a smooth transition from prompt to afterglow (Ronchi et
al., 2020). The late time VHE emission in this object also implies that the
blastwave has slowed down considerably.

Figure 1 represents the multi-wavelength light curves for this source from
optical (r-band), X-ray, to VHE. We focused on the explanation for the late-
time emission of each band at $10^{4}<t<10^{6.7}$ s, corresponding to the
late-time afterglow range of the H.E.S.S. observations. In this time range a
set of values are chosen such that the requirement in Equation (13) and
$Y\approx 1$ are satisfied, as well as to fit each wave band simultaneously,
with the exception of the photon energy. These values are: $E_{k}=4.7\times
10^{54}$ erg, $\Gamma$=318, $t_{\rm dec,w}=100$ s, $A_{\rm\star}=1\times
10^{-1}$, $p$=2.18, $\epsilon_{e}=3\times 10^{-2}$, $\epsilon_{B}=3\times
10^{-2}$ (with the latter two values chosen such that $Y\approx 1$ is fixed),
and $\phi$=10 (see Table 1). These parameters are used to determine the
quantities as described in Section 2.1 and 2.2 to calculate the luminosity and
therefore the EC energy flux. The optical data are described by SR at energies
$0.5$ eV, and the _Swift_ -XRT data by SR at energies $2.4$ keV (data from
figure (1) in Abdalla et al. (2019)). The set of parameters result in an
estimation for $t_{\rm dec,w}$ that excludes the early time data points. We
have chosen the parameter space as such that we can compare the EC emission to
the H.E.S.S. data at energies $100$ GeV. We included the SSC emission as well
since this has been used primarily to explain the VHE observations, however,
for our chosen set of model parameters the SSC energy flux is expected to be
significantly lower (Joshi & Razzaque, 2021). We applied Equation (17) to both
the SSC and EC since both these emission mechanisms can fit the data at VHE
energies, considering only the EBL models of Finke et al. (2010) and Inoue et
al. (2013) with the former model estimating the attenuation slightly larger
than the latter one. In the ISM case (not shown), the fit to the optical and
X-ray data requires a steeper slope, and the slope of the EC fit is too steep.
Therefore we could not find a set of parameters to fit the observations in the
different wavebands simultaneously, including the upper limit of the H.E.S.S.
data, compared to the wind case. We gave preference to X-ray over the optical
data.

Figure 2 shows the energy spectrum associated with Figure 1, using the same
physical parameters at a $t=3.6\times 10^{4}$ s ($T_{0}+10$ hrs). This set of
parameters result in a good prediction of the EC fit to the VHE data, as well
as the synchrotron to the X-ray and optical (Fraija et al., 2019). We included
a joint fit of SR with both SSC and EC respectively, although the contribution
of SSC to the joint fit will be $10\%$ which is negligible compared to that of
the EC for the chosen values to explain the VHE observations via EC. In the
ISM case the SED, for these same parameter values as for the light curves in
the same case, does not produce a reasonable fit through the X-ray data, with
the SSC flux relatively higher than in the wind case. However, the EC emission
still result in a good fit to the H.E.S.S. observations. In both cases the EBL
model of Finke et al. (2010) results in a slightly larger flux, however the
model of Inoue et al. (2013) results in a slightly extended energy tail.

It is important to note that the energy flux of the H.E.S.S. observations in
Figure 1 is slightly higher than those in Figure 2. This is due to the
different methods used to do the light curve and spectral model fitting
compatible with the instrument. This is also the case for the observations of
GRB 190829A discussed in the next section.

Figure 3: Same as Figure 1 for long GRB 190829A, without optical observations.
We include observations for more than one night. These include the temporal
evolution of the energy flux detected in X-rays by _Swift_ -XRT (grey
squares), including those that were simultaneous with the _H.E.S.S._
observations (200 GeV-4 TeV; Abdalla et al. 2021), with an additional night 3
for the latter. In this case the combined fit of EC (blue) with SSC (orange)
is indicated by the dark red lines.

Figure 4: Same as Figure 2, we show the SED for night 1 (top) and night 2
(bottom), with the signal for night 3 very low, for long GRB 190829A without
optical observations. The shaded grey regions show the spectrum and
uncertainty of the _Swift_ -XRT data (Abdalla et al., 2021), and the green
lines indicate the joint spectral fits of the synchrotron (grey) with the
combined SSC (orange) and EC (blue) components, indicated by the dark red
lines.

###  3.2 GRB 190829A

H.E.S.S. also detected VHE radiation from GRB 190829A, with a $z=0.0785$ and a
T${}_{90}=63$ s, above 3 TeV for the first night at $T_{0}+4.3$ hr, above 1
TeV for the second night at $T_{0}+27.2$ hr (Abdalla et al., 2021), and for
the third night no significant signal was detected. In the afterglow emission,
the X-ray and GeV-TeV (0.18 to 3.3 TeV) light curves follow similar decay
features (Abdalla et al., 2021).

In Figure 3 we present the multiwavelength light curves for this source in
X-ray and VHE. We consider the late-time emission of each band corresponding
to the H.E.S.S. observations for the afterglow in nights 1-3. In this time
range a set of values are chosen such that the requirement in Equation (13)
and $Y\approx 1$ are satisfied as well as to fit each wave band
simultaneously, with the exception of the photon energy. These values are:
$E_{k}=2\times 10^{52}$ erg, $\Gamma=42$, $t_{\rm dec,w}=1\times 10^{3}$ s,
$A_{\rm\star}=9\times 10^{-2}$, $p=2.073$, $\epsilon_{e}=9.1\times 10^{-2}$,
$\epsilon_{B}=9.1\times 10^{-2}$ (with the latter two values chosen such that
$Y\approx 1$ is fixed), and $\phi=10$ (see Table 1). The _Swift_ -XRT data are
taken from Abdalla et al. (2021) and include observations from night 1 and
night 2, and are described by SR at energies $3$ keV. There were no coinciding
observations for night 3 in the X-ray band. For this source a combined fit of
the SSC and EC components explains the H.E.S.S. data (including all three
nights) at energies $200$ GeV and found a reasonable fit to explain the VHE
observations. This is due to the strict X-ray data constraint. However, for
the chosen parameter space the SSC energy flux dominates the EC (Joshi &
Razzaque, 2021). We applied Equation (17) to both the SSC and EC since we fit
data at VHE energies, considering only the EBL models of Finke et al. (2010)
and Inoue et al. (2013) with both models estimating the attenuation at the
same flux level. In the ISM case (not shown), the fit to the X-ray data
requires a steeper slope, thereby missing the second and third night.
Therefore we could not find a set of parameters to fit the observations in the
different wavebands simultaneously, compared to the wind case.

Figure 4 shows the energy spectrum for night 1 (top panel) and night 2 (bottom
panel) associated with Figure 1, using the same physical parameters but with
$E_{k}=5.6\times 10^{51}$ erg, $p=2.05$, and $t_{n1}=1.55\times 10^{4}$ s
($T_{0}+4.3$ hrs) for night 1, and $E_{k}=5\times 10^{51}$ erg, $p=2.1$, and
$t_{n2}=9.79\times 10^{4}$ s ($T_{0}+27.2$ hrs) for night 2 (see Table 1
column Night 1 and Night 2, respectively). We varied the $E_{k}$ and $p$ to
fit the spectral data for each night, however the $E_{k}$ value chosen to fit
the light curve is an average of the values chosen for the spectra of both
nights. This resulted in a good prediction of a combination of the SSC and EC
components to the VHE data, as well as the synchrotron to the X-ray data. We
included a fit of SR with the combined SSC and EC components, although the
contribution of the EC flux to the combined fit is negligible compared to the
contribution of the SSC. In the ISM case (not shown) the SED, for these same
parameter values as for the light curves in the same case, does not produce a
reasonable fit through the VHE data, with the SSC spectrum not as hard as in
the wind case to fit the spectral observations. However, the SSC and EC
combination still results in a good fit to the H.E.S.S. observations. In both
cases of the EBL model of Inoue et al. (2013) result in a slightly larger
flux, however the model of Finke et al. (2010) result in a slightly extended
energy tail.

##  4 DISCUSSIONS AND CONCLUSIONS

Table 1:  Summary of the afterglow parameter choices (with $Y$ a fixed value) from fitting the multiwavelength light curves and spectra simultaneously.  Parameter | GRB 180720B | GRB 190829A
---|---|---
|  | Light curve | Night 1 | Night 2
$E_{k}$ (erg) | $4.7\times 10^{54}$ | $2\times 10^{52}$ | $5.6\times 10^{51}$ | $5\times 10^{51}$
$\Gamma$ | 318 | 42 | 42 | 42
$t_{\rm dec,w}$ (s) | $100$ | $1\times 10^{3}$ | $1\times 10^{3}$ | $1\times 10^{3}$
$t$ (s) | $3.6\times 10^{4}$ | – | $1.55\times 10^{4}$ | $9.79\times 10^{4}$
$A_{\rm\star}$ | $1\times 10^{-1}$ | $9\times 10^{-2}$ | $9\times 10^{-2}$ | $9\times 10^{-2}$
$p$ | 2.18 | 2.073 | 2.05 | 2.1
$\epsilon_{e}$ | $3\times 10^{-2}$ | $9.1\times 10^{-2}$ | $9.1\times 10^{-2}$ | $9.1\times 10^{-2}$
$\epsilon_{B}$ | $3\times 10^{-2}$ | $9.1\times 10^{-2}$ | $9.1\times 10^{-2}$ | $9.1\times 10^{-2}$

In this work, we have shown that late-time VHE afterglow of GRB 180720B and
GRB 190829A detected by H.E.S.S. can be produced by the interaction of
relativistic electrons in the GRB blastwave with the CMB radiation. We have
explained the light-curves in optical, X-ray, VHE and corresponding SED of
these bursts using a synchrotron, SSC and EC radiation model. For GRB 180720B
the contribution due to the SSC mechanism is negligible based on the model
parameters used, but this is not the case for GRB 190829A which result in a
combiantion fit of SSC and EC. In particular, we found that one needs a
relatively low environment density, although not unusual, and
$\epsilon_{e}\sim\epsilon_{B}$ for the EC emission to be dominant. However,
even with $Y\approx 1$ imposed the EC is not always dominant as shown for GRB
190829A where the EC flux is very low and SSC dominating the fit. A summary of
the model parameters are listed in Table 1. For both GRBs the wind scenario
was favoured over the ISM. Also, the EBL correction for the two different EBL
models do not deviate too much. For GRB 190829A we used parameters that fit
the light curves and SEDs simultaneously, but $E_{k}$ and $p$ of the blastwave
for both varies slightly due to the spectral shape for night 1 being harder
than night 2. In the light curve we tried to fit three nights at once, again
resulting in a slightly different $E_{k}$ and $p$. Thus, fitting GRB 190829A
is more complex since the observations are spread over three nights and with
limited observations available for GRB 180720B and GRB 190829A one can not get
a better fit of the model to the data. The fit of X-ray is preferred over
optical since the host galaxy contributes the vast majority of the optical
emission.

Our EC modelling of VHE emission differs from the SSC model used for fitting
VHE data from GRB 180720B where the blastwave traverses in a constant density
ISM (Wang et al., 2019; Fraija et al., 2019). The isotropic energy released in
the gamma rays from GRB 190829A during the prompt emission was $\sim 10^{50}$
erg (Chand et al., 2020) and termed this as an example of low-luminosity (LL)
GRB. We found that indeed the kinetic energy of the burst is lower than that
of GRB 180720B by more than an order of magnitude. SSC models were also used
to explain the VHE gamma-ray emission from GRB 190829A (Zhang et al., 2021a;
Salafia et al., 2022) and its low luminosity was also used as a proxy for off-
axis VHE emission (Sato et al., 2021). Further, EC models with target photons
from the early-afterglow flare was used to explain VHE radiation (Zhang et
al., 2021b). A photo-hadronic model was also used to produce VHE gamma rays in
this burst (Sahu et al., 2022).

In summary, we found that inverse-Compton scattering of CMB photons by
relativistic electrons in the GRB blastwave can explain VHE emission detected
for some GRBs at late time ($\sim$ hours) when the energy density in the CMB
radiation becomes comparable to that of the magnetic energy density in the
blastwave. We found that the energy densities in the magnetic field and in
relativistic electrons are both of the order of a few percent of the blastwave
kinetic energy.

## Acknowledgements

We thank the H.E.S.S. collaboration for sharing the data for the GRBs studied
in this paper and Maria Petropolou for discussion. We also thank the MNRAS
Referee for constructive suggestions that improved the paper. The research of
MB was supported by a GES fellowship of University of Johannesburg, where this
work was completed. SR acknowledges support from the National Research
Foundation, South Africa for a BRICS STI grant and other grants from the
University of Johannesburg Research Council. JCJ was partially supported by a
grant from the University of Johannesburg Research Council to visit CAPP.

## Data Availability

All data used in this article are publicly available from the H.E.S.S.
Collaboration (https://www.mpi-hd.mpg.de/HESS/) and from the _Swift_ -XRT
website: https://www.swift.ac.uk/analysis/xrt/.

## References

  * Abbasi et al. (2023) Abbasi R., et al., 2023, ApJ, 946, L26
  * Abdalla et al. (2019) Abdalla H., et al., 2019, Nature, 575, 464
  * Abdalla et al. (2021) Abdalla H., et al., 2021, Science, 372, 1081
  * Alves Batista (2022) Alves Batista R., 2022, arXiv e-prints, p. arXiv:2210.12855
  * Asano et al. (2009) Asano K., Inoue S., Mészáros P., 2009, ApJ, 699, 953
  * Baktash et al. (2022) Baktash A., Horns D., Meyer M., 2022, arXiv e-prints, p. arXiv:2210.07172
  * Barniol Duran & Kumar (2011) Barniol Duran R., Kumar P., 2011, MNRAS, 412, 522
  * Blandford & McKee (1976) Blandford R. D., McKee C. F., 1976, Physics of Fluids, 19, 1130
  * Böttcher & Dermer (1998) Böttcher M., Dermer C. D., 1998, ApJ, 499, L131
  * Chand et al. (2020) Chand V., et al., 2020, ApJ, 898, 42
  * Cheung (2022) Cheung K., 2022, arXiv e-prints, p. arXiv:2210.14178
  * Chiang & Dermer (1999) Chiang J., Dermer C. D., 1999, ApJ, 512, 699
  * Da Vela et al. (2023) Da Vela P., Martí-Devesa G., Saturni F. G., Veres P., Stamerra A., Longo F., 2023, Phys. Rev. D, 107, 063030
  * Das & Razzaque (2023) Das S., Razzaque S., 2023, A&A, 670, L12
  * Derishev & Piran (2021) Derishev E., Piran T., 2021, ApJ, 923, 135
  * Dichiara et al. (2022) Dichiara S., Gropp J. D., Kennea J. A., Kuin N. P. M., Lien A. Y., Marshall F. E., Tohuvavohu A., Williams M. A., 2022, GCN Circ., 32632
  * Fan et al. (2008) Fan Y.-Z., Piran T., Narayan R., Wei D.-M., 2008, MNRAS, 384, 1483
  * Finke & Razzaque (2023) Finke J. D., Razzaque S., 2023, ApJ, 942, L21
  * Finke et al. (2010) Finke J. D., Razzaque S., Dermer C. D., 2010, ApJ, 712, 238
  * Fraija et al. (2019) Fraija N., et al., 2019, ApJ, 885, 29
  * Fraija et al. (2021) Fraija N., Veres P., Beniamini P., Galvan-Gamez A., Metzger B. D., Barniol Duran R., Becerra R. L., 2021, ApJ, 918, 12
  * Galanti et al. (2022) Galanti G., Roncadelli M., Tavecchio F., 2022, arXiv e-prints, p. arXiv:2211.06935
  * Gehrels & Razzaque (2013) Gehrels N., Razzaque S., 2013, Frontiers of Physics, 8, 661
  * Georganopoulos et al. (2001) Georganopoulos M., Kirk J. G., Mastichiadis A., 2001, ApJ, 561, 111
  * Ghisellini et al. (2020) Ghisellini G., et al., 2020, A&A, 636, A82
  * Gill & Granot (2022) Gill R., Granot J., 2022, Galaxies, 10, 74
  * Huang et al. (2022a) Huang Z.-Q., Kirk J. G., Giacinti G., Reville B., 2022a, ApJ, 925, 182
  * Huang et al. (2022b) Huang Y., Hu S., Chen S., Zha M., Liu C., Yao Z., Cao Z., 2022b, GCN Circ., 32677
  * Inoue et al. (2013) Inoue Y., Inoue S., Kobayashi M. A. R., Makiya R., Niino Y., Totani T., 2013, ApJ, 768, 197
  * Isravel et al. (2022) Isravel H., Pe’er A., Begue D., 2022, arXiv e-prints, p. arXiv:2210.02363
  * Joshi & Razzaque (2021) Joshi J. C., Razzaque S., 2021, MNRAS, 505, 1718
  * Klinger et al. (2023) Klinger M., Tak D., Taylor A. M., Zhu S. J., 2023, MNRAS, 520, 839
  * Kumar & Zhang (2015) Kumar P., Zhang B., 2015, Phys. Rep., 561, 1
  * LHAASO Collaboration et al. (2023) LHAASO Collaboration et al., 2023, Science, 380, 1390
  * Laskar et al. (2023) Laskar T., et al., 2023, ApJ, 946, L23
  * Lesage et al. (2023) Lesage S., et al., 2023, arXiv e-prints, p. arXiv:2303.14172
  * Li & Ma (2023) Li H., Ma B.-Q., 2023, Astroparticle Physics, 148, 102831
  * Li & Wang (2023) Li L., Wang Y., 2023, arXiv e-prints, p. arXiv:2302.06116
  * Liu et al. (2023) Liu R.-Y., Zhang H.-M., Wang X.-Y., 2023, ApJ, 943, L2
  * MAGIC Collaboration et al. (2019a) MAGIC Collaboration et al., 2019a, Nature, 575, 455
  * MAGIC Collaboration et al. (2019b) MAGIC Collaboration et al., 2019b, Nature, 575, 459
  * Mao & Wang (2023) Mao J., Wang J., 2023, ApJ, 946, 89
  * Mészáros (2006) Mészáros P., 2006, Reports on Progress in Physics, 69, 2259
  * Nakagawa et al. (2023) Nakagawa S., Takahashi F., Yamada M., Yin W., 2023, Physics Letters B, 839, 137824
  * Nakar et al. (2009) Nakar E., Ando S., Sari R., 2009, ApJ, 703, 675
  * Piran (2004) Piran T., 2004, Reviews of Modern Physics, 76, 1143
  * Razzaque (2013) Razzaque S., 2013, Phys. Rev. D, 88, 103003
  * Razzaque et al. (2010) Razzaque S., Dermer C. D., Finke J. D., 2010, The Open Astronomy Journal, 3, 150
  * Ren et al. (2023) Ren J., Wang Y., Zhang L.-L., Dai Z.-G., 2023, ApJ, 947, 53
  * Ronchi et al. (2020) Ronchi M., et al., 2020, A&A, 636, A55
  * Sahu & Fortín (2020) Sahu S., Fortín C. E. L., 2020, ApJ, 895, L41
  * Sahu et al. (2022) Sahu S., Valadez Polanco I. A., Rajpoot S., 2022, ApJ, 929, 70
  * Salafia et al. (2022) Salafia O. S., et al., 2022, ApJ, 931, L19
  * Sari & Esin (2001) Sari R., Esin A. A., 2001, ApJ, 548, 787
  * Sari et al. (1998) Sari R., Piran T., Narayan R., 1998, ApJ, 497, L17
  * Sato et al. (2021) Sato Y., Obayashi K., Yamazaki R., Murase K., Ohira Y., 2021, MNRAS, 504, 5647
  * Sato et al. (2023) Sato Y., Murase K., Ohira Y., Yamazaki R., 2023, MNRAS, 522, L56
  * Smirnov & Trautner (2022) Smirnov A. Y., Trautner A., 2022, arXiv e-prints, p. arXiv:2211.00634
  * Veres et al. (2022) Veres P., Burns E., Bissaldi E., Lesage S., Roberts O., 2022, GCN Circ., 32636
  * Wang et al. (2015) Wang X.-G., et al., 2015, ApJS, 219, 9
  * Wang et al. (2019) Wang X.-Y., Liu R.-Y., Zhang H.-M., Xi S.-Q., Zhang B., 2019, ApJ, 884, 117
  * Wang et al. (2020) Wang Z.-R., Xi S.-Q., Liu R.-Y., Xue R., Wang X.-Y., 2020, Phys. Rev. D, 101, 083004
  * Williams et al. (2023) Williams M. A., et al., 2023, ApJ, 946, L24
  * Yamasaki & Piran (2022) Yamasaki S., Piran T., 2022, MNRAS, 512, 2142
  * Zhang & Mészáros (2001) Zhang B., Mészáros P., 2001, ApJ, 559, 110
  * Zhang et al. (2021a) Zhang L.-L., Ren J., Huang X.-L., Liang Y.-F., Lin D.-B., Liang E.-W., 2021a, ApJ, 917, 95
  * Zhang et al. (2021b) Zhang B. T., Murase K., Veres P., Mészáros P., 2021b, ApJ, 920, 55
  * Zhang et al. (2023) Zhang B. T., Murase K., Ioka K., Song D., Yuan C., Mészáros P., 2023, ApJ, 947, L14


# Interaction induced fractional Bloch and tunneling oscillations

Ramaz Khomeriki1,2 Dmitry O. Krimer1, Masudul Haque1 and Sergej Flach1 ${\
}^{1}$Max-Planck Institute for the Physics of Complex Systems, Nöthnitzer Str.
38, 01187 Dresden, Germany
${\ }^{2}$Physics Department, Tbilisi State University, Chavchavadze 3, 0128
Tbilisi, Georgia

###### Abstract

We study the dynamics of few interacting bosons in a one-dimensional lattice
with dc bias. In the absence of interactions the system displays single
particle Bloch oscillations. For strong interaction the Bloch oscillation
regime reemerges with fractional Bloch periods which are inversely
proportional to the number of bosons clustered into a bound state. The
interaction strength is affecting the oscillation amplitude. Excellent
agreement is found between numerical data and a composite particle dynamics
approach. For specific values of the interaction strength a particle will
tunnel from the interacting cloud to a well defined distant lattice location.

###### pacs:

67.85.-d, 37.10.Jk, 03.65.Ge

Bloch oscillations BO  in dc biased lattices are due to wave interference and
have been observed in a number of quite different physical systems: atomic
oscillations in Bose-Einstein Condensates (BEC) BEC , light intensity
oscillations in waveguide arrays OL , and acoustic waves in layered and
elastic structures kosevich , among others.

Quantum many body interactions can alter the above outcome. A mean field
treatment will make the wave equations nonlinear and typically nonintegrable.
For instance, for many atoms in a Bose-Einstein condensate, a mean field
treatment leads to the Gross-Pitaevsky equation for nonlinear waves. The main
effect of nonlinearity is to deteriorate Bloch oscillations, as recently
studied experimentally salerno  and theoretically DJ98 ; krimer ; kgk10 .

In contrast, we will explore the fate of Bloch oscillations for quantum
interacting few-body systems. This is motivated by recent experimental advance
zoller  in monitoring and manipulating few bosons in optical lattices. Few
body quantum systems are expected to have finite eigenvalue spacings,
consequent quasiperiodic temporal evolution and phase coherence. In a recent
report on interacting electron dynamics spectral evidence for a Bloch
frequency doubling was reported doubling . On the other hand, it has been also
recently argued that Bloch oscillations will be effectively destroyed for few
interacting bosons kolovsky .

In the present paper we show that for strongly interacting bosons a coherent
Bloch oscillation regime reemerges. If the bosons are clustered into an
interacting cloud at time $t=0$, the period of Bloch oscillations will be a
fraction of the period of the noninteracting case, scaling as the inverse
number of interacting particles (Fig.1). The amplitude (spatial extent) of
these fractional Bloch oscillations will decrease with increasing interaction
strength. For specific values of the interaction, one of the particles will
leave the interacting cloud and tunnel to a possibly distant and well defined
site of the lattice. For few particles the dynamics is always quasiperiodic,
and a decoherence similar to the case of a mean field nonlinear equation
krimer  will not take place.

Figure 1: Time evolution of the probability density function (PDF) $P_{j}(t)$
for the interaction constant $U=3$ and dc field $E=0.05$ and different
particle numbers initially occupying a single site at $t=0$. (a) shows one
particle Bloch oscillations with the conventional Bloch period $2\pi/E$, while
(b), (c) and (d) display two, three and four particle oscillations with the
periods $2\pi/(2E)$, $2\pi/(3E)$ and $2\pi/(4E)$, respectively.

We consider the Bose-Hubbard model with a dc field:

${\cal\hat{H}}=\sum\limits_{j}\left[t_{1}\left(\hat{b}_{j+1}^{+}\hat{b}_{j}+\hat{b}_{j}^{+}\hat{b}_{j+1}\right)+Ej\hat{b}_{j}^{+}\hat{b}_{j}+\frac{U}{2}\hat{b}_{j}^{+}\hat{b}_{j}^{+}\hat{b}_{j}\hat{b}_{j}\right]$
(1)

where $\hat{b}_{j}^{+}$ and $\hat{b}_{j}$ are standard boson creation and
annihilation operators at lattice site $j$; the hopping $t_{1}=1$; $U$ and $E$
are the interaction and dc field strengths, respectively. To study the
dynamics of $n$ particles we use the orthonormal basis of states $|{\bf
k}\rangle\equiv|k_{1},k_{2},...,k_{n}\rangle=b_{k_{1}}^{+}b_{k_{2}}^{+}...b_{k_{n}}^{+}|0\rangle$
where $|0\rangle$ is the zero particle vacuum state, and $k_{1}\leq
k_{2}...\leq k_{n}$ are lattice site indices (for instance, in the case of two
particles the state representation is mapped to the triangle). The
eigenvectors $|\nu\rangle$ of Hamiltonian (1) with eigenvalues $\lambda_{\nu}$
are then given by:

$|\nu\rangle=\sum\limits_{{\bf k}}A_{{\bf k}}^{\nu}|{\bf
k}\rangle\;,\qquad{\cal\hat{H}}|\nu\rangle=\lambda_{\nu}|\nu\rangle$ (2)

where the eigenvectors $A_{{\bf k}}^{\nu}\equiv\langle{\bf k}|\nu\rangle$ and
the time evolution of a wave function $|\Psi(t)\rangle$ is given by

$|\Psi(t)\rangle=\sum\limits_{\nu}\Phi_{\nu}e^{-i\lambda_{\nu}t}|\nu\rangle,\qquad\Phi_{\nu}\equiv\langle\nu|\Psi(0)\rangle.$
(3)

We monitor the probability density function (PDF)
$P_{j}(t)=\langle\Psi(t)|\hat{b}_{j}^{+}\hat{b}_{j}|\Psi(t)\rangle/n$, which
can be also computed using the eigenvectors and eigenvalues:

$P_{j}(t)=\dfrac{1}{n}\sum\limits_{\nu,\mu}\Phi_{\nu}\Phi_{\mu}^{*}e^{i(\lambda_{\mu}-\lambda_{\nu})t}\langle\mu|\hat{b}_{j}^{+}\hat{b}_{j}|\nu\rangle\;.$
(4)

In Fig.1 we show the evolution of $P_{j}(t)$ for $U=3$, $E=0.05$ and
$n=1,2,3,4$ with initial state $k_{1}=k_{2}=...=k_{n}\equiv p$, i.e. when all
particles are launched on the same lattice site $p$. For $n=1$ we observe the
usual Bloch oscillations with period $T=2\pi/E$ (Fig.1(a) and below). Due to
the small value of $E$, the amplitudes of oscillations are large. However,
with increasing number of particles, we find that the oscillation period is
reduced according to $2\pi/(nE)$, and at the same time the amplitude of
oscillations is also reduced.

## One particle case:

For $n=1$ the interaction term in (1) does not contribute. The eigenvalues
$\lambda_{\nu}=E\nu$ (with $\nu$ being an integer) form an equidistant
spectrum which extends over the whole real axis - the Wannier-Stark ladder.
The corresponding eigenfunctions obey the generalized translational invariance
$A_{k+\mu}^{\nu+\mu}=A_{k}^{\nu}$ BO  and are given by the Bessel function
$J_{k}(x)$ of the first kind book ; SL

$\displaystyle A_{k}^{\nu}=J_{k}^{\nu}\equiv J_{k-\nu}(2/E).$ (5)

All eigenvectors are spatially localized with an asymptotic decay
$|A_{k\rightarrow\infty}^{0}|\rightarrow\left(1/E\right)^{k}/k!$, giving rise
to the well-known localized Bloch oscillations with period $T_{B}=2\pi/E$. The
localization volume $\mathcal{L}$ of a single particle eigenstate
characterizes its spatial extent. It follows $\mathcal{L}\propto-[E\cdot\ln
E]^{-1}$ for $E\rightarrow 0$ and $\mathcal{L}\rightarrow 1$ for
$E\rightarrow\infty$ krimer . For $E=0.05$ the single particle oscillates with
amplitude of the order of $2\mathcal{L}\approx 160$ (Fig.1(a)). According to
Eqs. (4) and (5) the probability density function is given by:

$P_{j}(t)=\sum\limits_{\nu,\mu}J_{p}^{\nu}J_{p}^{\mu}J_{j}^{\nu}J_{j}^{\mu}e^{iE({\mu}-\nu)t}.$
(6)

Figure 2: Eigenvalue spectrum for $n=2$, $E=0.5$ and different interaction
constants $U$. The eigenvalues are displayed only for eigenvectors localized
in the center of the lattice (we select the 32 eigenstates which overlap most
strongly with the center of the lattice) (a): $U=0$, the spectrum is
equidistant with spacing $E$ and degenerate. (b): $U=2$, the degeneracy is
lifted. (c): $U=15$, the spectrum decomposes into two subspectra, with two
different equidistant spacings - $E$ and $2E$. Graph (d) displays the
eigenvalue spectrum of the 32 central eigenfunctions as a function of $U$.

## Two particle case (n=2):

For $U=0$ the eigenfunctions of the Hamiltonian (1) are given by tensor
products of the single particle eigenstates:

$|\mu,\nu\rangle=\sqrt{\frac{2-\delta_{\mu,\nu}}{2}}\sum\limits_{k,j}J_{k}^{\mu}J_{j}^{\nu}\hat{b}_{k}^{+}\hat{b}_{j}^{+}|0\rangle\;,\;\mu\leq\nu\;.$
(7)

The corresponding eigenvalues form an equidistant spectrum which is highly
degenerate:

${\cal\hat{H}}|\mu,\nu\rangle=(\mu+\nu)E|\mu,\nu\rangle$ (8)

For the above initial condition $k_{1}=k_{2}\equiv p$ the expression for the
PDF (6) is still valid (actually it is for any number of noninteracting
particles), with the same period $2\pi/E$ of Bloch oscillations as in the
single particle case.

Figure 3:  Upper plots: PDF for $E=0.1$, $n=2$, single site initial occupancy
and different interaction constants. For $U=0$ we find single particle Bloch
oscillations. For $U=4$ fractional Bloch oscillations take place, in agreement
with (11). Lower plots: probability density of the evolved wave function
(darker regions correspond to larger probablilities) after one half of the
respective Bloch period. For $U=0$ the two particles are with equal
probability close to each other and at maximal separation. For $U=4$ the two
particles avoid separation and form a composite particle which coherently
oscillates in the lattice. In lower graphs we use triangle $k<m$ mapping for
indistinguishable two particle state representation (index $m$ increases from
the right to the left).

For nonvanishing interaction the degeneracy of the spectrum is lifted, and the
eigenvalues of overlapping states are not equidistant any more (Fig.2).
Therefore we observe quasiperiodic oscillations which however are still
localizing the particles. For even larger values of $U$ the basis states with
two particles on the same site will shift their energies by $U$ exceeding the
hopping $2t_{1}$. Therefore for $U>2t_{1}$ the spectrum will be decomposed
into two nonoverlapping parts - a noninteracting one which excludes double
occupancy and has equidistant spacing $E$, and an interacting part which is
characterized by almost complete double occupancy and has corresponding
equidistant spacing $2E$, which is the cost of moving two particles from a
given site to a neighboring site. Some initial state can overlap strongly with
eigenstates from one or the other part of the spectrum, and therefore result
in different Bloch periods. In particular, when launching both particles on
the same site, one strongly overlaps with the interacting part of the spectrum
and observes a fractional Bloch period $2\pi/(2E)$.

In order to calculate the amplitude of these fractional Bloch oscillations, we
note that for $E=0$ there exists a two-particle bound state band of extended
states with band width $\sqrt{U^{2}+16}-U$ eilbeck . For large $U$ the bound
states are again almost completely described by double occupancy. Therefore we
can construct an effective Hamiltonian for a composite particle of two bound
bosons:

${\cal\hat{H}}\approx\sum\limits_{j}\left[t_{2}\left(\hat{R}_{j+1}^{+}\hat{R}_{j}+\hat{R}_{j}^{+}\hat{R}_{j+1}\right)+2Ej\hat{R}_{j}^{+}\hat{R}_{j}\right]$
(9)

where $\hat{R}_{j}^{+}$ and $\hat{R}_{j}$ are creation and annihilation
operators at lattice site $j$ of the composite particle (two bosons on the
same site) with the effective hopping

$t_{2}=\frac{\sqrt{U^{2}+16}-U}{4}.$ (10)

The corresponding PDF is given by

$P_{j}(t)=\sum\limits_{\nu,\mu}A_{p}^{\nu}A_{p}^{\mu}A_{j}^{\nu}A_{j}^{\mu}e^{i2E({\mu}-\nu)t}\;.$
(11)

The composite particle eigenvectors $A_{p}^{\nu}=J_{\nu-p}(2t_{2}/(2E))$ are
again expressed through Bessel functions, but with a modified argument as
compared to the single particle case. Bloch oscillations will evolve with
fractional period $2\pi/(2E)$ as observed in Fig.1(b). The amplitude of the
oscillations is reduced with increasing $U$ since the hopping constant $t_{2}$
is reduced (Fig.3). For $U=3$ it follows $t_{2}=0.5$, and together with the
doubled Bloch frequency the localization volume should be reduced by a factor
of 4 as compared to the single particle case. This is precisely what we find
when comparing Fig.1(a,b): for $n$=1 the amplitude is 160 sites, while for
$n=2$ it is 40 sites. In the lower plots in Fig.3 we show the probability
density of the wave functions $|\langle\Psi(t)|{\bf k}\rangle|^{2}$ after one
half of the respective Bloch period in the space of the two particle
coordinates with $k_{1}=k$ and $k_{2}=m$. For $U=0$ both particles are with
high probability at a large distance from each other. Therefore the density is
large not only for $k=m$ (the two particles are at the same site), but also
for $k=5$, $m=85$ (the two particles are at maximum distance). However, for
$U=4$ we find that the two particles, which initially occupy the site $p=45$,
do not separate, and the density is large only along the diagonal $k=m$ with
$35\leq k\leq 55$. (For $U=4$, the localization volume is $\sim 20$.)
Therefore, the two particles indeed form a composite state and travel
together.

##  $n$ particle case:

We proceed similar to the case $n=2$ and estimate perturbatively the effective
hopping constant for a composite particle of $n$ bosons. For that we use the
calculated width of the $n$-particle bound state band for $E=0$ eilbeck . In
leading order of $1/U$ it reads eilbeck :

$t_{n}\simeq\frac{n}{U^{n-1}(n-1)!}.$ (12)

For $n=2$ the above expression gives $t_{2}\simeq 2/U$, the first expansion
term of the exact relation for two bosons (10). The corresponding composite
particle Hamiltonian

${\cal\hat{H}}\approx\sum\limits_{j}\left[t_{n}\left(\hat{R}_{j+1}^{+}\hat{R}_{j}+\hat{R}_{j}^{+}\hat{R}_{j+1}\right)+nEj\hat{R}_{j}^{+}\hat{R}_{j}\right]\;.$
(13)

The PDF is given by

$P_{j}(t)=\sum\limits_{\nu,\mu}A_{p}^{\nu}A_{p}^{\mu}A_{j}^{\nu}A_{j}^{\mu}e^{inE({\mu}-\nu)t}\;,$
(14)

and the composite particle eigenvectors $A_{p}^{\nu}=J_{\nu-p}(2t_{n}/(nE))$.
Bloch oscillations will evolve with fractional period $2\pi/(nE)$ as observed
in Fig.1(c,d). The amplitude of the oscillations is reduced with increasing
$U$ since the hopping constant $t_{n}$ is reduced. For $U=3$ and $n=3$ it
follows $t_{3}=0.17$, and for $n=4$ we have $t_{4}=0.01$. This leads to a
reduction factor 18 and 400 respectively as compared to the single particle
amplitude and yields amplitudes of the order of 9 and 0.5 respectively, which
is in good agreement with the numerically observed amplitudes (10 and 2 sites
respectively) in Fig.1(c,d).

Figure 4: Time averaged and normalized localization volume $L$ of the
wavepacket which emerge from two initial distributions as a function of $U$
for $E=5$. Red (grey) curve: two particles are launched on the same site.
Orange (light grey) curve: two particles are launched on adjacent sites.
Inset: PDF for $U=19.79$, with clearly observed tunneling oscillations.

## Tunneling oscillations:

For $n=1$ the amplitude of Bloch oscillations is less than one site if $E\geq
10$ krimer . Thus, for $n\geq 2$ and increasing values of $U$, the amplitude
of fractional Bloch oscillations will be less than one site if
$EU^{n-1}(n-1)!\geq 10$. Then, $n$ particles launched on the same lattice site
$p$ will be localized on that site for all times. The energy of that state
will be $n((n-1)U/2+pE)$. If however one particle will be moved to a different
location with site $q$, then the energy would change to
$(n-1)((n-2)U/2+pE)+qE$. For specific values of $U$ these two energies will be
equal:

$(n-1)U=dE\;,\;d=q-p\;.$ (15)

In such a case, one particle will leave the interacting cloud at site $p$ and
tunnel to site $q$ at distance $d$ from the cloud, then tunnel back and so on,
following effective Rabi oscillation scenario between the states $p,p\rangle$
and $|p,q\rangle$. This process will appear as an asymmetric oscillation of a
fraction of the cloud either up or down the field gradient (depending on the
sign of $U$). We calculate the tunneling splitting of these two states using
higher order perturbation theory, for an example see Ref.afko96 . The
tunneling time is then obtained as

$\tau_{tun}\simeq\frac{\pi}{\sqrt{n}}E^{d-1}(d-1)!\;.$ (16)

In order to observe these tunneling oscillations, we compute the time averaged
second moment
$\overline{m_{2}}=\overline{\sum\limits_{j}j^{2}P_{j}(t)-\left(\sum_{j}jP_{j}(t)\right)^{2}}$
of the PDF $P$. Then an effective time-averaged volume of the interacting
cloud is taken to be $L=\sqrt{12\overline{m_{2}}}+1$. We launch $n=2$
particles at site $p=40$ and plot the ratio $L(U)/L(U=0)$ in Fig.4 (blue solid
line). We find pronounced peaks at $U=E,2E,3E,4E$ which become sharper and
higher with increasing value of $U$. As a comparison we also compute the same
ratio for the initial condition when both particles occupy neighbouring sites
(dashed red line), for which the resonant structures are absent. According to
the above, the resonant structures correspond to a tunneling of one of the
particles to a site at distance $d=1,2,3,4$. The width of the peaks is
inversely proportional to the tunneling time $\tau_{tun}$, and the height
increases linearly with the tunneling distance $d$. In the inset in Fig.4, we
plot the time evolution of the PDF $P_{j}$ for $U=19.79$. We observe a clear
tunneling process from site $p=40$ to site $q=44$. The numerically observed
tunneling time is approximately 1730 time units, while our above prediction
(16) yields $\tau_{tun}\approx 1666$, in very good agreement with the
observations.

## Conclusions.

The above findings can be useful for control of the dynamics of interacting
particles. They can be also used as a testbed of whether experimental studies
deal with quantum many body states. One such testbed is the observation of
fractional Bloch oscillations, another one is the resonant tunneling of a
particle from an interacting cloud. An intriguing question is the way these
quantum coherent phenomena will disappear in the limit of many particles,
where classical nonlinear and nonintegrable wave mechanics are expected to
take over.

## Acknowledgements.

R. Kh. acknowledges financial support of the Georgian National Science
Foundation (Grant No GNSF/STO7/4-197).

## References

  * (1) F. Bloch, Z. Phys. 52, 555 (1928); C. Zener, Proc. R. Soc. London, Ser. A 145, 523 (1934); G.H. Wannier, Rev. Mod. Phys. 34, 645 (1962).
  * (2) B.P. Anderson, M.A. Kasevich, Science 282, 1686 (1998); M. Ben Dahan et al, Phys. Rev. Lett. 76, 4508 (1996); O. Morsch et al, Phys. Rev. Lett. 87, 140402 (2001); G. Ferrari et al, Phys. Rev. Lett. 97, 060402 (2006).
  * (3) T. Pertsch et al, Phys. Rev. Lett. 83, 4752 (1999); H. Trompeter et al, Phys. Rev. Lett. 96, 023901 (2006); F. Dreisow et al, Phys. Rev. Lett. 102, 076802 (2009).
  * (4) H. Sanchis-Alepuz et al,, Phys. Rev. Lett. 98, 134301 (2007); L. Gutierrez et al, Phys. Rev. Lett. 97, 114301 (2006);
  * (5) M. Gustavsson et al, Phys. Rev. Lett. 100, 080404 (2008).
  * (6) P.K. Datta and A.M. Jayannavar, Phys. Rev. B 58, 8170 (1998).
  * (7) D.O. Krimer, R. Khomeriki, S. Flach, Phys. Rev. E, 80, 036201 (2009).
  * (8) A. R. Kolovsky, E. A. Gomez EA and H. J. Korsch, Phys. Rev. A 81, 025603 (2010)
  * (9) K. Winkler, et al, Nature, 441, 853 (2006); W.C. Bakr, et al, Nature, 462, 74 (2006);
  * (10) W. S. Dias, et al, Phys. Rev. B, 76, 155124 (2007).
  * (11) A. Buchleitner, A.R. Kolovsky, Phys. Rev. Lett. 91, 253002 (2003).
  * (12) M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions (1972).
  * (13) H. Fukuyama, R.A. Bari, H.C. Fogedby, Phys. Rev. B 8 5579 (1973); P. Feuer, Phys. Rev., 88, 92 (1952).
  * (14) A.C. Scott, J.C. Eilbeck, H. Gilhoj, Physica D, 78, 194 (1994)
  * (15) S. Aubry et al, Phys. Rev. Lett. 76, 1607 (1996).