# Hawking radiation in quantum Hall system with an expanding edge:
application of anomaly method

Riku Yoshimoto<EMAIL_ADDRESS>Yasusada Nambu
<EMAIL_ADDRESS>Department of Physics, Graduate School
of Science, Nagoya University, Chikusa, Nagoya 464-8602, Japan

###### Abstract

The relationship between gravitational anomalies and Hawking radiation of
black holes was revealed by Wilczek and Robinson. In this study, we apply
their method to an analogue de Sitter spacetime in the quantum Hall system
with an expanding edge. Because this system is chiral, there is no need to
impose the condition of ingoing modes near the horizon, which was assumed in
the original method. Moreover, this system is structured so that the de Sitter
space is sandwiched between two flat spaces, and although the effects of the
anomaly would not appear in an ordinal de Sitter spacetime, they manifest
themselves as boundary conditions between the de Sitter and the flat regions.
By performing calculations under these boundary conditions, we obtain the flux
of Hawking radiation in the outer flat region with the Gibbons-Hawking
temperature of the de Sitter horizon.

##  I Introduction

Hawking radiation has been studied extensively as a quantum property of black
holes. Because of its low temperature, observing it from actual black holes in
the universe is extremely challenging. However, in condensed matter systems,
it is possible to construct analog models with properties such as black hole
horizons and cosmological horizons. Among them, the analog de Sitter spacetime
with quantum Hall (QH) systems is discussed in [1, 2]. The usual experiments
on QH systems are performed in a static situation where the edge is not time-
varying, but by making it time-varying, the excitations moving along the
expanding edge can be a simulator of the chiral scalar field in a
(1+1)-dimensional curved spacetime. Also, in [1, 2], the expanding edge can
reproduce a two-dimensional de Sitter universe which has the spacetime
structure similar to that of a black hole formation via gravitational
collapse. It was shown that the entanglement between two spatial regions
decreases due to Hawking radiation from the de Sitter region, and when the
detection region is sufficiently large compared to the Hubble length of the de
Sitter region, the two regions become separable, and only classical
correlations remain as well as the inflationary universe.

In this paper, we focus on the chiral properties of edge excitations rather
than on quantum correlations in the QH system. This implies the presence of
chiral scalar fields in the analog de Sitter spacetime. In such a case, it is
known that there is a gravitational anomaly [3, 4, 5]. Hawking radiation has
been studied using various approaches beyond Hawking’s original calculation
[6]. These include approaches using Euclidean quantum gravity by Hawking and
Gibbons [7], calculations using the trace anomaly by Christensen and Fulling
[8, 9], and calculations based on quantum tunneling [10]. Robinson and Wilczek
provided a new calculation of Hawking radiation using gravitational anomaly
[11]. Classical ingoing modes do not affect the dynamics of the external
region near the black hole horizon. Thus, when they are ignored, the theory
becomes chiral and a gravitational anomaly appears. Since the original theory
is general covariant, this anomaly needs to be canceled by quantum effects,
leading to the interpretation of obtaining Hawking radiation. This method has
been applied not only to Schwarzschild black holes, but also to Kerr black
holes, Reissner-Nordström black holes, and others [12, 13, 14, 15].

However, it is known that this method provides incorrect results for the Unruh
effect and cosmological particle creation in the de Sitter spacetime [16, 17].
In this study, we evaluated the flux of Hawking radiation from the de Sitter
horizon in the analog de Sitter spacetime using the anomaly method. In the
usual de Sitter spacetime, covariant gravitational anomaly method fails
because the anomaly term becomes zero in de Sitter spacetime; because the
spacetime curvature is constant and provides no terms to cancel out the
anomaly. However, in the present QH system, the anomaly term does not become
zero at the boundary. Therefore, the boundary conditions between the Minkowski
and de Sitter regions are determined by the anomaly equation, and we can
obtain the flux with Gibbons-Hawking temperature.

The structure of the paper is as follows. In Section II, we review the
coordinate transformations of each region in the expanding edge of QH systems
introduced in [2, 1] and organize the coordinates to be used in later
sections. In Section III, we calculate the expectation value of the energy-
momentum tensor (EMT) in the analog de Sitter spacetime introduced in Section
II based on the anomalous conservation law and the boundary conditions
obtained from it. Section IV is devoted to summary and conclusions.

##  II Quantum Hall system with expanding edge

In this section, we review the coordinate setup for the expanding edge of the
QH system [2, 1]. We consider the edge of the QH system that has left-moving
edge excitation. Then we divide the entire edge into three regions as shown in
Fig. 1; Region I ($x>L/2$) and region III ($x<-L/2$) are static, while region
II ($-L/2<x<L/2$) expands with time.

Figure 1: Schematic picture of an expanding edge in QH systems. Regions I and
III are static and region II is the expanding region.

The analog metric of this system is written as

$\displaystyle\text{Region {{III}}:}~{}x_{\text{{III}}}\leq-L/2,\quad
ds^{2}=-dx_{\text{{III}}}^{+}\,dx_{\text{{III}}}^{-},$
$\displaystyle\text{Region {{II}}:}~{}|x|\leq L/2,\quad
ds^{2}=-e^{2\Theta(t)}dx^{+}dx^{-},$ $\displaystyle\text{Region I:}~{}L/2\leq
x_{\text{I}},\quad ds^{2}=-dx_{\text{I}}^{+}dx_{\text{I}}^{-},$

where $e^{\Theta(t)}$ is the scale factor in region II and
$x^{\pm}_{\text{I,{III}}}=t_{\text{I,{III}}}\pm
x_{\text{I,{III}}},~{}x^{\pm}=t\pm x$ are null coordinates in each region.

The conformal flat coordinates in region II become discontinuous at the
boundary $x=\pm L/2$, so it is not possible to smoothly extend them to regions
I and III. However, it is possible to introduce smooth coordinates
$x_{\rm{I}}$ ($x_{\text{{III}}}$) that cover both regions I and II (or II and
III), while keeping $x_{\rm{I}}=L/2$ ($x_{\text{{III}}}=-L/2$).

###  II.1 Region II and III

First, we perform the coordinate transformation between the regions II and III
and express the metric in the region II with $x_{\text{{III}}}^{\pm}$. The
matching condition at the boundary $x_{\text{{III}}}=-L/2$ is

$x_{\text{{III}}}^{+}\left[t-L/2\right]-x_{\text{{III}}}^{-}\left[t+L/2\right]=-L.$
(1)

By taking a derivative, we obtain

$\frac{dx_{\text{{III}}}^{+}}{dx^{+}}\left[t-L/2\right]=\frac{dx_{\text{{III}}}^{-}}{dx^{-}}\left[t+L/2\right].$
(2)

Since the conformal factor is unity at $x=x_{\text{{III}}}=-L/2$, we obtain

$\frac{dx_{\text{{III}}}^{+}}{dx^{+}}\left[t-L/2\right]=e^{\Theta(t)}.$ (3)

By shifting the argument of this function,

$\frac{dx_{\text{{III}}}^{+}}{dx^{+}}\left[x^{+}\right]=e^{\Theta(x^{+}+L/2)},$
(4)

and one finds

$\displaystyle x_{\text{{III}}}^{+}\left[x^{+}\right]$
$\displaystyle=\int^{x^{+}}_{0}dy\,e^{\Theta(y+L/2)}=\int^{x^{+}+L/2}_{L/2}dy\,e^{\Theta(y)}$
$\displaystyle=\Phi\left[x^{+}+L/2\right]-\Phi\left[L/2\right],$ (5)

where the function $\Phi\left[x\right]$ is introduced by

$\Phi\left[x\right]=\int^{x}_{0}dy\,e^{\Theta(y)}.$ (6)

The coordinate $x_{\text{{III}}}^{-}$ is determined by the matching condition:

$\displaystyle x_{\text{{III}}}^{-}\left[x^{-}\right]$
$\displaystyle=x_{\text{{III}}}^{+}\left[x^{-}-L\right]+L$
$\displaystyle=\int^{x^{-}-L}_{L/2}dy\,e^{\Theta(y+L/2)}+L$
$\displaystyle=\Phi\left[x^{-}-L/2\right]-\Phi\left[L/2\right]+L.$ (7)

From these equations, the metric in region II can be written as

$\displaystyle ds^{2}_{\text{{II}}}$
$\displaystyle=-e^{2\Theta(t)}dx^{+}dx^{-}$
$\displaystyle=-\exp\left[2\Theta(t)-\Theta(x^{+}+L/2)-\Theta(x^{-}-L/2)\right]dx_{\text{{III}}}^{+}dx_{\text{{III}}}^{-}.$
(8)

At the boundary, $x^{+}=t-L/2$ and $x^{-}=t+L/2$, so the conformal factor
becomes unity. Thus, the coordinates $x^{\pm}$ in region II can be extended to
region III.

###  II.2 Region I and II

Let us consider the matching conditions at $x=x_{\rm{I}}=L/2$ :

$x_{\rm{I}}^{+}\left[t+L/2\right]-x_{\rm{I}}^{-}\left[t-L/2\right]=L.$ (9)

By similar calculation as matching regions II and III, we obtain

$x_{\rm{I}}^{+}\left[x^{+}\right]=\int^{x^{+}}_{0}dy\,e^{\Theta(y-L/2)}=\Phi\left[x^{+}-L/2\right]-\Phi\left[-L/2\right],$
(10)

and

$\displaystyle x_{\rm{I}}^{-}\left[x^{-}\right]$
$\displaystyle=x_{\rm{I}}^{+}\left[x^{-}+L\right]-L$
$\displaystyle=\int^{x^{-}+L}_{0}dy\,e^{\Theta(y-L/2)}-L$
$\displaystyle=\Phi\left[x^{-}+L/2\right]-\Phi\left[-L/2\right]-L.$ (11)

The metric in region II becomes

$\displaystyle ds^{2}_{\text{{II}}}$
$\displaystyle=-e^{2\Theta(t)}dx^{+}dx^{-}$
$\displaystyle=-\exp\left[2\Theta(t)-\Theta(x^{+}-L/2)-\Theta(x^{-}+L/2)\right]dx_{\rm{I}}^{+}dx_{\rm{I}}^{-}.$
(12)

At the boundary, $x^{+}=t+L/2$ and $x^{-}=t-L/2$, so the conformal factor
becomes unity, and the coordinates $x^{\pm}$ can be extended into region I.

###  II.3 Region I and III

Using the function $\Phi\left[x\right]$ introduced above, we can express the
relationship between regions I and III. Coordinates $x_{\rm{I}}^{\pm}$ and
$x_{\text{{III}}}^{\pm}$ are written as

$\displaystyle x_{\text{{III}}}^{+}=$
$\displaystyle\Phi(x^{+}+L/2)-\Phi(L/2),$ (13) $\displaystyle x_{\rm{I}}^{+}=$
$\displaystyle\Phi(x^{+}-L/2)-\Phi(-L/2),$ (14) $\displaystyle
x_{\text{{III}}}^{-}=$ $\displaystyle\Phi(x^{-}-L/2)-\Phi(L/2)+L,$ (15)
$\displaystyle x_{\rm{I}}^{+}=$ $\displaystyle\Phi(x^{-}+L/2)-\Phi(-L/2)-L.$
(16)

After eliminating $x^{\pm}$, we obtain relations between $x_{\text{I}}^{\pm}$
and $x_{\text{{III}}}^{\pm}$ as

$\displaystyle
x_{\rm{I}}^{+}=\Phi\left[-L+\Phi^{-1}\left[x_{\text{{III}}}^{+}+\Phi\left[L/2\right]\right]\right]-\Phi\left[-L/2\right],$
(17) $\displaystyle
x_{\text{{III}}}^{-}=\Phi\left[-L+\Phi^{-1}\left[x_{\rm{I}}^{-}+L+\Phi\left[-L/2\right]\right]\right]-\Phi\left[L/2\right]+L.$
(18)

###  II.4 Example of analog spacetime : de Sitter case

For later usage in the analysis of the de Sitter spacetime, we calculate the
specific form of $\Phi\left[x\right]$ and $\Phi^{-1}\left[x\right]$ (see Fig.
2 for this setup). The conformal factor of the de Sitter spacetime is

$e^{\Theta(t)}=\frac{1}{\cos(Ht)}.$ (19)

Figure 2: Penrose diagram for the QH system with a de Sitter edge region II. A
future horizon and a past horizon exist (time-symmetric eternal case). The
vacuum condition is imposed on the past null infinity $\mathcal{I}^{-}$ in
region I.

Substituting this into the definition of $\Phi\left[x\right]$ (6), we obtain

$\displaystyle\Phi(x)=\frac{1}{2H}\ln\frac{1+\sin Hx}{1-\sin Hx},$ (20)
$\displaystyle\Phi^{-1}(x)=\frac{1}{H}\arcsin\tanh Hx.$ (21)

Then the metric in region II is written as

$\displaystyle
ds^{2}_{\text{{II}}}=-\frac{dx_{\rm{I}}^{+}dx_{\rm{I}}^{-}}{\cosh^{2}\left[H(x_{\rm{I}}-L/2)\right]},$
(22) $\displaystyle
ds^{2}_{\text{{II}}}=-\frac{dx_{\text{{III}}}^{+}\,dx_{\text{{III}}}^{-}}{\cosh^{2}\left[H(x_{\text{{III}}}+L/2)\right]}.$
(23)

In the following, we will refer to region II covered by $x_{\rm{I}}$ as A and
region covered by $x_{\text{{III}}}$ as B (Fig. 3).

Figure 3: Region A in the de Sitter region II is covered by
$x_{\text{I}}^{\pm}$ and region B in the de Sitter region II is covered by
$x_{\text{{III}}}^{\pm}$.

##  III Application of Anomaly method

This QH system reproduces a left-moving chiral field in de Sitter spacetime.
It is known that in such a situation, not only a trace anomaly but also a
gravitational anomaly exists. In $(1+1)$-dimensional spacetimes, they are
represented by the following expression which indicates a violation of the
conservation law for the expectation value of EMT [3, 4, 5]:

$\displaystyle T^{\alpha}_{\alpha}=\frac{R}{48\pi},$ (24)
$\displaystyle\nabla_{\mu}T^{\mu}_{\nu}=\frac{\epsilon_{\mu\nu}\partial^{\mu}R}{96\pi\sqrt{-g}},$
(25)

where $R$ denotes the spacetime scalar curvature and $\epsilon_{\mu\nu}$ is
the totally antisymmetric tensor in two dimensions. By solving the anomaly
equation (25), we can observe the emergence of Hawking radiation in region
III.

###  III.1 General solution of anomaly equation

First, we determine the general solution of the EMT in each region. The scalar
curvature of this system is given by

$R=2H^{2}\theta(-x_{\rm{I}}+L/2)\theta(x_{\text{{III}}}+L/2)$ (26)

where the unit step function defined by

$\theta(x)=\begin{cases}0&\quad\text{for}\quad x<0,\\\ 1&\quad\text{for}\quad
x\geq 0.\end{cases}$ (27)

The EMT is conserved in bulk regions beside boundaries. Thus in regions I and
III (Minkowski region), we solve

${T_{i}}_{\alpha}^{\alpha}=0,\quad\partial_{\mu}{T_{i}}^{\mu}_{\nu}=0,\quad(i=\text{I,
{III}}).$ (28)

We want to consider the general solution of Eq. (28). By using
$T^{+}_{+}=T^{-}_{-}=1/2(T^{t}_{t}+T^{x}_{x})=1/2\,T^{\alpha}_{\alpha}$, Eq.
(28) becomes

$\partial_{-}{T_{i}}^{-}_{+}=0,\quad\partial_{+}{T_{i}}^{+}_{-}=0.$ (29)

Thus, we obtain

${T_{i}}^{-}_{+}=-2{T_{i}}_{++}=F_{i}\left[x^{+}_{i}\right],\quad{T_{i}}^{+}_{-}=-2{T_{i}}_{--}=G_{i}\left[x^{-}_{i}\right],$
(30)

where $F_{i}\left[x^{+}_{i}\right]$ is an arbitrary function of $x^{+}_{i}$
and $G_{i}\left[x^{-}_{i}\right]$ is an arbitrary function of $x^{-}_{i}$.

In regions A and B (de Sitter region), the EMT obeys

${T_{j}}_{\alpha}^{\alpha}=\frac{R}{48\pi},\quad\nabla_{\mu}{T_{j}}^{\mu}_{\nu}=\frac{\epsilon_{\mu\nu}\partial^{\mu}R}{96\pi\sqrt{-g_{j}}},\quad(j=A,B).$
(31)

The non-zero Christoffel symbols are
${\Gamma_{j}}^{+}_{++}=h_{j}\partial_{+}(1/h_{j})$ and
${\Gamma_{j}}^{-}_{--}=h_{j}\partial_{-}(1/h_{j})$ with
$h_{A}:=\cosh^{2}\left[H(x_{\rm{I}}-L/2)\right]$ and
$h_{B}:=\cosh^{2}\left[H(x_{\text{{III}}}+L/2)\right]$. Therefore, we obtain
the following explicit form of the anomaly equation:

$\begin{split}\partial_{-}{T_{A}}^{-}_{+}+{\Gamma_{A}}^{-}_{--}{T_{A}}^{-}_{+}&=0,\\\
\partial_{+}{T_{A}}^{+}_{-}+{\Gamma_{A}}^{+}_{++}{T_{A}}^{+}_{-}&=-\frac{H^{2}}{48\pi}\delta(-x_{\rm{I}}+L/2).\end{split}$
(32)
$\begin{split}\partial_{-}{T_{B}}^{-}_{+}+{\Gamma_{B}}^{-}_{--}{T_{B}}^{-}_{+}&=0,\\\
\partial_{+}{T_{B}}^{+}_{-}+{\Gamma_{B}}^{+}_{++}{T_{B}}^{+}_{-}&=\frac{H^{2}}{48\pi}\delta(x_{\text{{III}}}-L/2).\end{split}$
(33)

The solution of these equations are given by

$\displaystyle{T_{A}}^{-}_{+}$
$\displaystyle=-2\cosh^{2}\left[H(x_{\rm{I}}-L/2)\right]{T_{A}}_{++}$
$\displaystyle=f_{A}\left[x^{+}_{\rm{I}}\right]\cosh^{2}\left[H(x_{\rm{I}}-L/2)\right],$
(34) $\displaystyle{T_{A}}^{+}_{-}$
$\displaystyle=-2\cosh^{2}\left[H(x_{\rm{I}}-L/2)\right]{T_{A}}_{--}$
$\displaystyle=\left[\frac{H^{2}}{24\pi}\theta(-x_{\rm{I}}+L/2)+a\right]$
$\displaystyle\qquad\times
g_{A}\left[x^{-}_{\rm{I}}\right]\cosh^{2}\left[H(x_{\rm{I}}-L/2)\right],$ (35)
$\displaystyle{T_{B}}^{-}_{+}$
$\displaystyle=-2\cosh^{2}\left[H(x_{\text{{III}}}+L/2)\right]{T_{B}}_{++}$
$\displaystyle=f_{B}\left[x^{+}_{\text{{III}}}\right]\cosh^{2}\left[H(x_{\text{{III}}}+L/2)\right],$
(36) $\displaystyle{T_{B}}^{+}_{-}$
$\displaystyle=-2\cosh^{2}\left[H(x_{\text{{III}}}+L/2)\right]{T_{B}}_{--}$
$\displaystyle=\left[\frac{H^{2}}{24\pi}\theta(x_{\text{{III}}}+L/2)+b\right]$
$\displaystyle\qquad\times
g_{B}\left[x^{-}_{\text{{III}}}\right]\cosh^{2}\left[H(x_{\text{{III}}}+L/2)\right],$
(37)

where $a,b$ are integration constants, and $f_{j}\left[x^{+}_{j}\right]$ are
arbitrary functions of $x^{+}_{j}$, $g_{j}\left[x^{-}_{j}\right]$ are
arbitrary functions of $x^{-}_{j}$.

In the later section, we will calculate ${T_{\text{{III}}}}_{++}$ as the left-
moving Hawking radiation. As this quantity depends only on the coordinate
$x_{\text{{III}}}^{+}$, it is determined by the initial condition in region I
and the boundary conditions at $x=\pm L/2$, which are provided by the anomaly
equation.

###  III.2 Boundary conditions

We write the total EMT in region I$\,\cup\,$A as

${T_{\rm{Itot}}}^{\mu}_{\nu}={T_{\text{I}}}^{\mu}_{\nu}\left[1-\theta(-x_{\rm{I}}+L/2)\right]+{T_{\rm{A}}}^{\mu}_{\nu}\,\theta(-x_{\rm{I}}+L/2).$
(38)

Under the diffeomorphism transformation, the change of effective action $W$ is

$-\delta_{\lambda}W=\int
d^{2}x\sqrt{-g}\,\lambda^{\nu}\nabla_{\mu}{T_{\rm{Itot}}}^{\mu}_{\nu},$ (39)

where $\lambda^{\nu}$ are arbitrary variational parameters. To restore the
diffeomorphism covariance, each coefficients of $\lambda^{\nu}$ have to
vanish. Each coefficient is divergence of the EMT and calculated as

$\displaystyle\nabla_{\mu}{T_{\rm{Itot}}}^{\mu}_{+}$
$\displaystyle=\nabla_{\mu}{T_{\text{I}}}^{\mu}_{+}\left[1-\theta(-x_{\rm{I}}+L/2)\right]$
$\displaystyle\quad+\nabla_{\mu}{T_{\rm{A}}}^{\mu}_{+}\,\theta(-x_{\rm{I}}+L/2)$
$\displaystyle\quad+({T_{\rm{A}}}^{\mu}_{+}-{T_{\text{I}}}^{\mu}_{+})\partial_{\mu}\theta(-x_{\rm{I}}+L/2)$
$\displaystyle=\left(-\frac{H^{2}}{48\pi}+\frac{1}{2}{T_{\rm{A}}}^{-}_{+}-\frac{1}{2}{T_{\rm{I}}}^{-}_{+}\right)\delta(-x_{\rm{I}}+L/2),$
(40)

and

$\displaystyle\nabla_{\mu}{T_{\rm{Itot}}}^{\mu}_{-}$
$\displaystyle=\nabla_{\mu}{T_{\text{I}}}^{\mu}_{-}\left[1-\theta(-x_{\rm{I}}+L/2)\right]$
$\displaystyle\quad+\nabla_{\mu}{T_{\rm{A}}}^{\mu}_{-}\,\theta(-x_{\rm{I}}+L/2)$
$\displaystyle\quad+({T_{\rm{A}}}^{\mu}_{-}-{T_{\text{I}}}^{\mu}_{-})\partial_{\mu}\theta(-x_{\rm{I}}+L/2)$
$\displaystyle=\left(-\frac{1}{2}{T_{\rm{A}}}^{+}_{-}+\frac{1}{2}{T_{\rm{I}}}^{+}_{-}\right)\delta(-x_{\rm{I}}+L/2).$
(41)

To make the total EMT anomaly free, the coefficient of the delta function must
vanish and boundary conditions for $T_{++}$, $T_{--}$ at $x=L/2$ become

$\begin{split}{T_{\rm{A}}}_{++}-{T_{\rm{I}}}_{++}&=-\frac{H^{2}}{48\pi},\\\
{T_{\rm{A}}}_{--}-{T_{\rm{I}}}_{--}&=0,\end{split}$ (42)

By performing the same calculation in region III$\,\cup\,$B, we obtain
following condition at $x=-L/2$:

$\begin{split}{T_{\text{{III}}}}_{++}-{T_{B}}_{++}&=\frac{H^{2}}{48\pi},\\\
{T_{\text{{III}}}}_{--}-{T_{B}}_{--}&=0.\end{split}$ (43)

###  III.3 Calculation of Hawking radiation

From the general solution of the EMT obtained above and the boundary
conditions, we can calculate the flux in region III. Since (30), (35) and (37)
show that $T_{++}$ in each region does not depend on $x^{-}$, $T_{++}$ in each
region is determined once its value is specified at a certain $x^{+}$. We
impose the in-vacuum condition $T_{++}=0$ at $\mathcal{I}^{-}$ in region I.
Then, one finds

${T_{\rm{I}}}_{++}=-\frac{1}{2}F_{\rm{I}}\left[x^{+}_{\rm{I}}\right]=0.$ (44)

With the EMT in region I, using the matching condition (42), we obtain

${T_{\rm{A}}}_{++}={T_{\rm{I}}}_{++}-\frac{H^{2}}{48\pi}=-\frac{H^{2}}{48\pi}.$
(45)

This boundary condition is same as to fix the arbitary function
$f_{A}[x_{\text{I}}^{+}]$ in Eq. (34) as
$f_{A}[x_{\text{I}}^{+}]=H^{2}/(24\pi)$.

Next, to determine the relationship between ${T_{\rm{A}}}_{++}$ and
${T_{\rm{B}}}_{++}$, we consider a coordinate transformation. In the overlaped
region A$\cap$B, they are related by the coordinate transformation

$T_{B++}=\left(\frac{\partial x_{\rm{I}}^{+}}{\partial
x_{\text{{III}}}^{+}}\right)^{2}T_{A++}.$ (46)

Using the formula (17) ,

$\displaystyle\frac{\partial x_{\rm{I}}^{+}}{\partial x_{\text{{III}}}^{+}}=$
$\displaystyle\csc{\left[HL+\arccos{\left[\frac{\cosh{Hx_{\text{{III}}}^{+}}\sin{\frac{HL}{2}}+\sinh{Hx_{\text{{III}}}^{+}}}{\cosh{Hx_{\text{{III}}}^{+}}+\sin{\frac{HL}{2}}\sinh{Hx_{\text{{III}}}^{+}}}\right]}\right]}$
$\displaystyle\times\sqrt{\frac{\cos^{2}{\frac{HL}{2}}}{(\cosh{Hx_{\text{{III}}}^{+}}+\sin{\frac{HL}{2}}\sinh{Hx_{\text{{III}}}^{+}})^{2}}}~{},$
(47)

and ${T_{\rm{B}}}_{++}$ is determined as

$\displaystyle{T_{\rm{B}}}_{++}$
$\displaystyle=-\frac{H^{2}}{48\pi}\csc^{2}{\left[HL+\arccos{\left[\frac{\cosh{Hx_{\text{{III}}}^{+}}\sin{\frac{HL}{2}}+\sinh{Hx_{\text{{III}}}^{+}}}{\cosh{Hx_{\text{{III}}}^{+}}+\sin{\frac{HL}{2}}\sinh{Hx_{\text{{III}}}^{+}}}\right]}\right]}$
$\displaystyle\quad\times\frac{\cos^{2}{\frac{HL}{2}}}{(\cosh{Hx_{\text{{III}}}^{+}}+\sin{\frac{HL}{2}}\sinh{Hx_{\text{{III}}}^{+}})^{2}}.$
(48)

Finally, by using the boundary condition at $x_{\text{{III}}}=-L/2$ (43), we
obtain the flux in region III as a function of $x_{\text{{III}}}^{+}$:

$\displaystyle{T_{\rm{{III}}}}_{++}$
$\displaystyle=-\frac{H^{2}}{48\pi}\csc^{2}{\left[HL+\arccos{\left[\frac{\cosh{Hx_{\text{{III}}}^{+}}\sin{\frac{HL}{2}}+\sinh{Hx_{\text{{III}}}^{+}}}{\cosh{Hx_{\text{{III}}}^{+}}+\sin{\frac{HL}{2}}\sinh{Hx_{\text{{III}}}^{+}}}\right]}\right]}$
$\displaystyle\quad\times\frac{\cos^{2}{\frac{HL}{2}}}{(\cosh{Hx_{\text{{III}}}^{+}}+\sin{\frac{HL}{2}}\sinh{Hx_{\text{{III}}}^{+}})^{2}}+\frac{H^{2}}{48\pi}.$
(49)

Figure 4: $x_{\text{{III}}}^{+}$ dependence of ${T_{\rm{{III}}}}_{++}$ in
region III with different values of $L$ $(H=1)$. For
$x_{\text{{III}}}^{+}\rightarrow\infty$, ${T_{\rm{{III}}}}_{++}$ approaches to
$H^{2}/(48\pi)$ which corresponds to the outgoing thermal flux with the
Gibbons-Hawking temperature $T_{\text{H}}=H/(2\pi)$. For
$x_{\text{{III}}}^{+}<0$, the left-moving flux is negative which represents
negative energy density. For a finite negative $x_{\text{{III}}}^{+}$, the
flux negatively diverges.

Figure 4 shows ${T_{\rm{{III}}}}_{++}$ as a function of $x_{\rm{{III}}}^{+}$.
For $x_{\text{{III}}}^{+}>0$, ${T_{\rm{{III}}}}_{++}>0$ and this represents
existence of the left propagating positive energy flux.
${T_{\rm{{III}}}}_{++}$ approaches to $H^{2}/(48\pi)$ as
$x_{\text{{III}}}^{+}\rightarrow+\infty$; If we assume thermality of the flux,
it reproduces Gibbons-Hawking temperature $T_{\text{H}}=H/(2\pi)$ which is
associated with the future event horizon (black hole horizon) in de Sitter
region II. On the other hand, for $x_{\text{{III}}}^{+}<0$, the left
propagating flux becomes negative and diverges at a finite negative value of
$x_{\text{{III}}}^{+}$. This behavior of the EMT reflects existence of the
negative energy flux associated with the past event horizon $\mathcal{H^{-}}$
(white hole horizon) in region II (see Fig. 2).

The flux formula (49) obtained in this study is correct for the parameter
range $LH<\pi/2$. Otherwise, the vacuum condition imposed on $\mathcal{I}^{-}$
in region I does not determine the state in region III. Actually, in such a
case, regions A and B do not have overlap and matching procedure fails. As the
limiting case $LH\rightarrow 0$, the formula (49) predicts
${T_{\rm{{III}}}}_{++}\rightarrow 0$ and consistently reproduces the Minkowski
case with no particle creations.

##  IV Summary and Conclusion

In this paper, Hawking radiation in an analog de Sitter spacetime in the
quantum Hall system was calculated by obtaining the general solution of the
anomalous conservation law in each region and considering the boundary
conditions between region I and II, and between region II and III. In this
system, the effect of gravitational anomaly does not appear in the
conservation law of the EMT in the bulk regions, but the effect of anomaly
appears as the boundary conditions between bulk regions. Furthermore, in the
original calculation of Hawking radiation using gravitational anomaly, the
assumption of the existence of anomaly was necessary. However, in the QH
system, the calculation becomes exact because of the inherent property that
only chiral fields that exist.

Figure 5: Penrose diagram for non-etnernal case that corresponds to black hole
formation via gravitational collapse.

It has been found that Hawking radiation can be obtained by calculating
gravitational anomaly, but there are several points needed to be discussed. We
found a negatively divergent behavior of ${T_{\text{{III}}}}_{++}$ at a
negative finite value of $x_{\text{{III}}}^{+}$. This divergent behavior of
the EMT is related to the existence of the past horizon $\mathcal{H}^{-}$ in
region II, but detailed investigation of this issue has not been done. We
assumed that the analog de Sitter spacetime is time-symmetric, which makes
boundary conditions simple. However, in order to extract only the effect of
Hawking radiation from the future horizon, an analysis using the non-eternal
setup [2] (Fig. 5) is required. We left this for our next work.

###### Acknowledgements.

The authors thank Yuuki Sugiyama for valuable discussion on this subject. YN
was supported by JSPS KAKENHI (Grant No. JP22H05257) and MEXT KAKENHI Grant-
in-Aid for Transformative Research Areas A “Extreme Universe”(Grant No.
24H00956).

## References

  * Hotta  _et al._ [2022] M. Hotta, Y. Nambu, Y. Sugiyama, K. Yamamoto, and G. Yusa, Expanding edges of quantum Hall systems in a cosmology language: Hawking radiation from de Sitter horizon in edge modes, Phys. Rev. D 105, 105009 (2022), arXiv:2202.03731 .
  * Nambu and Hotta [2023] Y. Nambu and M. Hotta, Analog de Sitter universe in quantum Hall systems with an expanding edge, Phys. Rev. D 107, 085002 (2023), arXiv:2301.09270 .
  * Alvarez-Gaume and Witten [1984] L. Alvarez-Gaume and E. Witten, Gravitational Anomalies, Nucl. Phys. B 234, 269 (1984).
  * Bertlmann [1996] R. A. Bertlmann,  _Anomalies in quantum field theory_ (1996).
  * Bertlmann and Kohlprath [2001] R. A. Bertlmann and E. Kohlprath, Two-dimensional gravitational anomalies, Schwinger terms and dispersion relations, Annals Phys. 288, 137 (2001), arXiv:hep-th/0011067 .
  * Hawking [1975] S. W. Hawking, Particle Creation by Black Holes, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * Gibbons and Hawking [1977] G. W. Gibbons and S. W. Hawking, Action Integrals and Partition Functions in Quantum Gravity, Phys. Rev. D 15, 2752 (1977).
  * Christensen and Fulling [1977] S. M. Christensen and S. A. Fulling, Trace Anomalies and the Hawking Effect, Phys. Rev. D 15, 2088 (1977).
  * Birrell and Davies [1984] N. D. Birrell and P. C. W. Davies,  _Quantum fields in curved space_ (Cambridge University Press, 1984).
  * Parikh and Wilczek [2000] M. K. Parikh and F. Wilczek, Hawking radiation as tunneling, Phys. Rev. Lett. 85, 5042 (2000), arXiv:hep-th/9907001 .
  * Robinson and Wilczek [2005] S. P. Robinson and F. Wilczek, A Relationship between Hawking radiation and gravitational anomalies, Phys. Rev. Lett. 95, 011303 (2005), arXiv:gr-qc/0502074 .
  * Iso  _et al._ [2006a] S. Iso, H. Umetsu, and F. Wilczek, Anomalies, Hawking radiations and regularity in rotating black holes, Phys. Rev. D 74, 044017 (2006a), arXiv:hep-th/0606018 .
  * Iso  _et al._ [2006b] S. Iso, H. Umetsu, and F. Wilczek, Hawking radiation from charged black holes via gauge and gravitational anomalies, Phys. Rev. Lett. 96, 151302 (2006b), arXiv:hep-th/0602146 .
  * Murata and Soda [2006] K. Murata and J. Soda, Hawking radiation from rotating black holes and gravitational anomalies, Phys. Rev. D 74, 044018 (2006), arXiv:hep-th/0606069 .
  * Jiang  _et al._ [2007] Q.-Q. Jiang, S.-Q. Wu, and X. Cai, Hawking radiation from the dilatonic black holes via anomalies, Phys. Rev. D 75, 064029 (2007), [Erratum: Phys.Rev.D 76, 029904 (2007)], arXiv:hep-th/0701235 .
  * Akhmedova  _et al._ [2010] V. E. Akhmedova, T. Pilling, A. de Gill, and D. Singleton, Tunneling/WKB and anomaly methods for Rindler and de Sitter space-times, Theor. Math. Phys. 163, 774 (2010).
  * Zampeli [2013] A. Zampeli, Radiation from horizons, chirality and the principle of effective theory, J. Phys. Conf. Ser. 453, 012024 (2013).


††institutetext: Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS,
The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

# On the importance of factorization for fast binned likelihood inference

C.Jesús-Valls<EMAIL_ADDRESS>
###### Abstract

Likelihood-based inference, central in modern particle physics data analysis
requires the extensive evaluation of a likelihood function that depends on set
of parameters defined by the statistical model under consideration. If an
analytical expression for the likelihood can be defined from first principles
the procedure is computationally straightforward. However, most experiments
require approximating the likelihood numerically using large statistical
samples of synthetic events generated using Monte Carlo methods. As a result,
the likelihood consists of a comparison of the expected versus the observed
event rates in a collection of histogram bins, defining binned likelihood
functions. When this occurs, evaluating the likelihood function involves, on
each occasion, recalculating the prediction in those bins, increasing the
computational load of these analysis drastically. In this text, I highlight
the importance of identifying which are the unique event configurations in the
binned likelihood definition and I provide an exact formula to update the
event rate predictions utilizing the minimum number of necessary calculations
by means of factorization. The aim of the discussion is to decrease the
computational load of widespread high-energy physics analyses, leading to
substantial speed improvements and reduced carbon footprints.

##  1 Introduction

Likelihood-based inference is a mainstream statistical framework for the
analysis of real and simulated data in contemporary high energy physics (HEP)
Cowan (1998). The foundational assumption is that the data observables,
represented by $\vec{x}$, stem from a known probability distribution,
$f_{\vec{\theta}}$, determined by a set of parameters, $\vec{\theta}$.
Bayesian and frequentist methods are used to draw conclusions through the
study of the likelihood function, $\mathcal{L}$, which defines the probability
of observing the collected data for any given value of $\vec{\theta}$,
treating the data as fixed. To achieve this, one needs to calculate the
likelihood of $\vec{\theta}$ given the observed data $\vec{x}$, conventionally
expressed as $\mathcal{L}(\vec{\theta}\,|\,\vec{x})$ = $f_{\vec{\theta}}$. In
some occasions, the system under study is simple enough as to allow for
$f_{\vec{\theta}}$ to be evaluated directly from the analytical expressions
characterizing it. In most cases, however, the accurate consideration of
nontrivial theoretical quantities and experimental settings results in the
necessity of using Monte Carlo (MC) methods. MC predictions are built upon
finite samples of discrete data, resulting in three major consequences: 1) the
extensive use of binned likelihoods, 2) the necessity of event reweighting and
3) and a largely increased computational load to evaluate the likelihood
function. These three aspects play a central role in many studies in HEP, each
imposing limitations that demarcate what physics studies are possible in
practice. However, whereas the two first points are mainstream knowledge in
most experiments and are often addressed in the scientific literature, e.g.
Refs. Cowan  _et al._ (2011); Buckley  _et al._ (2011); Campbell  _et al._
(2022), the latter has been left out of the main discussion. In this article,
I comment on the importance of using factorization in likelihood-based tests
involving event reweighting and I illustrate generically how to implement such
factorization with the goal of increasing the speed and reducing the carbon
footprint Allen (2022) of computationally heavy and widespread analysis in
HEP.

##  2 About the use of binned likelihoods

To commence, let’s review why binned likelihood analysis are ubiquitous HEP
analysis. In HEP, information is typically grouped in events, $\vec{n}$, each
resulting from the observation of a fundamental interaction and consisting of
the collection of the information measured by a detector. Observables
$\vec{x}$ are built performing operations on the events’ information, namely
$\vec{x}=\mathcal{Q}(\vec{n})$. Given a probability distribution for those
observables, $f_{\vec{\theta}}(\vec{x})$, determined by a set of parameters,
$\vec{\theta}$, then it is possible to do inference on those parameters by
studying the likelihood function $\mathcal{L}$, defined as:

$\mathcal{L}(\vec{\theta}\,|\,\vec{x})=f_{\vec{\theta}}(\vec{x})=\prod_{i}f_{\vec{\theta}}(x_{i}).$
(1)

In general, analysis often consist of looking for the value of $\vec{\theta}$
that maximizes the above function, determining how likely are other values of
$\vec{\theta}$ given the observed data, and studying the overall suitability
of $f_{\vec{\theta}}(\vec{x})$ to be a faithful model representation of the
data. When the probability of observing any given $\vec{x}$ can be calculated,
studying Eq. 1 is computationally straightforward. However, in HEP experiments
calculating such probability is oftentimes analytically impractical and Monte
Carlo (MC) methods are used to sample synthetic events according to the
underlying model probability density function. As a result, evaluating the
likelihood infinitesimally is no longer possible and to overcome this
limitation the likelihood is build through the comparison of the expected and
observed distributions in finite regions in observable space, corresponding to
the bin content of histograms of arbitrary dimensionality. Consequently, HEP
experiments very often base their physics studies on binned likelihood
analyses using:

$\mathcal{L}(\vec{\theta}\,|\,\vec{N}_{\text{obs}})=\text{P}(\vec{N}_{\text{obs}}\,|\,\vec{\lambda}(\vec{\theta}))=\prod_{i}\text{P}(N^{i}_{\text{obs}}\,|\,\lambda^{i}(\vec{\theta})),$
(2)

where $\vec{N}_{\text{obs}}$ and $\vec{\lambda}(\vec{\theta})$ denote
respectively the collection of observed and expected events in all the
observables bins $i$ and where P is a discrete probability distribution. If
the counts in every bin are independent, as it is normally the case in HEP,
then $P$ is known to be characterized, asymptotically, by the Poisson
distribution. Since the sample of synthetic data is finite, however, the error
of the MC sample plays a role that can be relevant in some conditions.
Consequently, modified Poisson likelihoods that include these corrections are
available Barlow and Beeston (1993); Argüelles  _et al._ (2019).
Remarkably, regardless of the likelihood choice, in Eq. 2
$\vec{N}_{\text{obs}}$ is treated as fixed and therefore the complexity of the
calculation depends entirely on the difficulty of evaluating
$\vec{\lambda}(\vec{\theta})$.

##  3 Event reweighting

The calculation of $\vec{\lambda}(\vec{\theta})$ using Monte Carlo techniques
is among the most demanding computational tasks in HEP, often requiring in
large experiments the use of hundreds or thousands of CPU hours to realize a
single simulation, and resulting in huge computational demands every year
Belyaev  _et al._ (2017). Therefore, re-calculating
$\vec{\lambda}(\vec{\theta})$ from scratch for every value of $\vec{\theta}$
in order to analyze $\mathcal{L}$ is computationally prohibitive. Instead,
modern experiments use event re-weighting. To explain it, let’s introduce some
definitions that will become key later.
Consider the nominal expected event rate in the $i$-th bin as described by:

$\lambda_{i}(\vec{\theta}_{\text{nom}})=\sum^{\text{events}}_{j}w^{ij}_{\text{nom}}$
(3)

where $w^{ij}_{\text{nom}}$ is the nominal event weight for the $j$-th event
in the $i$-th bin. As any event can only be present in one bin:

$w^{ij}_{\text{nom}}=\begin{cases}0&\text{if the $j$-th event is not present
in the i-$th$ bin.}\\\ w^{k}_{\text{nom}}&\text{if the $k$-th event belongs to
the $i$-th bin.}\end{cases}$ (4)

From this definition, we can simplify equation 3 into:

$\lambda_{i}(\vec{\theta}_{\text{nom}})=\sum^{\text{i-events}}_{k}w^{k}_{\text{nom}}$
(5)

where "i-events" indicates that the sum happens exclusively over the indices
of those events in the $i$-th bin.
Let’s consider now the effect of choosing a value of $\vec{\theta}$ other than
nominal. Instead of generating new MC samples, the solution consists in
calculating event-by-event response functions $\vec{g}$ that re-weight the
contribution of each event as a function $\vec{\theta}$ with respect to their
nominal value:

$\lambda_{i}(\vec{\theta})=\sum^{\text{i-events}}_{k}(w^{k}_{\text{nom}}\times
g_{k}(\vec{\theta})).$ (6)

Notice that, by definition, $g_{k}(\vec{\theta}_{\text{nom}})=1$, which allows
to recover Eq. 5 from Eq. 10. The reweight functions are connected to the fit
parameters through functions that act on subsets of disjoint parameters,
associated to parameter indices denoted by $\alpha$. Namely:

$g_{k}(\vec{\theta})=\prod_{\alpha}^{\text{all disjoint
$\alpha$}}h^{\alpha}_{k}(\theta_{\alpha}).$ (7)

where

$h^{\alpha}_{k}(\theta_{\alpha})=\begin{cases}1&\text{if the $k$-th event is
unaffected by $\theta_{\alpha}$.}\\\ w^{\alpha}_{k}&\text{if the $k$-th event
is affected by $\theta_{\alpha}$.}\end{cases}$ (8)

With this in mind, Eq. 7 reads:

$g_{k}(\vec{\theta})=\prod_{\alpha}^{\text{k-relevant disjoint
$\alpha$}}h^{\alpha}_{k}(\theta_{\alpha}).$ (9)

where "k-relevant" indicates that the product happens exclusively for those
$\alpha$ indices relevant for the $k$-th event.
Now that the notation has been introduced, let’s clarify the meaning of the
former equations considering a set of examples below.

### Fully disjoint example with one parameter

Let’s consider the simplest possible physical model where we have a single
parameter. Then, $\alpha$ is trivially disjoint, and in Eq. 9 $\alpha$ is just
an integer representing the $0$-th element of the 1-dimensional parameter
vector $\vec{\theta}$.
Imagine that the parameter under consideration, $\theta_{0}$, accounts for the
detector mass. In this case, all events are affected by this parameter:

$\lambda_{i}(\vec{\theta})=\sum^{\text{i-events}}_{k}(w^{k}_{\text{nom}}\times
h^{0}_{k}(\theta_{0})),$ (10)

and, $h^{0}_{k}(\theta_{0})$ is equal for all the events. Therefore, if we
define this identical function by $h^{0}$, follows:

$\lambda_{i}(\vec{\theta})=\sum^{\text{i-events}}_{k}(w^{k}_{\text{nom}}\times
h^{0}_{k}(\theta_{0}))=h^{0}(\theta_{0})\times\sum^{\text{i-events}}_{k}w^{k}_{\text{nom}}=h^{0}(\theta_{0})\times\lambda_{i}(\vec{\theta}_{\text{nom}}).$
(11)

In this example, Eq. 11 shows that full factorization is possible. Since
$\theta_{0}$ accounts for the detector mass, $\theta_{0}$ should play the role
of a normalization parameter. Namely, $\theta_{0}$=1 means that the detector
mass is nominal, and $\theta_{0}=1.2$ implies an increase of 20% with respect
to that nominal value. More in general, $h^{0}(\theta_{0})=\theta_{0}$. Thus,
in Eq. 11, we observe that the result aligns with the logic: if the detector
mass increases by 20%, the event rate in all bins goes up by 20%.

### Fully disjoint example with two parameters

Let’s consider a more generic case where we have two parameters, $\theta_{0}$
and $\theta_{1}$. Let $\theta_{0}$ be the same parameter as in the previous
example and we introduce a new normalization parameter, $\theta_{1}$,
representing the cross section of a specific interaction channel and that
affects only those events generated through it. In this case, Eq. 9 reads:

$g_{k}(\vec{\theta})=\prod_{\alpha}^{\\{0,1\\}}h^{\alpha}_{k}(\theta_{\alpha})=\begin{cases}h^{0}(\theta_{0})\times
h^{1}(\theta_{1})&\text{if the $k$-th event is associated to $\theta_{1}$.}\\\
h^{0}(\theta_{0})&\text{otherwise.}\end{cases}$ (12)

Where we have used the fact that all $h^{0}_{k}$ and $h^{1}_{k}$ functions are
identical regardless of $k$ and their value is represented by $h^{0}$ and
$h^{1}$ respectively. As in the previous example, the equivalence of these
functions would allow to factorize the calculation. Before focusing on this
topic on the next section, let’s conclude the explanation with another
example.

### General example with joint and disjoint parameters

Finally, we focus on the general case. Consider that, in addition to the
former parameters $\theta_{0}$ and $\theta_{1}$, we account in our model for
the role of two flavor neutrino oscillations. To do so, we include two new
parameters $\theta_{2}$ and $\theta_{3}$ playing the role of the physics
parameters $\phi$ and $\Delta m^{2}$ respectively. It is important to recall
that the physics equation governing this process is:

$f^{2\nu}_{\text{osc}}(\phi,\Delta
m^{2}_{23},E_{\nu},L)=\sin^{2}{2\phi}\sin^{2}\left(1.27\frac{\Delta
m^{2}L}{E_{\nu}}\right).$ (13)

Thus, this new example has two major novelties. First, in contrast with
$\theta_{0}$ and $\theta_{1}$ that can be considered independently,
$\theta_{2}$ and $\theta_{3}$ are intertwined. Therefore, it does not make
sense for $\alpha$ to be either 2 or 3, rather, we define this index as a
joint index that we express as [2,3]. Secondly, Eq. 13 is, additionally a
function of the length $L$ from the neutrino production to the detection point
and its energy $E_{\nu}$. Whereas in most practical cases $L_{\text{{fixed}}}$
is a constant for all the events that depends on the experimental settings,
the neutrino energy is rarely monochromatic such that $E_{\nu}$ needs to be
considered on an event-by-event basis. Consequently,

$h^{[2,3]}_{k}(\theta_{2},\theta_{3})=f^{2\nu}_{\text{osc}}(\theta_{2},\theta_{3},E^{k}_{\nu},L_{\text{{fixed}}})$
(14)

is different for every $k$, alike $h^{0}_{k}$ and $h^{1}_{k}$ that are equal
for all events as discussed in the examples above.
Following this considerations, we have that Eq. 12 becomes:

$g_{k}(\vec{\theta})=\prod_{\alpha}^{\\{0,1,[2,3]\\}}h^{\alpha}_{k}(\theta_{\alpha})=\begin{cases}h^{0}(\theta_{0})\times
h^{1}(\theta_{1})\times h^{[2,3]}_{k}(\theta_{2},\theta_{3})&\text{case 1}\\\
h^{0}(\theta_{0})\times h^{1}(\theta_{1})&\text{case 2}\\\
h^{0}(\theta_{0})\times h^{[2,3]}_{k}(\theta_{2},\theta_{3})&\text{case 3}\\\
h^{0}(\theta_{0})&\text{otherwise.}\end{cases}$ (15)

Where

  * •

Case 1: The $k$-th event is affected by all weights.

  * •

Case 2: The $k$-th event is affected by all but oscillation weights. That
would be the case for events produced by neutral currents.

  * •

Case 2: The $k$-th event is affected by all $\theta_{1}$-related weights, i.e.
those not generated through the interaction channel associated to
$\theta_{1}$.

An important observation is that, in this example two new effects converge:
the existence of joint indices "[2,3]" and the necessity of using event-by-
event response functions imposed by its dependence with $E^{k}_{\nu}$. This is
a specific characteristic of this example, but not an universal condition. In
other words, it is entirely possible to encounter scenarios where a disjoint
parameter necessitates unique event-by-event response functions and,
conversely, to have models with joint indices with identical response
functions for all events.

##  4 Event factorization

Now that the use of event reweighting in the calculation of the binned
likelihood has been formalized and illustrated by means of examples, we can
turn our attention to the computational aspect of the problem.
Firstly, note that in the generic Eq. 5, one iterates over all of the event is
one bin. Thus, to update the value of all bins, one needs to iterate over the
full set of events. Considering now Eq. 7, for each event is necessary to
iterate over all $\alpha$ relevant for that event. In the worse computational
case, every $\alpha$ is relevant for each event. This leads to an upper bound
of the number of calculations $N^{\text{max}}_{c}$ to update the event rate:

$N^{\text{max}}_{c}=N_{\text{sim}}\times N_{\text{alpha}},$ (16)

Namely, the maximum number of operations corresponds to the product of the
number of all simulated events $N_{\text{sim}}$ with the number of unique
$\alpha$ configurations $N_{\alpha}$.
Secondly, let’s look at Eq. 11. In the extreme case where there is only
parameter, and its associated reweighting function acts identically on all
events, we can achieve full factorization, such that to update the event rate
in every bin, one needs to perform one operation. This settles the lower bound
of the number of calculations $N^{\text{min}}_{c}$, that consist of one
operation per bin:

$N^{\text{min}}_{c}=N_{\text{bins}}.$ (17)

Therefore, we conclude that in general the number of calculation, $N_{c}$,
necessary to update $\vec{\lambda}(\vec{\theta})$ must be:

$N_{\text{bins}}\leq N_{c}\leq N_{\text{sim}}\times N_{\text{alpha}}.$ (18)

It is worth noting that, for the overall consistency of the likelihood
definition, at least one MC event needs to be expected per bin, such that
$N_{\text{bins}}\leq N_{\text{sim}}$, and since the system must have at least
one $\theta$, $1\leq N_{\text{alpha}}$.
To recap, event reweight factorization allows to reduce the number of
intermediate calculations to update $\vec{\lambda}(\vec{\theta})$. In
realistic HEP cases, $N_{\text{bins}}$ spans from 1 to several hundreds or
thousands, $N_{\text{sim}}$ is often in the range of $10^{5}$–$10^{8}$ events,
and $N_{\text{alpha}}$ can range from 1 to hundreds. When considering these
numbers, it becomes clear that, systems allowing a high degree of
factorization can benefit enormously from factorizing the calculation
$\vec{\lambda}(\vec{\theta})$, dramatically reducing the number of necessary
calculations.

### Maximum factorization

The goal is to derive a generic expression that allows to utilize the maximum
possible amount of factorization in every system. Let’s start by combining Eq.
10 and Eq. 9:

$\lambda_{i}(\vec{\theta})=\sum^{\text{i-events}}_{k}(w^{k}_{\text{nom}}\times\prod_{\alpha}h^{\alpha}_{k}).$
(19)

Where for simplicity, $h^{\alpha}_{k}(\theta_{\alpha})$ is shortened into
$h^{\alpha}_{k}$. We have seen there are two possible types of events in each
bin: those that share response functions with other events, and therefore can
be factorized; and those that get unique event-by-event weights and can’t be
factorized. Let’s make this explicit:

$\lambda_{i}(\vec{\theta})=\sum^{\text{i-fact}}_{k}(w^{k}_{\text{nom}}\times\prod_{\alpha}h^{\alpha}_{k})+\sum^{\text{i-no-
fact}}_{k}(w^{k}_{\text{nom}}\times\prod_{\alpha}h^{\alpha}_{k}).$ (20)

The second sum can not be simplified, so let’s focus on the first. We note
that:

$\prod_{\alpha}h^{\alpha}_{k}=\begin{cases}1&\text{Trivial Case}\\\
h^{\alpha}&\text{Case A}\\\ h^{\alpha}h^{\beta}&\text{Case B}\\\
h^{\alpha}h^{\beta}h^{\gamma}&\text{Case C}\\\ \text{etc.}&\end{cases}$ (21)

Where:

  * •

Trivial Case: Events that get no weights from any parameters.

  * •

Case A: Events that get weights from a single $\alpha$.

  * •

Case B: Events that get weights from two distinct parameter indices, $\alpha$
and $\beta$.

  * •

Case C: Events that get weights from three distinct parameter indices,
$\alpha$, $\beta$ and $\gamma$.

One can continue this list straightforwardly including the succeeding cases D,
E, F, etc, if necessary. However, for the illustration here that is not
necessary, and we limit ourselves to products involving up to three distinct
parameter indices.
With the above considerations follows:

$\displaystyle\sum^{\text{i-fact}}_{k}(w^{k}_{\text{nom}}\times\prod_{\alpha}h^{\alpha}_{k})$
$\displaystyle=\sum^{\text{Trivial}}_{k}w^{k}_{\text{nom}}$
$\displaystyle+\sum_{\\{\alpha\\}}\left(h^{\alpha}\sum^{\text{Case
A${}_{\\{\alpha\\}}$ }}_{k}w^{k}_{\text{nom}}\right)$
$\displaystyle+\sum_{\\{\alpha,\beta\\}}\left(h^{\alpha}h^{\beta}\sum^{\text{Case
B${}_{\\{\alpha,\beta\\}}$ }}_{k}w^{k}_{\text{nom}}\right)$
$\displaystyle+\sum_{\\{\alpha,\beta,\gamma\\}}\left(h^{\alpha}h^{\beta}h^{\gamma}\sum^{\text{Case
C${}_{\\{\alpha,\beta,\gamma\\}}$ }}_{k}w^{k}_{\text{nom}}\right)+\text{etc.}$
(22)

Where curly brackets indicate unique sets without permutations of parameter
indices. If we multiply and divide by
$\lambda_{i}(\vec{\theta}_{\text{nom}})$, follows:

$\displaystyle\sum^{\text{i-fact}}_{k}(w^{k}_{\text{nom}}\times\prod_{\alpha}h^{\alpha}_{k})=$
$\displaystyle=\lambda_{i}(\vec{\theta}_{\text{nom}})\left(C_{T}+\sum_{\\{\alpha\\}}h^{\alpha}C_{\alpha}+\sum_{\\{\alpha\beta\\}}h^{\alpha}h^{\beta}C_{\alpha\beta}+\sum_{\\{\alpha\beta\\}}h^{\alpha}h^{\beta}h^{\gamma}C_{\alpha\beta\gamma}+\text{etc.}\right)$
(23)

Where, for the $i$-th bin, $C_{T}$ denotes the fraction of nominal events that
are trivial, each Cα denotes the fraction of nominal events that only get
weights from one specific $\alpha$, Cαβ denotes the fraction of nominal events
that only get weights from one specific unique combination of $\alpha$ and
$\beta$, etc.
Finally, let’s plug Eq. 4 in Eq. 20 to provide the final generic expression:

$\displaystyle\lambda_{i}(\vec{\theta})$ $\displaystyle=\sum^{\text{i-no-
fact}}_{k}(w^{k}_{\text{nom}}\times\prod_{\alpha}h^{\alpha}_{k})+$
$\displaystyle+\lambda_{i}(\vec{\theta}_{\text{nom}})\left(C_{T}+\sum_{\\{\alpha\\}}h^{\alpha}C_{\alpha}+\sum_{\\{\alpha\beta\\}}h^{\alpha}h^{\beta}C_{\alpha\beta}+\sum_{\\{\alpha\beta\\}}h^{\alpha}h^{\beta}h^{\gamma}C_{\alpha\beta\gamma}+\text{etc.}\right)$
(24)

For any system, Eq. 4 maximally factorizes the expression of the reweighted
event rate prediction in every bin, minimizing the number of calculations
necessary to evaluate the likelihood.

##  5 A toy problem

For illustration, let’s consider the following toy problem. Paired with the
explanation, I provide a publicly available
implementation111https://github.com/cesarjesusvalls/factorization_demo meant
to demonstrate explicitly how to put the practice the concepts discussed in
this article.
Imagine that we have an accelerator neutrino experiment, that consists of only
one Cherenkov detector and a neutrino beam. Consider the beam to be identical
to that of T2K Abe  _et al._ (2011) with our hypothetical Cherenkov detector
placed at a distance short enough as to safely neglect neutrino oscillations.
Then, we simulate MC samples using GENIE Andreopoulos  _et al._ (2010)
together with a public T2K flux release Abe  _et al._ (2023). Although this is
by no means necessary for this toy problem, to mimic a pseudo-realistic
physics case, the output true muon angle and momentum from GENIE is smeared by
10% and a simplified 1 muon ring selection is applied using the following
criteria: 1) the reconstructed muon momentum $p^{\text{reco}}_{\mu}$ satisfies
$0.25\leq p_{\mu}\leq 2$GeV/c 2) there are no $\pi^{0}$ in the final state abd
3) there are no charged pions with momentum above 250 MeV/c.

Figure 1: Illustration of the toy problem sample described in the text. The
cyan lines illustrate the edges of its 57 bins.

We bin the reconstructed muon angle and momentum distributions in fifty seven
2D regions. The binning, presented in Fig. 1, is chosen as to not be far off
from a realistic choice for a real experiment. Nevertheless, in the context of
this toy problem the binning choice can be considered as arbitrary. We then
consider several choices for $\vec{\theta}$:

  * •

Test_0: Only one normalization parameter modifies charged-current quasielastic
(CCQE) events get.

  * •

Test_1: For parameters, divided in three regions in true neutrino energy are
used, to mimic flux reweighting.

  * •

Test_2: Corresponds to all parameters (4 in total) from Test_0 and Test_1.

  * •

Test_2: An intermediate number of parameters (12 in total): Two cross section
normalization CCQE and 2p2h channels, and 10 flux regions.

  * •

Test_3: A large number of parameters (51 in total): All cross section modes
are considered, there are 11 in the sample under study, and 40 flux regions.

After successfully verifying that both schemes –w/ and w/o factorization– lead
to identical binned event rates for any $\vec{\theta}$, a series of speed test
were done. The results are summarizes in Table 1.

Test ID | U. Conf |  $\sum$ U. Conf. | W. Calc. | T. w/ Fact | T. w/o Fact | Speed Up
---|---|---|---|---|---|---
0 | 2 | 114 | 200992 | 115 $\pm$ 13 µs | 349 $\pm$ 5 ms |  $\times$3044
1 | 3 | 162 | 285460 | 165 $\pm$ 18 µs | 405 $\pm$ 5 ms |  $\times$2454
2 | 6 | 278 | 486452 | 263 $\pm$ 23 µs | 451 $\pm$ 10 ms |  $\times$1713
3 | 29 | 1027 | 525824 | 810 $\pm$ 41 µs | 460 $\pm$ 7 ms |  $\times$568
4 | 308 | 6394 | 570920 | 5411$\pm$ 234 µs | 474 $\pm$ 10 ms |  $\times$88
Table 1: Summary metrics for the different tests in the toy problem. From left
to right: The number of unique parameter configurations (U. Conf), the sum of
the number of unique parameter configurations in each bin ($\sum$ U. Conf.),
the number of weight calculations without using factorization (W. Calc.), the
time to calculate $\vec{\lambda}(\vec{\theta})$ with (T. w/ Fact) and without
(T. w/o Fact) factorization, and the speed up improvement calculated as the
ratio of the last two columns.

Figure 2: Execution times associated to the tests in Table 1 for event
reweights without (left) and with (right) the application of factorization in
the calculation.

As expected, the reweight time without factorization, remains on the same
order of magnitude for all tests as the leading contributor to the execution
time are operations that are performed once for all events, such as searching
the parameter indices associated to each event and updating the event weight.
There is a linear sub-leading time dependence with the number of response
function evaluations, necessary to calculate the weight updates, which
translates into a slightly increasing the execution time from Test_0 to
Test_4, as observed in the left panel of Fig. 2. In contrast, the execution
time with factorization depends directly on the number of bins and the number
of unique configurations in each bin. The execution time grows linearly with
the sum of the unique parameter configurations in all bins, as illustrated in
the right panel of Fig. 2. This markedly different computational behavior
results in drastically improved execution times for all tests, illustrating
the importance of using event factorization.

##  6 Conclusions

In this article the importance of factorizing unique parameter configurations
has been explained and a general formula for its application in any system
involving event reweight has been presented. For illustration, a publicly
available toy model has been prepared and used to report in the article
various metrics that exemplify the advantages of incorporating factorization
in the implementation of likelihood-based analysis software. Such statistical
analyses are very common in HEP, where the computation time often slows down
measurements within large experimental collaborations. Therefore, this article
aims to guide readers on effectively utilizing factorization in calculations
related to event reweighting, thereby accelerating studies and reducing its
associated computational carbon footprint.

##  7 Acknowledgments

The author acknowledges fruitful discussions with A. Blanchet, L. Berns and A.
Muñoz.

## References

  * Cowan (1998) G. Cowan,  _Statistical Data Analysis_ (Oxford University Press, 1998) "ISBN 978-0-19-850156-5".
  * Cowan  _et al._ (2011) G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Eur. Phys. J. C 71, 1554 (2011), [Erratum: Eur.Phys.J.C 73, 2501 (2013)], arXiv:1007.1727 [physics.data-an] .
  * Buckley  _et al._ (2011) A. Buckley _et al._ , Phys. Rept. 504, 145 (2011), arXiv:1101.2599 [hep-ph] .
  * Campbell  _et al._ (2022) J. M. Campbell _et al._ , in  _Snowmass 2021_ (2022) arXiv:2203.11110 [hep-ph] .
  * Allen (2022) M. Allen, Physics World 35, 46 (2022).
  * Barlow and Beeston (1993) R. J. Barlow and C. Beeston, Comput. Phys. Commun. 77, 219 (1993).
  * Argüelles  _et al._ (2019) C. A. Argüelles, A. Schneider, and T. Yuan, JHEP 06, 030, arXiv:1901.04645 [physics.data-an] .
  * Belyaev  _et al._ (2017) N. Belyaev, D. Krasnopevtsev, R. Konoplich, V. Velikhov, and A. Klimentov,  _High performance computing system in the framework of the Higgs boson studies_ , Tech. Rep. (ATL-COM-SOFT-2017-089, 2017).
  * Abe  _et al._ (2011) K. Abe _et al._ (T2K), Nucl. Instrum. Meth. A 659, 106 (2011), arXiv:1106.1238 [physics.ins-det] .
  * Andreopoulos  _et al._ (2010) C. Andreopoulos _et al._ , Nucl. Instrum. Meth. A 614, 87 (2010), arXiv:0905.2517 [hep-ph] .
  * Abe  _et al._ (2023) K. Abe _et al._ (T2K), Eur. Phys. J. C 83, 782 (2023), arXiv:2303.03222 [hep-ex] .