# Layer Coherent Phase in Double Layer graphene at $\nu_{1}=\nu_{2}=0$

Amartya Saha<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Kentucky, Lexington, KY 40506, USA  Ankur Das
<EMAIL_ADDRESS>Department of Condensed Matter Physics, Weizmann
Institute of Science, Rehovot, 76100 Israel

###### Abstract

In the recent advancement in graphene heterostructures, it is possible to
create a double layer tunnel decoupled graphene system that has a strong
interlayer electronic interaction. In this work, we restrict the parameters in
the low energy effective Hamiltonian using simple symmetry arguments. Then, we
study the ground state of this system in the Hartree-Fock approximation at
$\nu_{1}=\nu_{2}=0$. In addition to the phases found in monolayer graphene, we
found an existence of layer coherent phase which breaks the layer $U(1)$
symmetry. At non-zero Zeeman coupling strength ($E_{z}$), this layer coherent
state has a small magnetization, that vanishes when $E_{z}$ tends to zero. We
discuss the bulk gapless modes using the Goldstone theorem. We also comment on
the edge structure for the layer coherent phase.

##  I Introduction

After the discovery of the quantum Hall effect in two dimensional electron gas
(2DEG) Klitzing et al. (1980), the interest started to build towards bilayer
2DEG’s Chakraborty and Pietiläinen (1987) (and for general systems with
internal quantum number Halperin (1983)). In the quantum Hall regime the
ground state of double layer 2DEG has been explored theoretically Jain (2007);
Scarola and Jain (2001); Chakraborty and Pietiläinen (1987) and later
experimentally Eisenstein (2014) (and the references therein). It has been
observed that in 2DEG (in this case GaAs), when the interlayer distance is
small, the system at total filling fraction
$\nu_{T}=\nu_{1}+\nu_{2}=\frac{1}{2}$ ($\nu_{1}$, $\nu_{2}$ being the filling
fraction of each layer) and $1$ forms an incompressible QHSKellogg et al.
(2004); Gramila et al. (1991); Eisenstein and MacDonald (2004). The individual
layers at $\nu_{T}=1$ and $\frac{1}{2}$ has even denominator filling fractions
of $\frac{1}{2}$ and $\frac{1}{4}$ respectively which are known to be
compressible states. This phase is a layer coherent phase in which the
electron of one layer forms a bound state with the holes of the other layer
forming excitons and sponteneously breaks the layer $U(1)$ symmetry. We can
think of layer coherent states as either easy plane layer pseudospin
ferromagnet or electron-hole bound exciton Eisenstein (2014). This arises from
the conservation of particle number in the individual layers. Developments in
graphene technology have boosted the interest in quantum Hall effect in
graphene Novoselov et al. (2005); Zhang et al. (2005). Some recent experiments
in the quantum Hall regime in double layer graphene systems has been found to
exhibit layer coherent statesLi et al. (2016, 2019a); Liu et al. (2019).
Experimentally it is possible to fabricate double layer graphene with very
small interlayer distance $d$ ($\sim$ 2nm) where $d/l<1$ ($l$ being the
magnetic length) Gorbachev et al. (2012); Titov et al. (2013); Li et al.
(2016); Liu et al. (2017); Li et al. (2017); Britnell et al. (2012); Lee et
al. (2016); Liu et al. (2017, 2019); Li et al. (2019a) which was earlier
difficult to achieve in GaAs systems. The separator between the graphene
layers is made out of stacked hBN layers. Thus by changing the number of
stacked hBN layers the interlayer interaction can be tuned from weak to
strong. This induced a huge interest in the understanding and testing of the
double layers of graphene, Bernal-stacked bilayer graphene Li et al. (2016,
2017), twisted magic angle bilayer graphene Liu et al. (2020a) etc. There has
been some theoretical Tse et al. (2007); Katsnelson (2011); Peres et al.
(2011); Narozhny et al. (2012); Kharitonov and Efetov (2010) and experimental
Kim et al. (2011) studies to understand the Coulomb drag in double layer
graphene in zero magnetic field. Ref. Min et al., 2008 predicted that at
higher temperature at the zero magnetic field there can be a normal to
superfluid transition in double layer graphene.

In the presence of a ultra short-range (compared to the magnetic length $l$)
interaction, the Hamiltonian projected to the $n=0$ Landau level manifold for
monolayer graphene (MLG) has
$SU(2)_{\text{spin}}\times\left(U(1)\times\mathbb{Z}_{2}\right)_{\text{valley}}$
symmetry in the absence of a Zeeman field Alicea and Fisher (2006). There
exists four different possible phases namely ferromagnet (F), charge density
wave (CDW), Kekulé distorted phase (KD) and anti-ferromagnet (AF), which
becomes canted anti-ferromagnet (CAF) in presence of Zeeman coupling
Kharitonov (2012a). The predicted phase transition from CAF to F Herbut
(2007a, b); Kharitonov (2012b); Murthy et al. (2014, 2016); Knothe and
Jolicoeur (2015) has been verified in the experiment Young et al. (2014). This
understanding of symmetry has been used to study the ground state at
fractional fillings as well Sodemann and MacDonald (2014).

In the case of double monolayer graphene electron fillings of each layer
($\nu_{1}$ and $\nu_{2}$) can be controlled independently. For this system
many interesting quantum states have been found. Some of those states can be
explained using interlayer Jain composite fermion states Jain (2007), proposed
for double layer two-dimensional electron gas Li et al. (2019a). In this
manuscript we propose the relevant symmetry in double monolayer graphene which
restricts the interacting Hamiltonian to three parameters. Within the scope of
this letter, we restrict ourselves to understand the mean field ground state
when two layers of graphene are at $\nu_{1}=\nu_{2}=0$. We show that for
certain values of the parameters, the system goes into a layer coherent phase,
which has a small magnetization in presence of a Zeeman field. Increasing the
Zeeman field strength one can drive a second order phase transition from
magnetized layer coherent phase to the ferromagnetic phase.

Here we would like to emphasize that we want to find a low energy Hamiltonian
that is restricted by symmetry. We also focus on the translation invariant
ground state solutions of this low energy Hamiltonian in the allowed parameter
space. Our method does not talk about the details of the microscopic model but
only the low energy effective model.

We describe our assumptions, method and findings in a few sections. In section
II we describe the assumptions and our Hamiltonian. After that we describe the
results and the Goldstone modes in section III and IV respectively. We also
have some discussion over possible lattice models, experimental signatures in
section V. Then in section VI we summarize our findings and describe possible
application of this work.

##  II Assumptions and Model

We restrict our calculations to the $n=0$ Landau level. When the interaction
strength is much smaller than the cyclotron energy gap the Landau level mixing
can be ignored. In the strong interaction strength regime the form of the
effective theory gets dictated by the symmetry (discussed below) when we
integrate out the higher Landau levels. When the layers are far enough from
each other, we should get two MLG with no interlayer interaction. The valley
$U(1)$ for each layer is conserved in order to conserve the translational
symmetry in each layer separately. Here we make an additional assumption that
the global spin $SU(2)$ symmetry can be enhanced to spin $SU(2)$ symmetry for
each layer separately. For this to be the symmetry of this theory we assume
that interlayer spin-spin interaction is zero (or negligible). In the absence
of the inter layer tunneling it is justified that the Heisenberg term
($\vec{S}\cdot\vec{S}$) will be absent. Other than the Heisenberg a long range
spin dipole-dipole interaction between layers can break the spin $SU(2)$
symmetry in each layer to a global $SU(2)$ symmetry. However, the spin dipole-
dipole interaction falls as $r^{-4}$. As the distance between the layers is a
few nanometers, we choose to ignore this interaction. Thus under these
assumptions the only term that is allowed is the $\vec{S}\cdot\vec{S}$ on each
layer which has the spin $SU(2)$ symmetry in each layer.

From this understanding and keeping in mind that the number of particles in
each layer is fixed we propose our symmetry of the continuum model to be (in
the absence of Zeeman coupling) $\left[SU(2)_{\text{spin}}\otimes
U(1)_{\text{valley}}\right]$ for each layer, a global
$\left(\mathbb{Z}_{2}\right)_{\text{valley}}$ and
$\left(U(1)\otimes\mathbb{Z}_{2}\right)_{\text{layer}}$ for the layers.
symetry for each layer.This restricts the interacting part of the Hamiltonian
to only three parameters. We can write the Hamiltonian as,

$\displaystyle H=H_{0}+H_{\text{int}}.$ (1)

where $H_{0}$ is the one body term coming from Zeeman coupling such that,

$\displaystyle
H_{0}=-E_{z}\left(\sigma^{z}\otimes\tau^{0}\otimes\gamma^{0}\right).$ (2)

$H_{\text{int}}$, the 2 body interaction term which obeys the above mentioned
symmetry, is given by

$\displaystyle H_{int}=\frac{\pi l^{2}}{A}\sum_{\begin{subarray}{c}\vec{q}\\\
k_{1},k_{2}\end{subarray}}e^{-iq_{x}(k_{1}-k_{2}-q_{y})l^{2}-\frac{q^{2}l^{2}}{2}}$
$\displaystyle\Bigg{[}K_{xy}\sum_{i=1,2}\mathrel{\mathop{\mathchar
58\relax}}\vec{c}^{\dagger}_{k_{1}-q_{y}}\left(\sigma^{0}\otimes\tau^{i}\otimes
P_{L}\right)\vec{c}^{\phantom{\dagger}}_{k_{1}}\vec{c}^{\dagger}_{k_{2}+q_{y}}\left(\sigma^{0}\otimes\tau^{i}\otimes
P_{L}\right)\vec{c}^{\phantom{\dagger}}_{k_{2}}\mathrel{\mathop{\mathchar
58\relax}}$ $\displaystyle+K_{z}\mathrel{\mathop{\mathchar
58\relax}}\vec{c}^{\dagger}_{k_{1}-q_{y}}\left(\sigma^{0}\otimes\tau^{3}\otimes
P_{L}\right)\vec{c}^{\phantom{\dagger}}_{k_{1}}\vec{c}^{\dagger}_{k_{2}+q_{y}}\left(\sigma^{0}\otimes\tau^{3}\otimes
P_{L}\right)\vec{c}^{\phantom{\dagger}}_{k_{2}}\mathrel{\mathop{\mathchar
58\relax}}$ (3) $\displaystyle+L_{z}\mathrel{\mathop{\mathchar
58\relax}}\vec{c}^{\dagger}_{k_{1}-q_{y}}\left(\sigma^{0}\otimes\tau^{0}\otimes\gamma^{3}\right)\vec{c}^{\phantom{\dagger}}_{k_{1}}\vec{c}^{\dagger}_{k_{2}+q_{y}}\left(\sigma^{0}\otimes\tau^{0}\otimes\gamma^{3}\right)\vec{c}^{\phantom{\dagger}}_{k_{2}}\mathrel{\mathop{\mathchar
58\relax}}\Bigg{]}$

Here
$\vec{c}_{k}=\big{(}c^{\phantom{\dagger}}_{k,\uparrow,K,1},c^{\phantom{\dagger}}_{k,\downarrow,K,1},c^{\phantom{\dagger}}_{k,\uparrow,K^{\prime},1},c^{\phantom{\dagger}}_{k,\downarrow,K^{\prime},1},c^{\phantom{\dagger}}_{k,\uparrow,K,2},\allowbreak
c^{\phantom{\dagger}}_{k,\downarrow,K,2},\allowbreak
c^{\phantom{\dagger}}_{k,\uparrow,K^{\prime},2},c^{\phantom{\dagger}}_{k,\downarrow,K^{\prime},2}\big{)}^{T}$
presents the column vector of fermionic annihilation operators, $A$ is the
area of the sample, and $l$ is the magnetic length. The index $k_{i}$
represents the guiding centers in the Landau gauge. We use the convention
where $\sigma^{i},\tau^{i},\gamma^{i}$ represents the Pauli matrices in spin,
valley, and layer respectively. Here,
$P_{L}=\frac{\left(\gamma^{0}-\left(-1\right)^{L}\gamma^{3}\right)}{2}$ is the
layer projection operator to layer $L$. The parameters $K_{z}$ and $K_{xy}$
arises from the intralayer interactions and are same as the parameters $u_{z}$
and $u_{\perp}$ respectively as defined by Kharitonov in the monolayer
graphene caseKharitonov and Efetov (2010). The parameter $L_{z}$ is a function
of the distance between the layers ($d$) which should go to zero as $d$
becomes very large (disjoint MLG limit). Here we would like to comment that we
also added a capacitance term which is zero when both layers have equal
fillings Wen and Zee (1992),

$H_{cap}=\frac{g_{es}\pi
l^{2}}{A}\left(\rho_{1}(\mathbf{q}=0)-\rho_{2}(\mathbf{q}=0)\right)^{2}.$ (4)

where $\rho_{L}(\vec{q})$ is the Fourier transformed electron density operator
of $L$th layer and $g_{es}$ is the coupling strength of the capacitance term.

We define an order parameter $\Delta$ matrix which specifies the HF states
$|HF\rangle$,

$\langle
HF|c^{\dagger}_{k,s,\alpha,L}c^{\phantom{\dagger}}_{k,s,^{\prime}\alpha^{\prime},L^{\prime}}|HF\rangle=\delta_{\vec{k},\vec{k}^{\prime}}\Delta_{s^{\prime}\alpha^{\prime}L^{\prime},s\alpha
L}.$ (5)

where $s$ being the spin, $\alpha$ being the valley and $L$ being the layer
index. This $\Delta$ matrix can also be thought of as a sum of projection
operators of the four filled states at each momemtum. The $\Delta$ matrix
completely determines the single Slater determinant states and any order
parameters e.g. electron density, magnetization etc. can be calculated using
it. We assume that the HF states preserve translation symmetry i.e. the
guiding centers are a good quantum number. Hence, we drop the guiding center
label from $\Delta$ matrix. Since the capacitance term is a classical term, we
only keep the Hartree term and drop the Fock term. In the next section, we
discuss the $\Delta$ matrix of the different HF states.

Figure 1: Here we represented the Layer coherent (LC) phase. The layers here
are color coded (green lines and red lines). As shown here the states are
linear combination of different layer indices. The corresponding states are
also color coded as blue and light red.

##  III Results

At $\nu_{1}=\nu_{2}=0$ there are four occupied single particle states in the
spin-valley-layer space. For $L_{z}\geq 0$ we find that the phase diagram is
exactly the same as the phase diagram found for MLG in Ref. Kharitonov, 2012a.
The energies of the phases (defined as $E_{gs}=\langle
HF|H_{\text{int}}|HF\rangle$ for the proposed HF ground state)) depend on
$L_{z}$. For these phases, the layer $U(1)$ is not broken and the $\Delta$
matrix is block diagonal in the layer index. Each of this block is a four
dimensional matrix in the space of valley and spin. The four phases are,

  1. 1.

Charge Density Wave (CDW): CDW breaks the valley $\mathbb{Z}_{2}$ symmetry. At
the zero Landau level, different valley indices are pinned to the sublattices.
In this phase in each layer the alternate sites (A) in the lattice are
occupied and the other sites (B) are left unoccupied. The $\Delta$ matrix for
this phase is

$\Delta_{\text{CDW}}=\sigma^{0}\otimes(\tau^{0}+\tau^{3})\otimes\gamma^{0},$
(6)

and the energy is

$E_{\text{CDW}}=2\left(K_{z}-L_{Z}\right).$ (7)

  2. 2.

Kekulé Distorted (KD): This is a bond order phase where the valley $U(1)$
symmetry is broken. In lattice limit the spontaneous breaking of the valley
$U(1)$ symmetry leads to the translation symmetry breaking in each layer. This
phase doesn’t has any Goldstone modes. The $\Delta$ matrix for this phase will
be

$\Delta_{\text{KD}}=\frac{1}{2}\sigma^{0}\otimes(\tau^{0}+\tau^{1})\otimes\gamma^{0},$
(8)

with energy,

$E_{\text{KD}}=2\left(K_{xy}-L_{z}\right).$ (9)

  3. 3.

Ferromagnet (F): This phase breaks the spin $SU(2)$ symmetry in each layer.
Similarly, the $\Delta$ matrix and energy will be

$\displaystyle\Delta_{\text{F}}$
$\displaystyle=\frac{1}{2}(\sigma^{0}+\sigma^{3})\otimes\tau^{0}\otimes\gamma^{0}$
(10) $\displaystyle E_{\text{F}}$
$\displaystyle=-4E_{z}-2\left(2K_{xy}+K_{z}+L_{Z}\right).$ (11)

  4. 4.

Canted Anti-Ferromagnet (CAF): This phase breaks the spin $U(1)$ symmetry in
each layer. The $\Delta$ matrix is,

$\displaystyle\Delta_{\text{CAF}}$
$\displaystyle=\frac{1}{2}\Big{[}\sin\phi\left(\sigma^{1}\otimes\tau^{3}\otimes\gamma^{0}\right)$
$\displaystyle+\cos\phi\left(\sigma^{3}\otimes\tau^{0}\otimes\gamma^{0}\right)+\sigma^{0}\otimes\tau^{0}\otimes\gamma^{0}\Big{]}$
(12)

where $\phi$ is given by

$\cos\phi=\frac{E_{z}}{2|{K_{xy}}|},$ (13)

with energy,

$E_{\text{CAF}}=-4E_{z}-2\left(2K_{xy}+K_{z}+L_{Z}\right).$ (14)

At $E_{z}=0$ the states becomes a pure anti-ferromagnetic state. Increasing
the Zeeman field $E_{z}$ beyond $2|K_{xy}|$ can drive a continuous phase
transition from canted anti-ferromagnetic phase to ferromagnetic phase.

Figure 2: Here we present the Phase diagram of for $L_{z}=-0.5$ in the absence
of Zeeman coupling. The phase LC Phase appears near $K_{xy}=K_{z}=0$. All the
phase transitions here are First order.

Next we come to the spacial phase of the double layer graphene. For $L_{z}<0$,
we find there exists a layer coherent phase which breaks the Layer $U(1)$
symmetry (see 3). We find the layer coherent phase both in the presence and
absence of the Zeeman energy. For a non-zero Zeeman coupling, there are two
parameters (and operators which are connected by the left over ground state
symmetry) that,

$\displaystyle\Phi_{L}=\sigma^{1}\times\tau^{1}\times\gamma^{1}$ (15a)
$\displaystyle S^{z}=\frac{\sigma^{3}\times\tau^{0}\times\gamma^{0}}{2}.$
(15b)

Figure 3: Here we present the Phase diagram of for $L_{z}=-0.5$ and
$E_{z}=0.1$. As we can see that the MLC phase appears and there is a second
order transition from MLC to F as marked by the broken line.

The magnetic layer coherent phase occurs at $E_{z}\neq 0$ with $\langle
S^{z}\rangle\neq 0$. At $E_{z}=0$ we find $\langle S^{z}\rangle=0$, we call it
layer coherent phase (LC) (see fig. 1). We can write the $\Delta$ matrix for
the MLC phase as,

$\displaystyle\Delta_{\text{MLC}}=$
$\displaystyle\frac{1}{2}\Big{[}\sin\theta(\sigma^{1}\otimes\tau^{1}\otimes\gamma^{1})$
$\displaystyle+\cos\theta(\sigma^{3}\otimes\tau^{0}\otimes\gamma^{0})+\sigma^{0}\otimes\tau^{0}\otimes\gamma^{0}\Big{]}$
(16)

with $\cos\theta$ defined as,

$\cos\theta=\frac{2E_{z}}{\mathinner{\\!\left\lvert
2K_{xy}+K_{z}+2L_{z}\right\rvert}}.$ (17)

The energy of the phase is

$\displaystyle
E_{\text{MLC}}=-2K_{xy}-K_{z}-\frac{4E^{2}_{z}}{\mathinner{\\!\left\lvert
2K_{xy}+K_{z}+2L_{z}\right\rvert}}.$ (18)

Here $\langle\Phi_{L}\rangle=4\sin\theta$ and $\langle
S^{z}\rangle=2\cos\theta$. For $2E_{z}\geq\mathinner{\\!\left\lvert
2K_{xy}+K_{z}+2L_{z}\right\rvert}$ with $\theta=0$, this $\Delta$ will
represent a ferromagnetic ground state and for any other value of $\theta$ the
$\Delta$ matrix represents a magnetic layer coherent state (MLC). For the zero
Zeeman coupling, we have $\cos\theta=1\Rightarrow\langle S^{z}\rangle=0$, a
purely layer coherent state. The phase transition from MLC to F is a second
order transition (see Fig. 3). However, similar to the AF to F phase
transition, the LC to F phase transition is a first order transition. Thus all
phase transitions are first order at $E_{z}=0$ (see Fig. 2). The Phase
boundary between MLC and F changes, as we change the total Zeeman couplings at
a fixed $L_{z}$. Here in Table 1 we represent all different phase boundaries.

Phases | Boundary equation
---|---
KD,CAF | $K_{z}=-K_{xy}+E^{2}_{z}/K_{xy}$
KD,CDW | $K_{z}=K_{xy}$
F,CDW | $K_{z}=-K_{xy}+E_{z}$
F,MLC | $K_{z}=-2K_{xy}-2L_{z}-2E_{z}$
MLC,CDW | $-\frac{3K_{z}}{2}=\left(L_{z}+2K_{xy}-\sqrt{\left(2L_{z}+K_{xy}\right)^{2}+3E^{2}_{z}}\right)$
CAF,MLC | $K_{z}=2\left(K_{xy}-L_{z}\right)$
KD,MLC | $K_{z}=-3K_{xy}-\sqrt{\left(K_{xy}-2L_{z}\right)^{2}+4E^{2}_{z}}$
F,CAF | $K_{xy}=-E_{z}/2$
Table 1: Phase boundary equations as a function of the parameters

##  IV Goldstone modes

The Hamiltonian in the presence of Zeeman term has five different $U(1)$
symmetries coming from $U(1)_{\text{spin}}\otimes U(1)_{\text{valley}}$ for
the two layers and a layer $U(1)\times\mathbb{Z}_{2}$ symmetry. For the layer
diagonal phases, the presence of gapless bulk Goldstone mode is known. The CDW
phase has no gapless bulk mode. In the continuum limits it seems that KD phase
breaks a continuous symmetry but as valley indices are momenta, it breaks
lattice symmetry. Hence in this phase we will have no Goldstone modes. The F
phase has spin wave mode and at long wavelength, it’s gap is proportional to
the Zeeman coupling strength ($E_{z}$). As the CAF phase breaks the spin
$U(1)$ symmetry there will be a pair of gapless neutral modes in the bulk
Murthy et al. (2014).

Next we discuss the new layer coherent phase and its bulk modes. From Eq. 16
one can easily see that the ground state has the two leftover $U(1)$
symmetries defined by operators $\sigma^{3}\otimes\tau^{3}\otimes\gamma^{0}$
and $\sigma^{3}\otimes\tau^{0}\otimes\gamma^{3}$. These operations can be
understood as opposite spin rotations at different valleys or different
layers. In other words, these are relative valley and layer spin twists
respectively. Thus out of five different continuous symmetries, three are
broken by the ground states giving rise to three different Goldstone modes in
the bulk. However, these modes will be neutral as there is a charge gap in the
bulk and these excitations are similar to spin waves. We remind the readers
here that the breaking of the valley part of the symmetry breaks the lattice
$\mathcal{C}_{3}$ rotation about a site. This happens as the $n=0$ manifold
the $\mathbf{K},\mathbf{K}^{\prime}$ of each layer maps to the $A,B$
sublattice of each layerKharitonov (2012a). This means we will count one extra
Goldstone mode in the continuum analysis.

##  V Discussion

In this manuscript, we constructed the Hamiltonian using symmetry principles
without discussing the nature and details of the interaction at the lattice
scale. The model only assumes the lattice interactions are local and thus
their Fourier transform is a function independent of momentum.

In principle we can reproduce the interactions in the continuum model by
projecting the microscopic Hamiltonian to the lowest Landau Level. We present
a simplified example which includes the onsite Hubbard interaction ($U_{1}$),
nearest neighbor interlayer electron density-density interaction ($U_{2}$) and
a intralayer nearest neighbor spin-spin interaction ($J$),

$\displaystyle H_{\text{lat}}=$ $\displaystyle
U_{1}\sum_{\begin{subarray}{c}s_{1},s_{2}\\\
r,L\end{subarray}}\mathrel{\mathop{\mathchar
58\relax}}n_{s_{1},r,L}n_{s_{2},r,L}\mathrel{\mathop{\mathchar
58\relax}}+U_{2}\sum_{\begin{subarray}{c}s_{1},s_{2}\\\ r,r^{\prime}\\\
L_{1}\neq L_{2}\end{subarray}}\mathrel{\mathop{\mathchar
58\relax}}n_{s_{1},r,L_{1}}n_{s_{2},r^{\prime},L_{2}}\mathrel{\mathop{\mathchar
58\relax}}$ $\displaystyle+J\sum_{\begin{subarray}{c}\langle
r,r^{\prime}\rangle,L\end{subarray}}\mathrel{\mathop{\mathchar
58\relax}}\vec{S}_{r,L}\cdot\vec{S}_{r^{\prime},L}\mathrel{\mathop{\mathchar
58\relax}}.$ (19)

Here $n$ is the fermion number operator and $\vec{S}$ is the local spin
operator. Though we assumed ultra-short interactions, adding finite range to
these interactions do not change the symmetry of the continuum model when we
project the hamilitonian in the zero Landau level manifold. We can write the
relation between the continuum parameters in Eq. 3 in terms of the parameters
of Eq. 19 as $K_{xy}\propto-J$, $K_{z}\propto U_{1}$ and $L_{z}\propto
U_{1}/2-U_{2}$.

We would like to emphasise here that the Hamiltonian presented here is just an
example to show that even at the simplest model at the lattice level, we can
achieve the Hamiltonian in Eq. 3. Here we are not concerned with the
values/signs of different parameters $K_{xy},K_{z},L_{z}$ but showing the
phases that is determined by these parameters. Modeling a generic lattice
theory with physically motivated parameters and their values is an interesting
study but out of the scope for the current manuscript. We hope to study a
lattice model of a double layer graphene in the future.

From Eq. 16 we can see that as the states are mixture of
$\vec{K},\vec{K}^{\prime}$ of different layers, near the edge the dispersion
will contain two pairs of particle-like bands and another two pairs of hole-
like bands due to the breaking of the translation symmetry Murthy et al.
(2017). There will be a pair of counter propagating modes only if we have
identical layers near the edge of the system. Near the edge, both the valley
$U(1)$ symmetry and the layer $U(1)$ will be broken generically. Thus the edge
of a double layer graphene will be gapped if the bulk is in MLC phase. As the
states are superposition of two different layers, there will be a drag in the
two terminal measurements at least for finite temperature Titov et al. (2013).
However, we know that there will be two pair of counter propagating modes at
the edge for each layer in the F phase Kharitonov (2012b); Murthy et al.
(2017); Young et al. (2014). Thus by changing the Zeeman energy with respect
to the interaction energies, one can make a transition from MLC to F. This
should show up in the two terminal conductance measurementYoung et al. (2014).
It will also be interesting to measure the lattice scale structure using both
spin resolved and spin unresolved Li et al. (2019b) tunneling electron
microscope to confirm the phase directly.

##  VI Summary and outlook

Here we argued that the continuum limit of the double layered graphene at
$\nu_{1}=\nu_{2}=0$ has a big symmetry group that restricts the interacting
part of the Hamiltonian severely to only three parameters at $n=0$ Landau
Level. Further, we find a candidate ground state using the HF approximation
that breaks the layer $U(1)$ symmetry. We also find a second order phase
transition from MLC to F as a function of Zeeman energy. We argued for a
general system, the edge of the MLC will be gapped. This leads to the possible
experiment to find two terminal conductance that will change when we go from
the MLC to F.

This study is just the beginning of understanding the double layer graphene
ground state in the quantum Hall regime. It was previously shown that the
phase transition from CAF to F connects the bulk gapless modes of the CAF to
the gapless edge modes of F Murthy et al. (2014). We hope to study the edge
theory of double layer graphene in future to answer the question of the phase
transition from MLC to F.

It has been shown that if we have finite range interactions in monolayer
graphene then we can have co-existence of phases Das et al. (2021). Similarly,
in a lattice model co-existence can also be shown by doing HF calculation in
the lattice limit Das (2020). This explains the experimental results Li et al.
(2019c), where bond order was observed using Scanning Tunnelling Microscope.
This question may also be important in double layer graphene case, as we might
have a similar coexistence. To understand that possibility one needs to study
that lattice Hamiltonian similar to the one mentioned in Eq. 19. Furthermore,
this theory can be used to explore the phase diagrams at other filling
fractions in the parameter space of $K_{xy},K_{z}$ and $L_{z}$ similar to the
MLG case Sodemann and MacDonald (2014). There is also a surge of interest in
understanding the BCS/BEC condensation Liu et al. (2020b); Eisenstein (2014);
Liu et al. (2017); Li et al. (2017) in double layer graphene systems. As
previously mentioned this state breaks the layer $U(1)$ symmetry just like a
superfluid state. At low enough temperatures these excitons can form a
superfluid state where the interaction between the electron and hole can be
tuned by tuning the $L_{z}$ parameter (which depends on the interlayer
separation $d$).

###### Acknowledgements.

We thank G. Murthy, E. Berg, and J. S. Hofmann for useful discussions. AS
would like to thank the US-Israel Binational Science Foundation for its
support via grant no. 2016130 and the University of Kentucky Center for
Computational Sciences and Information Technology Services Research Computing
for their support and use of the Lipscomb Compute Cluster and associated
research computing resources. AD was supported by the German-Israeli
Foundation (GIF) Grant No. I-1505-303.10/2019 and the Minerva Foundation. AD
also thanks Weizmann Institute of Science, Israel Deans fellowship, and Israel
planning and budgeting committee for financial support.

##  Appendix A Details of the technique

The interacting Hamiltonian in Eq. 3 can be written in a simplified form as,

$\displaystyle H_{int}=\frac{\pi
l^{2}}{A}\sum_{\begin{subarray}{c}\mathbf{q},k_{1},k_{2}\\\
a,b,c,d\end{subarray}}e^{-iq_{x}(k_{1}-k_{2}-q_{y})l^{2}-\frac{q^{2}l^{2}}{2}}V_{a,b,c,d}\mathrel{\mathop{\mathchar
58\relax}}c^{\dagger}_{k_{1}-q_{y},a}c_{k_{1},b}c^{\dagger}_{k_{2}+q_{y},c}c_{k_{2},b}\mathrel{\mathop{\mathchar
58\relax}}$ (20)

where

$\displaystyle V_{a,b,c,d}=K_{xy}$
$\displaystyle\sum_{i=1,2}\left(\sigma^{0}\otimes\tau^{i}\otimes
P_{L}\right)_{a,b}\left(\sigma^{0}\otimes\tau^{i}\otimes
P_{L}\right)_{c,d}+K_{z}\left(\sigma^{0}\otimes\tau^{3}\otimes
P_{L}\right)_{a,b}\left(\sigma^{0}\otimes\tau^{3}\otimes P_{L}\right)_{c,d}$
$\displaystyle+L_{z}\left(\sigma^{0}\otimes\tau^{0}\otimes\gamma^{3}\right)_{a,b}\left(\sigma^{0}\otimes\tau^{0}\otimes\gamma^{3}\right)_{c,d}.$
(21)

To calculate the total HF energy $E=\langle HF|H_{0}+H_{int}|HF\rangle$ we
write the average of the four-fermion term that arises in the interacting
Hamiltonian $H_{int}$ as

$\langle
HF|c^{\dagger}_{k_{1}-q_{y},a}c^{\phantom{\dagger}}_{k_{1},b}c^{\dagger}_{k_{2}+q_{y},c}c^{\phantom{\dagger}}_{k_{2},d}|HF\rangle=\Delta_{b,a}\Delta_{d,c}\delta_{q_{y},0}-\Delta_{d,a}\Delta_{b,c}\delta_{q_{y},k_{1}-k_{2}}.$
(22)

The first term gives the Hartree term and the second term is the Fock term.
Using this, we can calculate the energy from the electron-electron interaction
given by $E_{int}=\langle HF|H_{int}|HF\rangle$

$\displaystyle E_{int}=$ $\displaystyle\frac{\pi
l^{2}}{A}\sum_{\begin{subarray}{c}a,b\\\
c,d\end{subarray}}\sum_{\begin{subarray}{c}\mathbf{q}\\\
k_{1},k_{2}\end{subarray}}e^{-iq_{x}(k_{1}-k_{2}-q_{y})l^{2}-\frac{q^{2}l^{2}}{2}}V_{a,b,c,d}\left(\Delta_{b,a}\Delta_{d,c}\delta_{q_{y},0}-\Delta_{d,a}\Delta_{b,c}\delta_{k_{1}-q_{y},k_{2}}\right)$
$\displaystyle=$ $\displaystyle\frac{\pi
l^{2}}{A}\sum_{\begin{subarray}{c}a,b,c,d\end{subarray}}V_{a,b,c,d}\left(\sum_{\begin{subarray}{c}q_{x}\\\
k_{1},k_{2}\end{subarray}}e^{-iq_{x}(k_{1}-k_{2})l^{2}-\frac{q_{x}^{2}l^{2}}{2}}\Delta_{b,a}\Delta_{d,c}-\sum_{\begin{subarray}{c}\mathbf{q},k_{1}\end{subarray}}e^{-\frac{q^{2}l^{2}}{2}}\Delta_{d,a}\Delta_{b,c}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2N_{\Phi}}\sum_{\begin{subarray}{c}i\\\
a,b,c,d\end{subarray}}V_{a,b,c,d}\left(N_{\Phi}^{2}\Delta_{b,a}\Delta_{d,c}-\frac{N_{\Phi}A}{(2\pi)^{2}}\int_{\mathbf{q}}d\mathbf{q}e^{-\frac{q^{2}l^{2}}{2}}\Delta_{d,a}\Delta_{b,c}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2N_{\Phi}}\sum_{\begin{subarray}{c}a,b,c,d\end{subarray}}V_{a,b,c,d}\left(N_{\Phi}^{2}\Delta_{b,a}\Delta_{d,c}-\frac{N_{\Phi}A}{2\pi
l^{2}}\Delta_{d,a}\Delta_{b,c}\right)$ (23a) $\displaystyle E_{int}=$
$\displaystyle\frac{N_{\Phi}}{2}\sum_{\begin{subarray}{c}i\\\
a,b,c,d\end{subarray}}V_{a,b,c,d}\left(\Delta_{b,a}\Delta_{d,c}-\Delta_{d,a}\Delta_{b,c}\right)$
(23b)

where $A$ is the area of the system and $N_{\Phi}=A/(2\pi l^{2})$ is the
number of guiding centers in the system. Hence the total energy of the system
per guiding center is

$\displaystyle\frac{E}{N_{\Phi}}=E_{z}(\sigma_{3}\otimes\tau_{0}\otimes\gamma_{0})_{ab}\Delta_{b,a}+\frac{1}{2}\sum_{a,b,c,d}V_{a,b,c,d}\left(\Delta_{b,a}\Delta_{d,c}-\Delta_{d,a}\Delta_{b,c}\right).$
(24)

The first term is the Zeeman contribution and the second term comes from the
elecron-electron interaction.

## References

  * Klitzing et al. (1980) K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980), URL https://link.aps.org/doi/10.1103/PhysRevLett.45.494.
  * Chakraborty and Pietiläinen (1987) T. Chakraborty and P. Pietiläinen, Phys. Rev. Lett. 59, 2784 (1987), URL https://link.aps.org/doi/10.1103/PhysRevLett.59.2784.
  * Halperin (1983) B. I. Halperin, Helv. Phys. Acta 56, 75 (1983).
  * Jain (2007) J. K. Jain, _Composite Fermions_ (Cambridge University Press, 2007).
  * Scarola and Jain (2001) V. W. Scarola and J. K. Jain, Phys. Rev. B 64, 085313 (2001), URL https://link.aps.org/doi/10.1103/PhysRevB.64.085313.
  * Eisenstein (2014) J. Eisenstein, Annual Review of Condensed Matter Physics 5, 159 (2014), eprint https://doi.org/10.1146/annurev-conmatphys-031113-133832, URL https://doi.org/10.1146/annurev-conmatphys-031113-133832.
  * Kellogg et al. (2004) M. Kellogg, J. P. Eisenstein, L. N. Pfeiffer, and K. W. West, Phys. Rev. Lett. 93, 036801 (2004), URL https://link.aps.org/doi/10.1103/PhysRevLett.93.036801.
  * Gramila et al. (1991) T. J. Gramila, J. P. Eisenstein, A. H. MacDonald, L. N. Pfeiffer, and K. W. West, Phys. Rev. Lett. 66, 1216 (1991), URL https://link.aps.org/doi/10.1103/PhysRevLett.66.1216.
  * Eisenstein and MacDonald (2004) J. P. Eisenstein and A. H. MacDonald, Nature 432, 691 (2004), ISSN 1476-4687, URL https://doi.org/10.1038/nature03081.
  * Novoselov et al. (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005), ISSN 1476-4687, URL https://doi.org/10.1038/nature04233.
  * Zhang et al. (2005) Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim, Nature 438, 201 (2005), ISSN 1476-4687, URL https://doi.org/10.1038/nature04235.
  * Li et al. (2016) J. I. A. Li, T. Taniguchi, K. Watanabe, J. Hone, A. Levchenko, and C. R. Dean, Phys. Rev. Lett. 117, 046802 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.046802.
  * Li et al. (2019a) J. I. A. Li, Q. Shi, Y. Zeng, K. Watanabe, T. Taniguchi, J. Hone, and C. R. Dean, Nature Physics 15, 898 (2019a), ISSN 1745-2481, URL https://doi.org/10.1038/s41567-019-0547-z.
  * Liu et al. (2019) X. Liu, Z. Hao, K. Watanabe, T. Taniguchi, B. I. Halperin, and P. Kim, Nature Physics 15, 893 (2019), ISSN 1745-2481, URL https://doi.org/10.1038/s41567-019-0546-0.
  * Gorbachev et al. (2012) R. V. Gorbachev, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. Tudorovskiy, I. V. Grigorieva, A. H. MacDonald, S. V. Morozov, K. Watanabe, T. Taniguchi, et al., Nature Physics 8, 896 (2012), ISSN 1745-2481, URL https://doi.org/10.1038/nphys2441.
  * Titov et al. (2013) M. Titov, R. V. Gorbachev, B. N. Narozhny, T. Tudorovskiy, M. Schütt, P. M. Ostrovsky, I. V. Gornyi, A. D. Mirlin, M. I. Katsnelson, K. S. Novoselov, et al., Phys. Rev. Lett. 111, 166601 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.111.166601.
  * Liu et al. (2017) X. Liu, K. Watanabe, T. Taniguchi, B. I. Halperin, and P. Kim, Nature Physics 13, 746 (2017), ISSN 1745-2481, URL https://doi.org/10.1038/nphys4116.
  * Li et al. (2017) J. I. A. Li, T. Taniguchi, K. Watanabe, J. Hone, and C. R. Dean, Nature Physics 13, 751 (2017), ISSN 1745-2481, URL https://doi.org/10.1038/nphys4140.
  * Britnell et al. (2012) L. Britnell, R. V. Gorbachev, R. Jalil, B. D. Belle, F. Schedin, M. I. Katsnelson, L. Eaves, S. V. Morozov, A. S. Mayorov, N. M. R. Peres, et al., Nano Letters 12, 1707 (2012), ISSN 1530-6984, URL https://doi.org/10.1021/nl3002205.
  * Lee et al. (2016) K. Lee, J. Xue, D. C. Dillen, K. Watanabe, T. Taniguchi, and E. Tutuc, Phys. Rev. Lett. 117, 046803 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.046803.
  * Liu et al. (2020a) X. Liu, Z. Wang, K. Watanabe, T. Taniguchi, O. Vafek, and J. I. A. Li, _Tuning electron correlation in magic-angle twisted bilayer graphene using coulomb screening_ (2020a), eprint arXiv:2003.11072.
  * Tse et al. (2007) W.-K. Tse, B. Y.-K. Hu, and S. Das Sarma, Phys. Rev. B 76, 081401 (2007), URL https://link.aps.org/doi/10.1103/PhysRevB.76.081401.
  * Katsnelson (2011) M. I. Katsnelson, Phys. Rev. B 84, 041407 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.84.041407.
  * Peres et al. (2011) N. M. R. Peres, J. M. B. L. dos Santos, and A. H. C. Neto, EPL (Europhysics Letters) 95, 18001 (2011), URL https://doi.org/10.1209/0295-5075/95/18001.
  * Narozhny et al. (2012) B. N. Narozhny, M. Titov, I. V. Gornyi, and P. M. Ostrovsky, Phys. Rev. B 85, 195421 (2012), URL https://link.aps.org/doi/10.1103/PhysRevB.85.195421.
  * Kharitonov and Efetov (2010) M. Y. Kharitonov and K. B. Efetov, Semiconductor Science and Technology 25, 034004 (2010), URL https://doi.org/10.1088/0268-1242/25/3/034004.
  * Kim et al. (2011) S. Kim, I. Jo, J. Nah, Z. Yao, S. K. Banerjee, and E. Tutuc, Phys. Rev. B 83, 161401 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.83.161401.
  * Min et al. (2008) H. Min, R. Bistritzer, J.-J. Su, and A. H. MacDonald, Phys. Rev. B 78, 121401 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.78.121401.
  * Alicea and Fisher (2006) J. Alicea and M. P. A. Fisher, Phys. Rev. B 74, 075422 (2006), URL https://link.aps.org/doi/10.1103/PhysRevB.74.075422.
  * Kharitonov (2012a) M. Kharitonov, Phys. Rev. B 85, 155439 (2012a), URL https://link.aps.org/doi/10.1103/PhysRevB.85.155439.
  * Herbut (2007a) I. F. Herbut, Phys. Rev. B 75, 165411 (2007a), URL https://link.aps.org/doi/10.1103/PhysRevB.75.165411.
  * Herbut (2007b) I. F. Herbut, Phys. Rev. B 76, 085432 (2007b), URL https://link.aps.org/doi/10.1103/PhysRevB.76.085432.
  * Kharitonov (2012b) M. Kharitonov, Phys. Rev. B 86, 075450 (2012b), URL https://link.aps.org/doi/10.1103/PhysRevB.86.075450.
  * Murthy et al. (2014) G. Murthy, E. Shimshoni, and H. A. Fertig, Phys. Rev. B 90, 241410 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.90.241410.
  * Murthy et al. (2016) G. Murthy, E. Shimshoni, and H. A. Fertig, Phys. Rev. B 93, 045105 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.045105.
  * Knothe and Jolicoeur (2015) A. Knothe and T. Jolicoeur, Phys. Rev. B 92, 165110 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.92.165110.
  * Young et al. (2014) A. F. Young, J. D. Sanchez-Yamagishi, B. Hunt, S. H. Choi, K. Watanabe, T. Taniguchi, R. C. Ashoori, and P. Jarillo-Herrero, Nature 505, 528 (2014), ISSN 1476-4687, URL https://doi.org/10.1038/nature12800.
  * Sodemann and MacDonald (2014) I. Sodemann and A. H. MacDonald, Phys. Rev. Lett. 112, 126804 (2014), URL https://link.aps.org/doi/10.1103/PhysRevLett.112.126804.
  * Wen and Zee (1992) X.-G. Wen and A. Zee, Phys. Rev. Lett. 69, 1811 (1992), URL https://link.aps.org/doi/10.1103/PhysRevLett.69.1811.
  * Murthy et al. (2017) G. Murthy, E. Shimshoni, and H. A. Fertig, Phys. Rev. B 96, 245125 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.245125.
  * Li et al. (2019b) S.-Y. Li, Y. Zhang, L.-J. Yin, and L. He, Phys. Rev. B 100, 085437 (2019b), URL https://link.aps.org/doi/10.1103/PhysRevB.100.085437.
  * Das et al. (2021) A. Das, R. K. Kaul, and G. Murthy, _Coexistence of canted antiferromagnetism and bond-order in $\nu=0$ graphene_ (2021), eprint 2109.07515.
  * Das (2020) A. Das, Ph.D. thesis, University of Kentucky, Lexington (2020), URL https://uknowledge.uky.edu/physastron_etds/70/.
  * Li et al. (2019c) S.-Y. Li, Y. Zhang, L.-J. Yin, and L. He, Phys. Rev. B 100, 085437 (2019c), URL https://link.aps.org/doi/10.1103/PhysRevB.100.085437.
  * Liu et al. (2020b) X. Liu, J. I. A. Li, K. Watanabe, T. Taniguchi, J. Hone, B. I. Halperin, P. Kim, and C. R. Dean, _Crossover between strongly-coupled and weakly-coupled exciton superfluids_ (2020b), eprint 2012.05916.


# Nonlinear Hawkes Process with Gaussian Process Self Effects

Noa Malem-Shinitski
Institute of Mathematics
University of Potsdam
Potsdam, Germany
<EMAIL_ADDRESS>
&Cèsar Ojeda
Artificial Intelligence Group
Technische Universität Berlin
Berlin, Germany
&Manfred Opper
Artificial Intelligence Group
Technische Universität Berlin
Berlin, Germany

###### Abstract

Traditionally, Hawkes processes are used to model time–continuous point
processes with history dependence. Here we propose an extended model where the
self–effects are of both excitatory and inhibitory type and follow a Gaussian
Process. Whereas previous work either relies on a less flexible
parameterization of the model, or requires a large amount of data, our
formulation allows for both a flexible model and learning when data are
scarce. We continue the line of work of Bayesian inference for Hawkes
processes, and our approach dispenses with the necessity of estimating a
branching structure for the posterior, as we perform inference on an
aggregated sum of Gaussian Processes. Efficient approximate Bayesian inference
is achieved via data augmentation, and we describe a mean–field variational
inference approach to learn the model parameters. To demonstrate the
flexibility of the model we apply our methodology on data from three different
domains and compare it to previously reported results.

##  1 Introduction

Sequences of self exciting, or inhibiting, temporal events are frequent
footmarks of natural phenomena: Earthquakes are known to be temporally
clustered as aftershocks are commonly triggered following the occurrence of a
main event (Ogata, 1988); in social networks, the propagation of news can be
modeled in terms of information cascades over the edges of a graph (Zhao et
al., 2015); and in neuronal activity, the occurrence of one spike may increase
or decrease the probability of the occurrence of the next spike over some time
period (Dayan and Abbott, 2001).

Traditionally, sequences of events in continuous time are modeled by Point
processes, of which Cox processes (Cox, 1955), or doubly stochastic processes,
use a stochastic process for the intensity function, which depends only on
time and is not effected by the occurrences of the events. The Hawkes process
(Hawkes and Oakes, 1974) extends the Cox process to capture phenomena in which
the past events affects future arrivals, by introducing a memory dependence
via a memory kernel. When incorporating dependence of the process on its own
history, due to the superposition theorem of point process, new events will
depend on either an exogenous rate, which is independent of the history, or an
endogenous rate from past arrivals. This results in a branching structure,
where new events that originate from previous events can be seen as "children"
of the past events.

Originally, the dependence on the history in the Hawkes process is assumed to
be self excitatory, and the memory kernel is parameterized by an exponential
or power law decay, which results in a model with low flexibility.
Furthermore, assuming only excitatory relation between the events, does not
hold for other phenomena we wish to model. For example, inhibitory effects
between neurons (Maffei et al., 2004), and even self–inhibition (Smith and
Jahr, 2002), are crucial for regulating the neuronal activity. Thus, the
memory kernel should also include inhibitory relations between the events and
by doing so the intensity may become negative. To ensure that the intensity
function is non–negative, a nonlinear link function is applied on the memory
kernel, and the resulting model is often referred to as a Nonlinear Hawkes
process (Brémaud and Massoulié, 1996; Zhu, 2013; Truccolo, 2016).

In this work we present a Nonlinear Hawkes process with Gaussian Process
Self–effects (NH-GPS) which extends the class of Nonlinear Hawkes processes.
We choose a non–parametric approach which avoids the limiting parameterization
of the memory kernel and the background rate. We assume a Gaussian Process
(GP) prior on the exogenous events intensity and on the memory kernel, which
allows also for an inhibitory effect between the events. To ensure that the
intensity function is non–negative we use the Sigmoid link function. This
modeling approach is not only descriptive, but also allows us to obtain a fast
inference procedure. The history of self–effects defines an aggregated
Gaussian process, and we perform the inference directly on this aggregation
rather than obtaining a posterior over each self effect.

While highly flexible approaches to modeling the intensity function of
nonlinear Hawkes processes have been presented before, they mainly rely deep
neural network solutions (Jia and Benson, 2019; Xiao et al., 2017). These
approaches are hindered by the necessity of large datasets. However, our
methodology retains the modeling flexibility due to the non–parametric nature
of Gaussian processes, while being able to perform well when data are scarce.

Our main contributions are the following

  * •

We present a new flexible non–parametric approach to the nonlinear Hawkes
process, which allows us to perform inference when data are scarce.

  * •

We derive an efficient variational inference scheme to the model.

##### Outline

In Section 2 we discuss related work and describe how it relates to ours. In
Section 3 we describe the NH–GPS model and emphasize how its structure allows
for efficient Bayesian inference. In Section 4 we describe the augmentation
scheme and derive the mean–field variational inference algorithm. In Section 5
we present the results of our experiments both on synthetic data and different
real world examples, and compare to existing work when possible. In Section 6
we conclude by discussing our work and future research directions.

##  2 Related Work

Bayesian approaches to Cox Processes model the intensity with a Gaussian
process prior, which is then passed through a link function to ensure its
positivity. A common choice of the link function is the exponential or the
quadratic functions (Hensman et al., 2015; Lloyd et al., 2015). Another
choice, which is more relevant to our wok, is the sigmoid link function,
resulting in the _sigmoidal Gaussian Cox process_. Inference in this model was
first done via Markov Chain Monte Carlo (Adams et al., 2009) and later via
variational inference (Donner and Opper, 2018a).

As for the Hawkes process, first attempts to perform Bayesian inference relied
on the definition in terms of a marked Poisson cluster process and identifying
the branching structure of the self–excitation (Rasmussen, 2013). Alternative
approaches to inference in temporal point processes and Hawkes processes
include directly maximizing the Likelihood function via stochastic gradient
descent (Du et al., 2016) or re–modeling the intensity function as a Recurrent
Neural Network (RNN) (Mei and Eisner, 2017; Jing and Smola, 2017; Zuo et al.,
2020). The latter often require training over very large datasets to produce
reliable results.

Other than RNN, another highly flexible approach to estimating the intensity
function of the Hawkes process relies on GP priors (Zhang et al., 2018; Zhou
et al., 2019; Zhang et al., 2020). Recent adaptation of this approach is the
model described in Zhou et al. (2020). Similarly to the model described in our
work, the authors avoid the limiting parameterization of the memory kernel by
using GPs. Differently to our work, Zhou et al. (2020) remain in the linear
Hawkes process regime and assume that the effects of past events are assumed
to be only excitatory, whereas our approach allows both excitatory and
inhibitory effects.

Noticeable variations of the nonlinear Hawkes process are the Isotonic Hawkes
Process (Wang et al., 2016) and the Mutually Regressive Point Process
(Apostolopoulou et al., 2019). In the first variation, although the model
allows for both inhibitory and excitatory self–effects, there is a relation of
either–or between them. Thus, the entire history may be either inhibiting or
exciting. In contrast, the GP prior in our approach allows for the self-
effects to change over time. Hence, one particular event may have inhibiting
effect on the next one, but have an exciting effect on the one after that.

The second variation is more closely related to our work. The Mutually
Regressive Point Process is designed to model neuronal spike trains. In this
work, the classical self–excitatory Hawkes Process intensity function is
augmented by a probability term. This term induces inhibition when it is close
to zero. In a sense, this model includes two memory kernels – one excitatory
only which appears in the intensity function and another which can also induce
inhibition in the augmenting probability term. In the current work, we achieve
such flexibility of the self–effects in a simpler fashion by assuming the GP
prior on the self effects. As mentioned before this also allows for the type
of effect to change over time, which dows not appear in the work of
Apostolopoulou et al. (2019).

##  3 Proposed Model

###  3.1 Classical Hawkes Process

Let $\mathcal{T}_{T}=[0,t]\in\mathbb{R}$. We define the counting measure
$N(\mathcal{T}_{t})$ as the number of arrivals in the sequence
$\mathcal{H}_{t}=\\{T_{1},...,T_{N(\mathcal{T}_{t})}:T_{i}\in\mathcal{T}_{t}\wedge
T_{i-1}<T_{i}\\}$ where $\mathcal{H}_{t}$ defines the history of the process
until time $t$, and $T_{i}$ corresponds to the time of arrival $i$. For a
temporal point process, the counting measure $N(\cdot)$ has an associated
intensity defined as

$\displaystyle\Lambda(t)=\lim_{\Delta t\rightarrow
0}\frac{\mathbb{E}[N(\mathcal{T}_{t+\Delta
t})-N(\mathcal{T}_{t})|\mathcal{H}_{t}]}{\Delta t}\mbox{.}$

The intensity function may depend on the history of the process. An example of
such a process is the Hawkes process, or self exciting point process,
(Kingman, 1993) which defines self excitations (Daley and Vere-Jones, 2007)
around exogenous events.

Following Hawkes and Oakes (1974), the intensity of the Hawkes process is
defined by

$\displaystyle\Lambda(t|\mathcal{H}_{t})=s\left(t\right)+\sum_{t_{n}<t}g\left(t-t_{n}\right)\mbox{,}$
(1)

where $s(t)$ is the base intensity of exogenous arrivals and
$g\left(t-t_{n}\right)$ is the memory kernel defining the change in the
excitation value for each arrival. In the classical Hawkes process, only
excitations are allowed and the memory kernel is usually of the form
$g(t-t_{n})=\beta e^{-\alpha(t-t_{n})}$ for an exponentially decaying memory.

###  3.2 Nonlinear Hawkes process with Gaussian Process Self–Effects

In the classical Hawkes process, the memory kernel $g$ in Equation 1 must be
non–negative, to prevent the intensity function from being negative. As a
result the history of the model has only excitatory effect on future events.
We are interested in a model that includes inhibition between events, and we
release the constraint over $g$ so it can be negative, and define the
following nonlinear intensity function

$\displaystyle\Lambda$
$\displaystyle\left(t\right)=\lambda\sigma\left(\phi(t)\right)$ (2)
$\displaystyle\sigma\left(\phi(t)\right)$
$\displaystyle=\frac{1}{1+\exp{\left(-\phi\left(t\right)\right)}}$ (3)
$\displaystyle\phi(t)=s\left(t\right)+$
$\displaystyle\sum_{t_{n}<t}g\left(t-t_{n}\right)\exp\left(-\alpha\left(t-t_{n}\right)\right).$
(4)

Here, we choose the sigmoid function to ensure that the intensity function
$\Lambda\left(\cdot\right)$ is non–negative. $\lambda$ is the intensity bound
and we refer to $\phi\left(\cdot\right)$ as the linear intensity function.

We explicitly add the exponential decay to enforce the forgetting constraint
which is essential for most realistic processes. Although we choose here a
specific parameterization of the memory decay, one can choose other forms of
memory decay with minimal adaptation to the learning procedure of the model
parameters.

To maximize the flexibility of the model we avoid any specific
parameterization of the background rate or the memory kernel. Thus, rather
than specifying a functional form for $s\left(t\right)$ and $g\left(t\right)$,
we assume the following GP priors

$\displaystyle s$ $\displaystyle\sim GP\left(0,K^{s}\right)$ (5)
$\displaystyle g$ $\displaystyle\sim GP\left(0,K^{g}\right)$ (6)
$\displaystyle K$
$\displaystyle\left(t_{1},t_{2}\right)=a\cdot\exp{\left(-\frac{\|t_{1}-t_{2}\|^{2}}{\sigma^{2}}\right)}.$
(7)

In this work we use the Radial Basis Function (RBF) kernel for the GP priors.
This choice is not a constraint of the model – one can choose any other
kernel, and it will not effect the augmentation and inference processes
described bellow.

Finally, we assume a prior distribution also on the upper intensity bound

$\displaystyle\lambda\sim Gamma\left(\alpha_{0},\beta_{0}\right).$

and we identify the hyperparameters of the model as
$\\{\sigma_{g},a_{g},\alpha,\sigma_{s},a_{s}\\}$.

In this work we propose Bayesian inference for fitting the model to data. Due
to the non–linearity over $\phi\left(\cdot\right)$ we are no longer able to
easily utilize the branching structure of the Hawkes process which allowed for
the estimation of $s\left(\cdot\right)$ and $g\left(\cdot\right)$ (Rasmussen,
2013; Zhou et al., 2020). Thus, a natural solution is to perform the inference
directly on $\phi\left(\cdot\right)$.

Next, we identify the prior over the entire linear intensity
$p\left(\phi\right)$. From Equation 4 we see that the linear intensity
function $\phi$ is nothing but the sum of GPs, and as such it is also a GP

$\displaystyle\phi\sim GP\left(0,\tilde{K}\right)$
$\displaystyle\tilde{K}_{lk}=K^{s}_{lk}+\sum_{t_{i}<t_{l}}\sum_{t_{j}<t_{k}}K^{g}_{t_{l}-t_{i},t_{k}-t_{j}}\exp\left(-\alpha\left(t_{l}-t_{i}+t_{k}-t_{j}\right)\right).$

####  3.2.1 Multivariate Model

We propose an extension of the model to multiple dimensions. This is useful in
applications where different types of events are observed, or the events
originate from different processes that effect each other. We define an
$R$–dimensional point process with intensity in dimension $r$

$\displaystyle\Lambda^{r}\left(t\right)=\lambda^{r}\sigma\left(\phi^{r}\left(t\right)\right)$
$\displaystyle\phi^{r}\left(t\right)=s^{r}\left(t\right)+\sum_{m=1}^{R}\sum_{t^{m}_{n}<t}g_{r,m}\left(t-t^{m}_{n}\right)\exp\left(-\alpha_{r,m}\left(t-t^{m}_{n}\right)\right)$

where $t_{n}^{m}$ is the time of event number $n$ of type $m$. We assume that
every dimension has its own intensity bound $\lambda^{k}$ and background rate
$s^{k}\left(\cdot\right)$. The different dimension interact with each other
via the self effects term. $g_{k,l}\left(\cdot\right)$ defined the effects of
events of type $l$ on events of type $k$. As in the univariate case, this
effect may change over time.

Given the observations, the different dimensions are independent of each other
and we can learn their parameters separately. Thus, in the following section
we present the inference for the univariate model, and the extension to the
multivariate case is straight forward.

##  4 Inference

Conditioned on the intensity function $\Lambda\left(\cdot\right)$, the
likelihood of observations $\\{t_{1},...t_{n}\\}$ from a Hawkes process is
(Daley and Vere-Jones, 2003)

$\displaystyle\ell\left(\\{t_{1},...t_{n}\\}|\Lambda\left(\cdot\right)\right)=\exp\left\\{-\int^{T}_{0}\Lambda\left(t^{\prime}\right)dt^{\prime}\right\\}\prod^{N}_{i=1}\Lambda(t_{i}).$

Looking at the likelihood defined above, Equations 4 and 4, implementing
Bayesian inference for the model is not straightforward, due to the non-
conjugate structure of the likelihood and prior. Similarly to previous work on
Cox and Hawkes processes (Donner and Opper, 2018a; Apostolopoulou et al.,
2019; Zhou et al., 2020), we augment the model with auxiliary variables, which
leads to a conditionally conjugated model with closed form solutions for Gibbs
sampler and variational inference.

###  4.1 Model Augmentation

The first step we take in treating the likelihood function is using the Pólya-
Gamma (PG) augmentation scheme. Following Theorem $1$ in Polson et al. (2013),
we can rewrite the nonlinear intensity function as

$\displaystyle\sigma\left(\phi\left(t\right)\right)=\int_{0}^{\infty}e^{f\left(w,\phi\left(t\right)\right)}PG\left(w;1,0\right)dw$
(8) $\displaystyle
f\left(w,\phi\left(t\right)\right)=-\frac{\phi\left(t\right)^{2}w}{2}+\frac{\phi\left(t\right)}{2}-\ln{2}.$
(9)

As we augment each observation with a variable $w_{n}$ from a PG distribution,
the joint likelihood of the observed events $\\{t_{n}\\}$ and PG variables
$\\{w_{n}\\}$ is

$\displaystyle
p\left(\\{t_{n}\\}_{n=1}^{N},\\{w_{n}\\}_{n=1}^{N}|\phi,\lambda\right)=$ (10)
$\displaystyle\exp\left(-\int^{T}_{0}\lambda\sigma\left(\phi\left(t\right)\right)dt\right)\cdot\prod^{N}_{n=1}\lambda
e^{f\left(w_{n},t_{n}\right)}PG\left(w_{n};1,0\right)$

with

$\displaystyle\exp\left\\{-\int^{T}_{0}\lambda\sigma\left(\phi\left(t\right)\right)dt\right\\}=$
(11) $\displaystyle\exp\left(-\int_{0}^{T}\int_{0}^{\infty}\lambda
PG\left(w;1,0\right)\left(1-e^{f\left(w,-\phi\left(t\right)\right)}\right)dwdt\right).$

Where we used $\sigma(t)=1-\sigma(-t)$.

Next, we utilize the Campbell’s theorem (Kingman, 1993) which states that for
a Poisson process $\Pi$ with intensity $\varphi$

$\displaystyle\mathbb{E}_{\varphi}$
$\displaystyle\left(\prod_{x\in\Pi}\exp\left(h\left(x\right)\right)\right)=$
$\displaystyle\exp\left(-\int\left(1-\exp\left(h\left(x\right)\right)\right)\varphi\left(x\right)dx\right).$

Looking at Equation 11 we identify $x=\left(t,w\right)$ and
$\varphi\left(t,w\right)=\lambda PG\left(w|1,0\right)$ is the intensity of a
marked Poisson process in $\mathcal{T}$ with marks $w\sim PG\left(0,1\right)$.
Further, we determine $h\left(x\right)=f\left(w,-\phi\left(t\right)\right)$.
We can now rewrite the exponential in Equation 10 as

$\displaystyle\exp\left\\{-\int^{T}_{0}\lambda\sigma\left(\phi\left(t\right)\right)dt\right\\}=\mathbb{E}_{\varphi}\left(\prod_{m=1}^{M}e^{f\left(\hat{w}_{m},\hat{t}_{m}\right)}\right)$
(12)

for realizations $\\{\hat{t}_{m},\hat{w}_{m}\\}_{m=1}^{M}$.

We substitute Equation 12 into Equation 10 which results in the full augmented
likelihood. Given the prior distribution over $\phi$ and $\lambda$, we can now
write the model’s posterior distribution as

$\displaystyle
p\left(\\{\hat{t}_{m},\hat{w}_{m}\\},\\{w_{n}\\},\phi,\lambda|\\{t_{n}\\}\right)\propto\exp\left(-\lambda
T\right)\times$ (13) $\displaystyle\prod_{m=1}^{M}\lambda
e^{f\left(\hat{w}_{m},-\phi\left(\hat{t}_{m}\right)\right)}PG\left(w_{m};1,0\right)\times$
$\displaystyle\prod_{n=1}^{N}\lambda
e^{f\left(w_{n},\phi\left(t_{n}\right)\right)}PG\left(w_{n};1,0\right)\times
p\left(\phi\right)p\left(\lambda\right).$

To summarize, we augment the model with two sets of variables – the PG
variables $\\{w_{n}\\}$ which augment the actual realizations and the tuples
$\\{\hat{t}_{m},\hat{w}_{m}\\}$ which are the realizations and marks of the
auxiliary marked Poisson process.

As mentioned above, we intend to learn directly the linear intensity function
$\phi\left(\cdot\right)$. This allows us to utilize the efficient mean–field
variational inference previously introduced in Donner and Opper (2018a) and
Donner and Opper (2018b). Next, we go through the steps of the algorithm, and
we refer the reader to the two papers mentioned above for further details. As
a baseline we compare the performance of the variational inference algorithm
to a Gibbs sampler. The details of the Gibbs sampler can be found in the
Supplementary Material.

###  4.2 Variational Inference

In variational inference (Jordan et al., 1999; Bishop, 2006) we define a
tractable distribution family and adapt it to approximate the posterior by
maximizing the lower bound $\mathcal{L}(Q)$ defined below. This procedure
minimizes the Kullback–Leibler divergence between the unknown posterior and
the proposed approximating distribution. The posterior density is approximated
by

$\displaystyle
p\left(\\{\hat{t}_{m},\hat{w}_{m}\\},\\{w_{n}\\},\phi,\lambda|\\{t_{n}\\}\right)$
$\displaystyle\approx
q_{1}\left(\phi,\lambda\right)q_{2}\left(\\{w_{n}\\}_{n=1}^{N},\\{\hat{t}_{m},\hat{w}_{m}\\}_{m=1}^{M}\right).$

This leads to the following lower bound on the evidence

$\displaystyle\mathcal{L}(Q)=\mathbb{E}_{Q}\left[\log\left\\{\frac{p\left(\\{\hat{t}_{m},\hat{w}_{m}\\},\\{w_{n}\\},\phi,\lambda|\\{t_{n}\\}\right)}{q_{1}\left(\phi,\lambda\right)q_{2}\left(\\{w_{n}\\}_{n=1}^{N},\\{\hat{t}_{m},\hat{w}_{m}\\}_{m=1}^{M}\right)}\right\\}\right]\mbox{.}$

Here $Q$ refers to the probability measure of the variational posterior. We
can maximize the bound by alternating the maximization over each of the
factors (Bishop, 2006). The optimal solution for each factor is

$\displaystyle\log q^{*}_{1}\left(\phi,\lambda\right)=$ (14)
$\displaystyle\mathbb{E}_{q_{2}\left(\\{w_{n}\\}_{n=1}^{N},\\{\hat{t}_{m},\hat{w}_{m}\\}_{m=1}^{M}\right)}[\log
P(\\{\hat{t}_{m},\hat{w}_{m}\\},\\{w_{n}\\},\phi,\lambda,\\{t_{n}\\})]$
$\displaystyle\log
q^{*}_{2}\left(\\{w_{n}\\}_{n=1}^{N},\\{\hat{t}_{m},\hat{w}_{m}\\}_{m=1}^{M}\right)=$
(15) $\displaystyle\mathbb{E}_{q_{1}\left(\phi,\lambda\right)}[\log
P(\\{\hat{t}_{m},\hat{w}_{m}\\},\\{w_{n}\\},\phi,\lambda,\\{t_{n}\\})].$

Thus, to obtain the optimal distribution of one of the factors, one must
calculate expectations of the logarithm of the joint distribution over the
remaining factors in the approximation, resulting in an iterative algorithm.

In the following subsections, we explicitly express the functional form of the
optimal distributions, and obtain the corresponding expectations required for
updating the factors. The hyperparameters
($\\{\sigma_{g},a_{g},\alpha,\sigma_{s},a_{s}\\}$) are learned via gradients
update of the lower bound, we present details in the Supplementary Material.

###  4.3 Optimal $q_{1}$

We find that the optimal $q_{1}$ is factorized as

$\displaystyle
q_{1}\left(\phi,\lambda\right)=q_{1}\left(\lambda\right)q_{1}\left(\phi\right)$

The first factor is identified as a Gamma distribution

$\displaystyle q_{1}\left(\lambda\right)=Gamma\left(\alpha,\beta\right)$ (16)
$\displaystyle\alpha=\alpha_{0}+N+\int_{\mathcal{T}x\mathcal{W}}\Lambda_{q_{2}}\left(t,w\right)dtdw$
$\displaystyle\beta=\beta_{0}+T$

with known expectations.

Figure 1: (a) Comparison of the ground truth predictive intensity and the one
sampled from the VI and Gibbs inference. (b) Comparison of the ground truth
linear intensity $\phi\left(\cdot\right)$ and the ones learned by the VI and
Gibbs sampler. (c) Comparison of the Ground truth intensity bound and the one
learned by the inference, and the prior distribution. (d) The autocorrelation
of the intensity bound Gibbs samples. (e) The variational lower bound as a
function of the algorithm iteration. (f) Comparison of the test log–likelihood
of the Gibbs sampler and the VI.

The optimal distribution for the second factor is of the form

$\displaystyle q_{1}^{\star}$ $\displaystyle\propto
e^{-U\left(\phi\right)+\log p\left(\phi\right)}$ $\displaystyle U(\phi)$
$\displaystyle=\frac{1}{2}\int A(t)\phi^{2}(t)dt-\int b(t)\phi(t)dt$
$\displaystyle A(t)$
$\displaystyle=\sum_{n}\langle\omega_{n}\rangle_{q_{2}^{\star}}\delta(t-t_{n})+\langle\omega(t)\rangle_{q_{2}^{\star}}\Lambda_{q_{2}^{\star}}(t)$
$\displaystyle b(t)$
$\displaystyle=\sum_{n}\frac{1}{2}\delta(t-t_{n})-\frac{1}{2}\Lambda_{q_{2}^{\star}}\left(t\right).$

Generally, the integrals above cannot be evaluated analytically. Thus, we
resort to another variational approximation, where we approximate the
likelihood term, by a distribution that depends only on a finite set of
inducing point $\\{c\\}$,
$\tilde{q}\left(\phi_{c},\phi\right)=p\left(\phi|\phi_{c}\right)q\left(\phi_{c}\right)$
and the ELBO is

$\displaystyle\left\langle\log\frac{e^{-\log\langle
U(\phi)\rangle_{p(\phi|\phi_{c})}}p(\phi_{c})}{\tilde{q}(\phi_{c})}\right\rangle_{\tilde{q}}$

and we use the notation $\langle p\rangle_{q}=\mathbb{E}_{q}\left(p\right)$.
The optimal $\tilde{q}\left(\phi_{c}\right)$ is given by

$\displaystyle\tilde{q}^{\star}(\boldsymbol{\phi}_{c})\propto e^{-\log\langle
U(\boldsymbol{\phi})\rangle_{p(\boldsymbol{\phi}|\boldsymbol{\phi}_{c})}}p(\boldsymbol{\phi}_{c}).$

From here, using known results of conditional GPs and sparse variational GPs
(Csató et al., 2001; Titsias, 2009) we have

$\displaystyle\tilde{q}^{\star}(\boldsymbol{\phi}_{c})$
$\displaystyle=\mathcal{N}(\boldsymbol{\phi}_{c}|\boldsymbol{\mu}_{c},\Sigma_{c})$
(17) $\displaystyle\Sigma_{c}$
$\displaystyle=\left[\int\boldsymbol{\kappa}(t)^{\top}A(t)\boldsymbol{\kappa}(t)dt+K_{c}^{-1}\right]^{-1}$
$\displaystyle\boldsymbol{\mu}_{c}$ $\displaystyle=\Sigma_{c}\left(\int
b(t)\boldsymbol{\kappa}(t)dt\right)$

with $K_{c}$ the covariance kernel between the inducing points,
$\boldsymbol{\kappa}(t)=\boldsymbol{k}_{c}(t)^{\top}K_{c}^{-1}$ and
$\boldsymbol{k}_{c}(t)$ is the kernel between the inducing points and another
set of points (either the real data or the integration points). The mean and
the variance of the sparse approximated GP are

$\displaystyle\langle g(t)\rangle=$
$\displaystyle\boldsymbol{\kappa}(t)\boldsymbol{\mu}_{c}$ (18)
$\displaystyle\sigma^{2}(t)=$ $\displaystyle
K(t,t)-\boldsymbol{\kappa}(t)^{\top}\boldsymbol{k}_{c}(t)+\boldsymbol{\kappa}(t)^{\top}\Sigma_{c}\boldsymbol{\kappa}(t)$

###  4.4 Optimal $q_{2}$

Similarly to the previous section, we find that the optimal $q_{2}$ is
factorized as

$\displaystyle
q_{2}\left(\\{w_{n}\\}_{n=1}^{N},\Pi\right)=q_{2}\left(\\{w_{n}\\}_{n=1}^{N}\right)q_{2}\left(\\{\hat{t}_{m},\hat{w}_{m}\\}\right)$

Given Equation 13, we define the first factor as

$\displaystyle
q_{2}^{\star}(w_{n})\propto\exp\left(-\frac{\langle\phi_{n}^{2}\rangle_{q_{1}^{\star}}}{2}w_{n}\right)PG(w_{n}|1,0),$

which corresponds to a tilted PG distribution

$\displaystyle
q_{2}^{\star}(w_{n})=PG\left(w_{n}|1,\sqrt{\langle\phi_{n}^{2}\rangle_{q_{1}^{\star}}}\right).$
(19)

with known expectations (Polson et al., 2013).

The second factor takes the form

$\displaystyle q_{2}^{\star}(\\{\hat{t}_{m},\hat{w}_{m}\\}_{m=1}^{M})$
$\displaystyle\propto\prod_{m=1}^{M}\exp\left(-\frac{\langle\phi_{m}\rangle_{q_{1}^{\star}}}{2}-\frac{\langle\phi_{m}^{2}\rangle_{q_{1}^{\star}}}{2}w_{m}\right)\cdot\exp\left(\left\langle\ln\lambda^{\star}\right\rangle_{q_{1}^{\star}}\right).$

It can be shown that this distribution corresponds to a Poisson process with
intensity function

$\displaystyle\Lambda_{q_{2}}\left(\hat{t},\hat{w}\right)$ (20)
$\displaystyle=\exp\left(\langle\ln\lambda\rangle_{q_{1}^{\star}}\right)\frac{\exp\left(-\frac{\langle\phi\rangle_{q_{1}^{\star}}}{2}\right)}{2\cosh\left(\langle\phi^{2}\rangle_{q_{1}^{\star}}\right)}PG\left(w_{m}|1,\sqrt{\langle\phi^{2}\rangle_{q_{1}^{\star}}}\right)$

where to simplify the notation we write $\phi$ instead of
$\phi\left(\hat{t}\right)$.

We summarize the algorithm and discuss the hyper parameters learning in the
Supplementary Material.

##  5 Experiments

All the algorithm and experiments for this work are implemented in Python and
are available online. For further implementation details please see the
Supplementary Material.

Figure 2: (a) raster plot of a neuron from a monkey cortex. (b) Data generated
from the learned model. (c) Results of the Kolmogorov-Smirnov test. The NHGSE
generates data that resembles the real data, and passes the goodness of fit
test.

###  5.1 Synthetic Data

To assess the performance of the inference algorithms presented in Section 4,
we learn the parameters of data generated by the model, and compare the
learned parameters to the ground truth. To generate data we start by sampling
the memory GP and the background GP, based on Equations 6 and 7. We generate
events from the model using Poisson thinning (Lewis and Shedler, 1979). First
we sample the number of candidates $J\sim Poisson\left(\lambda T\right)$, and
sample candidate events $\\{t_{1},\dots t_{J}\\}$ uniformly. Next we
chronologically iterate through the candidates and accept them with
probability $\frac{\Lambda\left(t_{j}|\\{t_{1},\dots
t_{j-1}\\}\right)}{\lambda}$.

The results for the synthetic data are included in Figure 1. The time window
used was one second, and the dataset includes $18$ events. In panel (a), the
comparison between the ground truth predictive intensity and the ones inferred
by the learning algorithms demonstrates the accuracy of the inference methods.
We compare the ground truth to the mean of the Gibbs samples, and the mean of
the approximating distribution of the VI. Due to space limitation we do not
include the standard deviation of the inferred predictive intensity, and we
refer the reader to the Supplementary Material for further results.

Panel (b) compares the ground truth value of the linear intensity and the one
inferred by the two learning algorithms. Unlike for the predictive density,
the Gibbs sampler is more accurate than the VI. Panel (c) present the
inference results for the intensity bound. As expected the approximated
distribution by the VI is much more narrow than the distribution of the Gibbs
samples.

Panels (d) and (e) show the autocorrelation of the intensity bound, and the
ELBO through the Gibbs samples and VI iterations. In this example the
convergence of the ELBO is very fast, the autocorrelation of the Gibbs sample
vanish only after a few thousands iterations.

We use the test log–likelihood per data point, averaged over ten datasets, to
quantify the performance of the two inference algorithms. The Gibbs sampler
and the VI achieve very similar results. Thus, in the next section we present
the results only from the VI.

###  5.2 Real Data

####  5.2.1 Crime Report Data

Our model assumes both inhibitory and excitatory self effects, but it should
also be able to capture phenomena where only one of the two types of effects
exist. To test this, we fit our model to crime report data, where it is
assumed that past events have excitatory effect on future events (Mohler et
al., 2011). We use the same two datasets described in Zhou et al. (2020), and
follow their data processing procedure. Each dataset contains one type of
crime and so we use the univariate version of the model. The work of Zhou et
al. (2020) includes several inference methods and we compare our results to
the results of their reported mean–field variational inference approach, as it
is the closest to our inference procedure.

Table 1 compares the test log–likelihood of our NH-GPS model to the one
reported in Zhou et al. (2020). We perform the experiment five times and
report the mean and variance of the test log–likelihood. As expected, our
model performs similarly to the non–parametric Hawkes process presented by
Zhou et al. (2020).

Table 1: Crime Report Data Test Log–Likelihood. Dataset | Zhou et al. (2020) | NH–GPS
---|---|---
Vancouver | $453.11\pm 8.94$ | $453.8\pm 12.2$
NYPD | $-200.7\pm 3.32$ | $-202.8\pm 7.54$

####  5.2.2 Neuronal Activity Data

One of the motivating real world phenomena behind our work is the spiking
activity of neurons, where it is known that the process has both
self–excitatory and self–inhibitory effects. As an example for our model’s
ability to capture neuronal activity we use the datasets that were first
presented in Gerhard et al. (2017) (Figure 2.c and 2.b). These data were
further analyzed in Apostolopoulou et al. (2019) (Figure 5) where the Mutually
Regressive Point Process (MR-PP) is introduced. One dataset includes ten
recordings from a single neuron in a monkey cortex, with the duration of one
second each, and the other includes ten recordings from single neuron in a
human cortex for a duration of ten seconds each.

In Figure 2 we assess the ability of the model to capture the data for the
recordings from the monkey cortex. The results for the human cortex are
included in the Supplementary Material. Panel a and b present the raster plot
of the real data and the raster plot generated from the fitted model
respectively. Similarly to the real data, the generated data displays both
excitation, in the form of clustered events, and inhibition.

To quantify how suitable the model is to the data, we apply the random time
change theorem (Daley and Vere-Jones, 2003) to the inferred intensity and the
experimental data. The theorem states that realizations from a general point
process can be transformed to realizations from a homogeneous Poisson process
with unit rate. Similarly to the work of Apostolopoulou et al. (2019), we
further transform the exponential realizations to those from a uniform
distribution, following Brown et al. (2002). We then use the Kolmogorov-
Smirnov test to compare the quantiles of the distribution of the transformed
realizations to the quantiles of the uniform distribution. The results of this
test are displayed in the in Figure 2 c. The comparison relies between the
$95\%$ confidence bounds, which are indicated by the dashed lines. The model
passes the goodness of fit test (p value $>0.05$), and the p value is higher
than the one achieved by the MR-PP. We compare the reported p–value achieved
by the MR–PP model for the two datasets in Table 2. For both datasets our
model achieves higher p–value than the MR–PP.

Table 2: P–Value of the KS Test on Neuronal Activity Data. Dataset | MR–PP | NH–GPS
---|---|---
Monkey Cortex | 0.103 | 0.23
Human Cortex | 0.096 | 0.175

####  5.2.3 Retweets Data

Figure 3: Comparison of the theoretical quantiles and experimental quantiles
for each type of event in a retweets cascade dataset.

To demonstrate the necessity of the multivariate variation of the model we use
a subset of a retweets dataset. We used the processed data that was published
in Zuo et al. (2020). Each time series in the dataset contains the cascade of
retweets of a specific tweet and the event time is relative to the time of the
initial tweet. Each tweet is tagged according to the amount of followers the
user has, and we use the same grouping in to three groups as in Zuo et al.
(2020), which results in three event types.

It can be shown that the change time theorem can be applied separately to
every dimension of the multivariate intensity function (Daley and Vere-Jones,
2003), which in our case corresponds to each of the event types. As our model
is best fit to use cases were the data are scarce, we fit a model to each
retweets cascade separately. Figure 3 presents the QQ-plot for each type of
event for one retweet cascade dataset. All types pass the goodness–of–fit
test.

As a baseline we fit the univariate version of the model to the same data,
without differentiating between different event types. Table 3 presents the
comparison between the fit of the univariate and multivariate versions of the
model, averaged over five retweets cascades. It includes the p–values as
calculated by the Kolmogorov-Smirnov test, and the test log–likelihood per
data point. To calculate the test log–likelihood we trained the model on an
hour of retweets and tested it on the following half hour. Both the univariate
and the multivariate versions of the model pass the goodness–of–fit test, but
the multivariate model achieves a higher p–value, and a higher test
log–likelihood.

Table 3: Retweets Data P–Value and Test Log–Likelihood. Model | P–Value | Test Log–Likelihood
---|---|---
Univariate | $0.29\pm 0.15$ | $-1.99\pm 0.36$
Multivariate | $0.58\pm 0.25$ | $-0.8\pm 0.6$

##  6 Conclusion

In this work we presented the nonlinear Hawkes model with Gaussian process
self–effects (NH-GPS). We motivated the development of the new model with the
need for a flexible model that can capture both exciting and inhibiting
interactions between events, while maintaining the ability to learn also when
data are scarce.

Due to the structure of the model, we dispense with the branching structure
that is commonly used for Bayesian inference in Hawkes processes. We propose
an efficient mean–field variational inference algorithm which relies on a data
augmenting scheme. We show that the results of the variational inference are
comparable with those of a Gibbs sampler.

We demonstrate the performance of our model in three different real world
applications. Due to the flexibility of our model, it achieves good results on
data where events have only excitatory effects and on data where events have
both excitatory and inhibitory effects.

Last, we present the multivariate variation of the model. Here, we use it to
fit data where events are of different types. This direction can be extended
further. The multivariate version of the model can be used to describe events
in a network, and even learn the connections in the network. Another future
research direction is spatio–temporal processes. The model can be directly
extended to capture events in time and space, as the memory kernel is reduced
to a kernel function of a GP which can be applied also to multi–dimensional
data.

##### Acknowledgments

The authors would like to thank Dr. Christian Donner for topical discussions
which contributed to the research presented in this paper. The research of NMS
and MO have been partially funded by the Deutsche Forschungsgemeinschaft
(DFG)- Project-ID 318763901 - SFB1294. The research of CO was funded by the
BIFOLD-Berlin Institute for the Foundations of Learning and Data (ref.
01IS18025A and ref 01IS18037A).

## References

  * Ogata [1988] Yosihiko Ogata.  Statistical models for earthquake occurrences and residual analysis for point processes.  _Journal of the American Statistical association_ , 83(401):9–27, 1988.
  * Zhao et al. [2015] Qingyuan Zhao, Murat A Erdogdu, Hera Y He, Anand Rajaraman, and Jure Leskovec.  Seismic: A self-exciting point process model for predicting tweet popularity.  In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 1513–1522. ACM, 2015.
  * Dayan and Abbott [2001] Peter Dayan and Laurence F Abbott.  _Theoretical neuroscience: computational and mathematical modeling of neural systems_.  Computational Neuroscience Series, 2001.
  * Cox [1955] David R Cox.  Some statistical methods connected with series of events.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 17(2):129–157, 1955.
  * Hawkes and Oakes [1974] Alan G Hawkes and David Oakes.  A cluster process representation of a self-exciting process.  _Journal of Applied Probability_ , 11(3):493–503, 1974.
  * Maffei et al. [2004] Arianna Maffei, Sacha B Nelson, and Gina G Turrigiano.  Selective reconfiguration of layer 4 visual cortical circuitry by visual deprivation.  _Nature neuroscience_ , 7(12):1353–1359, 2004\.
  * Smith and Jahr [2002] T Caitlin Smith and Craig E Jahr.  Self-inhibition of olfactory bulb neurons.  _Nature neuroscience_ , 5(8):760–766, 2002.
  * Brémaud and Massoulié [1996] Pierre Brémaud and Laurent Massoulié.  Stability of nonlinear hawkes processes.  _The Annals of Probability_ , pages 1563–1588, 1996.
  * Zhu [2013] Lingjiong Zhu.  Central limit theorem for nonlinear hawkes processes.  _Journal of Applied Probability_ , 50(3):760–771, 2013.
  * Truccolo [2016] Wilson Truccolo.  From point process observations to collective neural dynamics: Nonlinear hawkes process glms, low-dimensional dynamics and coarse graining.  _Journal of Physiology-Paris_ , 110(4):336–347, 2016.
  * Jia and Benson [2019] Junteng Jia and Austin R Benson.  Neural jump stochastic differential equations.  _arXiv preprint arXiv:1905.10403_ , 2019.
  * Xiao et al. [2017] Shuai Xiao, Mehrdad Farajtabar, Xiaojing Ye, Junchi Yan, Le Song, and Hongyuan Zha.  Wasserstein learning of deep generative point process models.  In _Advances in Neural Information Processing Systems_ , pages 3247–3257, 2017.
  * Hensman et al. [2015] James Hensman, Alexander G de G Matthews, Maurizio Filippone, and Zoubin Ghahramani.  Mcmc for variationally sparse gaussian processes.  _arXiv preprint arXiv:1506.04000_ , 2015.
  * Lloyd et al. [2015] Chris Lloyd, Tom Gunter, Michael Osborne, and Stephen Roberts.  Variational inference for gaussian process modulated poisson processes.  In _International Conference on Machine Learning_ , pages 1814–1822. PMLR, 2015.
  * Adams et al. [2009] Ryan Prescott Adams, Iain Murray, and David JC MacKay.  Tractable nonparametric bayesian inference in poisson processes with gaussian process intensities.  In _Proceedings of the 26th Annual International Conference on Machine Learning_ , pages 9–16. ACM, 2009.
  * Donner and Opper [2018a] Christian Donner and Manfred Opper.  Efficient bayesian inference of sigmoidal gaussian cox processes.  _The Journal of Machine Learning Research_ , 19(1):2710–2743, 2018a.
  * Rasmussen [2013] Jakob Gulddahl Rasmussen.  Bayesian inference for hawkes processes.  _Methodology and Computing in Applied Probability_ , 15(3):623–642, 2013.
  * Du et al. [2016] Nan Du, Hanjun Dai, Rakshit Trivedi, Utkarsh Upadhyay, Manuel Gomez-Rodriguez, and Le Song.  Recurrent marked temporal point processes: Embedding event history to vector.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 1555–1564. ACM, 2016.
  * Mei and Eisner [2017] Hongyuan Mei and Jason M Eisner.  The neural hawkes process: A neurally self-modulating multivariate point process.  In _Advances in Neural Information Processing Systems_ , pages 6754–6764, 2017.
  * Jing and Smola [2017] How Jing and Alexander J Smola.  Neural survival recommender.  In _Proceedings of the Tenth ACM International Conference on Web Search and Data Mining_ , pages 515–524. ACM, 2017.
  * Zuo et al. [2020] Simiao Zuo, Haoming Jiang, Zichong Li, Tuo Zhao, and Hongyuan Zha.  Transformer hawkes process.  In _International Conference on Machine Learning_ , pages 11692–11702. PMLR, 2020.
  * Zhang et al. [2018] Rui Zhang, Christian Walder, Marian-Andrei Rizoiu, and Lexing Xie.  Efficient non-parametric bayesian hawkes processes.  _arXiv preprint arXiv:1810.03730_ , 2018.
  * Zhou et al. [2019] Feng Zhou, Zhidong Li, Xuhui Fan, Yang Wang, Arcot Sowmya, and Fang Chen.  Efficient em-variational inference for hawkes process.  _arXiv preprint arXiv:1905.12251_ , 2019.
  * Zhang et al. [2020] Rui Zhang, Christian Walder, and Marian-Andrei Rizoiu.  Variational inference for sparse gaussian process modulated hawkes process.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 6803–6810, 2020.
  * Zhou et al. [2020] Feng Zhou, Zhidong Li, Xuhui Fan, Yang Wang, Arcot Sowmya, and Fang Chen.  Efficient inference for nonparametric hawkes processes using auxiliary latent variables.  _Journal of Machine Learning Research_ , 21(241):1–31, 2020.
  * Wang et al. [2016] Yichen Wang, Bo Xie, Nan Du, and Le Song.  Isotonic hawkes processes.  In _International conference on machine learning_ , pages 2226–2234, 2016.
  * Apostolopoulou et al. [2019] Ifigeneia Apostolopoulou, Scott Linderman, Kyle Miller, and Artur Dubrawski.  Mutually regressive point processes.  In _Advances in Neural Information Processing Systems_ , pages 5115–5126, 2019.
  * Kingman [1993] John Frank Charles Kingman.  _Poisson processes_.  Wiley Online Library, 1993.
  * Daley and Vere-Jones [2007] Daryl J Daley and David Vere-Jones.  _An introduction to the theory of point processes: volume II: general theory and structure_.  Springer Science & Business Media, 2007.
  * Daley and Vere-Jones [2003] Daryl J Daley and David Vere-Jones.  _An introduction to the theory of point processes, volume 1: Elementary theory and methods_.  Verlag New York Berlin Heidelberg: Springer, 2003.
  * Polson et al. [2013] Nicholas G Polson, James G Scott, and Jesse Windle.  Bayesian inference for logistic models using pólya–gamma latent variables.  _Journal of the American statistical Association_ , 108(504):1339–1349, 2013.
  * Donner and Opper [2018b] Christian Donner and Manfred Opper.  Efficient bayesian inference for a gaussian process density model.  In _Conference on Uncertainty in Artificial Intelligence_. PMLR, 2018b.
  * Jordan et al. [1999] Michael I Jordan, Zoubin Ghahramani, Tommi S Jaakkola, and Lawrence K Saul.  An introduction to variational methods for graphical models.  _Machine learning_ , 37(2):183–233, 1999.
  * Bishop [2006] Christopher Bishop.  Pattern recognition and machine learning.  _Pattern Recognition and Machine Learning_ , 2006.
  * Csató et al. [2001] Lehel Csató, Manfred Opper, and Ole Winther.  Tap gibbs free energy, belief propagation and sparsity.  In _NIPS_ , pages 657–663, 2001.
  * Titsias [2009] Michalis Titsias.  Variational learning of inducing variables in sparse gaussian processes.  In _Artificial Intelligence and Statistics_ , pages 567–574, 2009\.
  * Lewis and Shedler [1979] PA W Lewis and Gerald S Shedler.  Simulation of nonhomogeneous poisson processes by thinning.  _Naval research logistics quarterly_ , 26(3):403–413, 1979.
  * Mohler et al. [2011] George O Mohler, Martin B Short, P Jeffrey Brantingham, Frederic Paik Schoenberg, and George E Tita.  Self-exciting point process modeling of crime.  _Journal of the American Statistical Association_ , 106(493):100–108, 2011.
  * Gerhard et al. [2017] Felipe Gerhard, Moritz Deger, and Wilson Truccolo.  On the stability and dynamics of stochastic spiking neuron models: Nonlinear hawkes process and point process glms.  _PLoS computational biology_ , 13(2):e1005390, 2017.
  * Brown et al. [2002] Emery N Brown, Riccardo Barbieri, Valérie Ventura, Robert E Kass, and Loren M Frank.  The time-rescaling theorem and its application to neural spike train data analysis.  _Neural computation_ , 14(2):325–346, 2002.
  * Ogata [1981] Yosihiko Ogata.  On lewis’ simulation method for point processes.  _IEEE transactions on information theory_ , 27(1):23–31, 1981.
  * Rasmussen and Williams [2006] Carl Edward Rasmussen and Christopher KI Williams.  _Gaussian process for machine learning_.  MIT press, 2006.
  * Gilks and Wild [1992] Walter R Gilks and Pascal Wild.  Adaptive Rejection Sampling for Gibbs Sampling.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 41(2):337–348, 1992.
  * Martino et al. [2015] Luca Martino, Hanxue Yang, David Luengo, Juho Kanniainen, and Jukka Corander.  A fast universal self-tuned sampler within Gibbs sampling.  _Digital Signal Processing_ , 47:68–83, 2015.
  * Hastings [1970] W Keith Hastings.  Monte carlo sampling methods using markov chains and their applications.  1970\.
  * Duane et al. [1987] Simon Duane, Anthony D Kennedy, Brian J Pendleton, and Duncan Roweth.  Hybrid monte carlo.  _Physics letters B_ , 195(2):216–222, 1987.
  * Bradbury et al. [2018] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang.  JAX: composable transformations of Python+NumPy programs, 2018\.  URL http://github.com/google/jax.
  * Linderman [2017] Scott Linderman.  Pypólyagamma.  https://github.com/slinderman/pypolyagamma, 2017.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.

Supplementary Materials

##  1 Gibbs Sampler

We use a blocked Gibbs sampler, which groups two or more variables and sample
them at once. Thus, we need to identify the conditional posterior distribution
of all the relevant groups.

###  1.1 Conditional Distribution of the Upper Intensity Bound

The conditional distribution of the upper intensity bound is

$\displaystyle
p\left(\lambda|\\{\hat{t}_{m},\hat{w}_{m}\\},\\{w_{n}\\},\phi,\\{t_{n}\\}\right)\propto
e^{-\lambda T}\lambda^{N+M}p\left(\lambda\right)$

which we identify as a Gamma distribution

$\displaystyle p$
$\displaystyle\left(\lambda|\\{\hat{t}_{m},\hat{w}_{m}\\},\\{w_{n}\\},\phi,\\{t_{n}\\}\right)\propto
Gamma\left(\alpha,\beta\right)$ (21) $\displaystyle\alpha$
$\displaystyle=\alpha_{0}+N+M$ $\displaystyle\beta$
$\displaystyle=\beta_{0}+T.$

###  1.2 Conditional Distribution of the Linear Intensity Function

The conditional distribution of the linear intensity function in the observed
and augmenting events is

$\displaystyle
p\left(\phi_{N+M}|\\{\hat{t}_{m},\hat{w}_{m}\\},\\{w_{n}\\},\phi,\\{t_{n}\\},\lambda\right)$
$\displaystyle\propto\exp\left(\sum_{n=1}^{N}f\left(w_{n},\phi_{n}\right)+\sum_{m=1}^{M}f\left(\hat{w}_{m},-\phi_{m}\right)\right)p\left(\phi_{N+M}\right),$

where we use the shortened notation $\phi_{n}$ instead of
$\phi\left(t_{n}\right)$ and $\phi_{N+M}$ instead of
$\\{\\{\phi\left(t_{n}\right)\\}_{n=1}^{N},\\{\phi\left(\hat{t}_{m}\right)\\}_{m=1}^{M}\\}\\}$.
Given Equation 9 in the main text, the likelihood term in the posterior is a
GP with mean

$\displaystyle\mu=\left(\frac{1}{2w_{1}},...,\frac{1}{2w_{N}},-\frac{1}{2\hat{w}_{1}},-\frac{1}{2\hat{w}_{M}}\right)^{\top}$

and diagonal covariance matrix

$\displaystyle\Sigma^{-1}=Diag\left(w_{1},...,w_{N},\hat{w}_{1},\hat{w}_{M}\right).$

Given the GP prior over $\phi$ the conditional posterior is also a GP

$\displaystyle\phi_{N+M}\sim GP\left(\mu_{M+N},\Sigma_{M+N}\right)$ (22)
$\displaystyle\mu_{M+N}=\Sigma_{N+M}\Sigma^{-1}\mu$
$\displaystyle\Sigma_{N+M}^{-1}=\Sigma^{-1}+\tilde{K}^{-1}.$

###  1.3 Conditional Distribution of the PG variables

The conditional posterior distribution of the augmenting PG variables is

$\displaystyle p\left(\\{w_{n}\\},\\{\hat{w}_{m}\\}\right)$
$\displaystyle\propto\prod_{n=1}^{N}\exp\left(-\frac{\phi_{n}^{2}}{2}w_{n}\right)$
$\displaystyle\times
PG\left(w_{n};1,0\right)\prod_{m=1}^{M}\exp\left(-\frac{\phi_{m}^{2}}{2}\hat{w}_{m}\right)PG\left(w_{m};1,0\right).$

Using the definition of the tilted PG distribution [Polson et al., 2013] the
posterior distribution for these parameters is

$\displaystyle w_{n}\propto PG\left(1,\phi_{n}\right)$
$\displaystyle\hat{w}_{m}\propto PG\left(1,\phi_{m}\right).$

###  1.4 Conditional Distribution of the augmenting events

The conditional posterior of the augmenting events is proportional to

$\displaystyle
p\left(\\{\hat{t}_{m}\\}_{m=1}^{M}|\\{t_{n}\\}_{n=1}^{N},\\{\hat{w}_{m}\\}_{m=1}^{M},\\{w_{n}\\}_{n=1}^{N},\phi,\lambda\right)$
$\displaystyle\propto\prod_{m=1}^{M}\lambda
e^{f\left(\hat{w}_{m},-\phi_{m}\right)}PG\left(\hat{w}_{m};1,0\right).$

We recognize the conditional posterior with the unnormalized density of a
marked inhomogeneous Poisson process with intensity

$\displaystyle\Pi_{0}\left(\hat{w},\hat{t}\right)=\lambda
e^{f\left(\hat{w}_{m},-\phi_{m}\right)}PG\left(\hat{w}_{m};1,0\right).$ (23)

The underlying density of the inhomogeneous Poisson process in the
realizations space $\mathcal{T}$ is given by

$\displaystyle\int_{\mathcal{W}}\Pi_{0}\left(\hat{w},\hat{t}\right)d\hat{w}=\lambda\sigma\left(-\phi_{m}\right).$

To sample the augmenting realizations we use the thinning algorithm [Ogata,
1981]. First, we sample the expected number of events

$\displaystyle J\sim Poisson\left(\lambda T\right).$

Next, $J$ candidates are sampled uniformly over the realizations space
$\mathcal{T}$. To evaluate the intensity at the candidate points we need to
evaluate the linear intensity $\phi$ in these points, given its values in the
real events and the previously sampled augmenting events. This can be done
using results from GP regression [Rasmussen and Williams, 2006]

$\displaystyle\phi_{J|N+M}\sim GP\left(\mu,\tilde{K}\right)$
$\displaystyle\mu=\tilde{K}_{J,N+M}\tilde{K}_{N+M}^{-1}\phi_{N+M}$
$\displaystyle\tilde{K}=\tilde{K}_{J}-\tilde{K}_{J,N+M}\tilde{K}_{N+M}^{-1}\tilde{K}_{J,N+M}^{\top}.$

Once we have $\phi_{J}$ we perform thinning - for each candidate $\hat{t}_{j}$
we generate a random number $r_{j}$ between $0$ and $1$. If
$r_{j}<\sigma\left(-\phi_{j}\right)$ we accept the candidate $\hat{t}_{j}$ and
otherwise we discard it.

Input: Observed events $\\{t_{n}\\}_{n=1}^{N}$, hyper parameters
$\\{\sigma_{g},\alpha,a_{s},\sigma_{s},\alpha_{0},\beta_{0}\\}$

Output: R samples of $\\{\hat{t}_{m}\\}_{m=1}^{M_{r}}$, $\lambda$,
$\phi_{N+M_{r}}$

1 Initialize - $M,\phi_{N+M},\\{\hat{t}_{m}\\}_{m=1}^{M_{r}}$ randomly. for
_$r\leftarrow 0$ to $R$_ do

2 Sample $w\sim PG\left(1,\phi_{N+M}\right)$

3 Sample $\lambda\sim\text{Gamma}\left(\alpha,\beta\right)$ as in Equation 21

4 Sample $\phi_{N+M}\sim GP\left(\mu_{N+M},\Sigma_{N+M}\right)$ as in Equation
22

5 Sample $\\{\hat{t}\\}_{m=1}^{M_{r}}$ as in Section 1.4

6 end for

Algorithm 1 NH-GPS Gibbs Sampler

###  1.5 Hyper parameters Learning

The augmented model is not conditionally conjugated with respect to the kernel
hyperparameters. This is usually solved by using MCMC within Gibbs sampler
approach [Gilks and Wild, 1992, Martino et al., 2015]. This method applies
rejection sampling, such as Metropolis–Hastings (MH) [Hastings, 1970],
Hamiltonian Monte Carlo (HMC) [Duane et al., 1987], to sample the
hyperparameters, and relies heavily on some design choice. A wrong choice of
the proposal distribution (for MH) or the mass matrix (for HMC) may result in
a very slow convergence, or prevent the sampler from converging at all.

We choose the less traditional approach of taking a gradient step within the
Gibbs sampler. This is implemented in the following way – after sampling all
the model parameters from the conditional posterior distributions described
above, we derive the negative model log–posterior with respect to the
hyperparameters and take a step in the direction of the negative gradient.

This approach can be developed further in the spirit of Stochastic Gradient
Descent (SGD). Meaning, rather than updating the hyper parameters after each
iteration of the Gibbs sampler, we perform several steps of sampling, take the
gradient of the averaged posterior and update the hyper parameters following
the averaged gradient.

We include bellow the derivatives with respect for the hyper parameters of the
model $\theta={\sigma_{g},\alpha,a_{g},\sigma_{s}}$.

To learn the hyperparameters we derive the posterior of the model which
appears in Equation 13 in the main text. First, we notice that all of the
hyper parameters appear in the prior over the linear intensity

$\displaystyle\log
p\left(\phi\right)\propto-\frac{1}{2}\log\det\left(\tilde{K}\right)-\frac{1}{2}\phi^{\top}\tilde{K}^{-1}\phi$

and all of the hyperparameters appear in the prior kernel.

We next derive an entry in the kernel with respect to the different hyper
parameters.

$\displaystyle\frac{\partial\tilde{K_{l,k}}}{\partial a_{s}}$
$\displaystyle=\exp\left(-\frac{\parallel
t_{l}-t_{k}\parallel^{2}}{\sigma_{s}^{2}}\right)$
$\displaystyle\frac{\partial\tilde{K_{l,k}}}{\partial\sigma_{s}}$
$\displaystyle=a_{s}\exp\left(-\frac{\parallel
t_{l}-t_{k}\parallel^{2}}{\sigma_{s}^{2}}\right)\frac{\parallel
t_{l}-t_{k}\parallel^{2}}{\sigma_{s}^{3}}$
$\displaystyle\frac{\partial\tilde{K_{l,k}}}{\partial\alpha}$
$\displaystyle=\sum_{t_{i}<t_{l}}\sum_{t_{j}<t_{k}}K_{t_{i}-t_{l},t_{j}-t_{k}}^{g}$
$\displaystyle\times\exp\left(-\alpha\left(t_{i}-t_{l}+t_{j}-t_{k}\right)\right)\left(t_{l}-t_{i}+t_{k}-t_{j}\right)$
$\displaystyle\frac{\partial\tilde{K_{l,k}}}{\partial\sigma_{g}}$
$\displaystyle=\sum_{t_{i}<t_{l}}\sum_{t_{j}<t_{k}}K_{t_{i}-t_{l},t_{j}-t_{k}}^{g}$
$\displaystyle\times\exp\left(-\alpha\left(t_{i}-t_{l}+t_{j}-t_{k}\right)\right)\frac{\parallel\left(t_{l}-t_{i}\right)-\left(t_{k}-t_{j}\right)\parallel^{2}}{\sigma_{g}^{3}}.$

We can plug these results to the chain rule and we get

$\displaystyle\nabla\log
p\left(\phi\right)=-\frac{1}{2}\text{trace}\left(\tilde{K}^{-1}\nabla\tilde{K}\right)+\frac{1}{2}\phi^{\top}\tilde{K}^{-1}\nabla\tilde{K}\tilde{K}^{-1}\phi.$

##  2 Variational Inference

Input: Observed Events $\\{t_{n}\\}_{n=1}^{N}$, hyper parameters
$\\{\sigma_{g},\alpha,a_{s},\sigma_{s},\alpha_{0},\beta_{0}\\}$

Output: Mean and variance of $\phi$

1 Initialize. while  _$\mathcal{L}$ not converged_ do

2 Estimate the mean and covariance of the GP over the inducing points as in
Equation 17

3 Estimate the mean and variance of the GP over the observations as in
Equation 18 Estimate the mean intensity bound as in Equation 16

4 Estimate the mean and variance of the augmenting PG variables as in Equation
19

5 Estimate the intensity function over the augmenting events as in Equation 20

6 end while

Algorithm 2 NH-GPS Variational Inference

###  2.1 Hyper Parameters Learning

In order to learn the hyperparameters of the model we perform one step of
gradient descent with respect to the ELBO at each iteration. The derivatives
of the ELBO are give in Donner and Opper [2018a] Appendix F.

##  3 Implementation Details

All the experiments were implemented in Python. To parallelize the computation
over the available computing resources we used the JAX package [Bradbury et
al., 2018]. In the Gibbs sampler - the sampling of the PG variables was done
using the PyPólyaGamma package [Linderman, 2017].

##### Gradient step within the Gibbs sampler

We update the hyper parameters of the model using SGD. As an optimizer we use
the ADAM algorithm [Kingma and Ba, 2014]. We found that updating the
hyperparameters every second iteration yields the best results in terms of
running time and convergence time. The ADAM optimizer was also used for the
gradient step in the Variational Inference algorithm. We used the recommended
default variables of the ADAM optimizer.

###  3.1 Hyper parameters values

The parameters values that were used to generate the synthetic dataset
presented in Figure 1 are presented in Table 4. The hyper parameters for the
Variational inference, for all the real world datasets, are presented in Table
5. The noise level refers to the noise added the diagonal of the relevant
covariance matrices.

Table 4: Model Parameters Values for Synthetic Data. Parameter | Value
---|---
Time limit | 1
$\lambda$ | 320
$\sigma_{g}$ | 0.07
$a_{g}$ | 1
$\alpha$ | 20
$\sigma_{s}$ | 0.5
$a_{s}$ | 1

##  4 Further Results

###  4.1 Synthetic Data

Figure 4 includes three examples of a generated ground truth intensity
function and the inference results of the Gibbs sampler and variational
inference. The upper panel includes the same data as Figure 1 in the main
text, but includes also the standard deviation of the inference results.

###  4.2 Neuronal Activity Data

Figure 5 presents the results for the dataset recorded from a single neuron in
a human cortex. Similarly to the results in the main text, the data generated
by the model resembles the recordings. Furthermore, the fitted model passes
the goodness–of–fit test and the comparison between the theoretical and
experimental quantiles are within the $95\%$ confidence bounds.

Table 5: Hyper Parameters Values for the VI Algorithm. Parameter |  Synthetic Data |  Vancouver Crime |  NYPD Data |  Neuronal Data |  Retweets Data
---|---|---|---|---|---
$\alpha_{0}$ | $2.5$ | $1.5$ | $1.3$ | $2$ | $1.5$
$\beta_{0}$ | $0.007$ | $0.1$ | $0.1$ | $0.03$ | $0.1$
inducing points | $100$ | $200$ | $300$ | $100$ | $200$
integration points | $1000$ | $5000$ | $3000$ | $1000$ | $3600$
Noise | $10^{-4}$ | $10^{-4}$ | $10^{-4}$ | $10^{-4}$ | $10^{-4}$
Figure 4: Three examples for the ground truth intensity function and the one
inferred by the Gibbs and VI algorithms. The dashes lines are the standard
deviation of the inference results. Figure 5: (a) raster plot of a neuron from
a monkey cortex. (b) Data generated from the learned model. (c) Results of the
Kolmogorov-Smirnov test. The NHGSE generates data that resembles the real
data, and passes the goodness of fit test.