Pair Production of the Open Superstrings from the Parallel-dressed D3-branes
in the Compact Spacetime

Niloufar Barghi-Janyar and Davoud Kamani

Department of Physics, Amirkabir University of Technology (Tehran Polytechnic)
P.O.Box: 15875-4413, Tehran, Iran
e-mails<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

We employ the boundary state formalism to compute the pair creation rate of
the open superstrings from the interaction of two parallel D3-branes. The
branes live in the partially compact spacetime. In addition, they have been
dressed with the internal gauge potentials and the Kalb-Ramond field.

Keywords: Background field; Compact; Boundary state; D-brane; Amplitude; Pair
production; Decay rate.

##  1 Introduction

The D-branes are effective tools to study the string theory [1]. For
investigating the D-branes interaction an appropriate approach is the boundary
state formalism [2]-[14]. Since the boundary state encodes all properties of
the corresponding D-brane, it is a powerful tool for studying the various
configurations of the D-branes. For example, a phenomenon similar to the
Casimir effect, which is known as the open string pairs creation from the
branes system, is one of these setups. The pair production of the open strings
in the bosonic and superstring theories has been studied in various papers
[15]-[25]. The fundamental strings creation through the branes interaction
demonstrates the attractive force between the D-branes [19]-[22].

Using the boundary state formalism, we shall compute the rate of the open
superstring pair creation from the interaction of two parallel D3-branes. The
branes have been embedded in the partially compact spacetime. Besides, the
electric and magnetic fields live on them. Therefore, at first we shall
compute the closed superstring amplitude for a general branes dimension $p$,
and then, we rewrite the interaction amplitude for the case $p=3$. For the
small separation of the branes, the nature of the interaction (repulsion or
attraction) completely is ambiguous. However, the annulus open superstring
amplitude is an appropriate tool for describing such cases. Thus, from the
cylinder amplitude we shall obtain the latter one. The one-loop annulus
amplitude enables us to extract the rate of the open superstring pair
creation. In fact, production of the superstring pairs obviously indicates the
decay of the underlying system. As a special case, the pair production rate
for the particular electromagnetic fields on the branes will be calculated. We
shall see that the compactification, similar to the magnetic field,
drastically increases the pair production rate.

This paper is organized as follows. In Sec. 2, the boundary state, associated
with a partially compact D$p$-brane in the presence of the background fields,
will be introduced. In Sec. 3, the interaction amplitude of a system, which
consists of two parallel-dressed D$p$-branes, will be calculated. In Sec. 4,
the decay rate of a system of two D3-branes via the open superstring pair
production will be obtained. Sec. 5 is devoted to the conclusions.

##  2 The associated boundary state to the D$p$-brane

We begin with the following sigma-model action for the closed string

$\displaystyle S=$ $\displaystyle-$
$\displaystyle\frac{1}{4\pi\alpha^{\prime}}\int_{\Sigma}{\rm
d}^{2}\sigma\left(\sqrt{-h}h^{ab}G_{\mu\nu}\partial_{a}X^{\mu}\partial_{b}X^{\nu}+\varepsilon^{ab}B_{\mu\nu}\partial_{a}X^{\mu}\partial_{b}X^{\nu}\right)$
(2.1) $\displaystyle+$
$\displaystyle\frac{1}{2\pi\alpha^{\prime}}\int_{\partial\Sigma}{\rm
d}\sigma\left(A_{\alpha}\partial_{\sigma}X^{\alpha}\right).$

We consider a constant Kalb-Ramond background field $B_{\mu\nu}$. For the
internal $U(1)$ gauge potential $A_{\alpha}$, we apply the Landau gauge
$A_{\alpha}=-\frac{1}{2}F_{\alpha\beta}X^{\beta}$ whit the constant field
strength $F_{\alpha\beta}$. We employ the flat worldsheet $\Sigma$ with the
metric $h_{ab}=\eta_{ab}=\rm{diag}{(-1,1)}$, which lives in the flat spacetime
with the metric $G_{\mu\nu}=\eta_{\mu\nu}=\rm{diag}{(-1,1,\cdots,1)}$. Note
that the devoted indices to the D$p$-brane worldvolume directions are
$\alpha,\beta\in\\{0,\alpha_{\rm n},\beta_{\rm n},\alpha_{\rm c},\beta_{\rm
c}\\}$ with $\\{\alpha_{\rm n},\beta_{\rm n}\\}\cup\\{\alpha_{\rm
c},\beta_{\rm c}\\}=\\{1,2,\cdots,p\\}$. Some of the brane directions and some
of the normal directions to it are compacted on tori. Therefore, the subscript
“n” and “c” refer to the non-compact and compact directions, respectively. For
the perpendicular directions to the brane we shall use $i\in\\{i_{\rm
n},i_{\rm c}\\}$ with $\\{i_{\rm n}\\}\cup\\{i_{\rm c}\\}=\\{p+1,\cdots,9\\}$.

The bosonic portions of the boundary state equations are obtained by the
variation of above action with respect to $X^{\mu}(\sigma,\tau)$,

$\displaystyle\left(\partial_{\tau}X^{\alpha}+\mathcal{F}^{\alpha}\;_{\beta}\partial_{\sigma}X^{\beta}-B^{\alpha}\;_{i}\partial_{\sigma}X^{i}\right)_{\tau=0}|B_{x}\rangle=0,$
(2.2) $\displaystyle\left(X^{i}-y^{i}\right)_{\tau=0}|B_{x}\rangle=0,$

where
$\mathcal{F}_{\alpha\beta}=\partial_{\alpha}A_{\beta}-\partial_{\beta}A_{\alpha}-B_{\alpha\beta}$
shows the total field strength. On the closed string boundary we applied
$\left(\delta X^{i}\right)_{\tau=0}=\left(X^{i}-y^{i}\right)_{\tau=0}$. The
parameters $y^{i}$s represent the location of the D$p$-branes. The second
equation implies that $\partial_{\sigma}X^{i}$ vanishes on the boundary.
Hence, the last term of the first equation is removed.

The solution of the equation of motion of closed string is

$\displaystyle
X^{\mu}(\sigma,\tau)=x^{\mu}+2\alpha^{\prime}p^{\mu}\tau+2L^{\mu}\sigma+\frac{i}{2}\sqrt{2\alpha^{\prime}}\sum_{m\neq
0}\frac{1}{m}\left(\alpha^{\mu}_{m}\;e^{-2im\left(\tau-\sigma\right)}+\tilde{\alpha}^{\mu}_{m}\;e^{-2im\left(\tau+\sigma\right)}\right).$
(2.3)

For the non-compact and compact directions $L^{\mu}$ is zero and
$N^{\mu}R^{\mu}$, respectively. If the $x^{\mu}$-direction is compact it has
the radius $R^{\mu}$. Thus, around this direction, the closed string possesses
a winding number $N^{\mu}$. Note that we assume the time direction $x^{0}$ is
non-compact.

By combining Eqs. (2.2) and (2.3) the boundary state equations take the
features

$\displaystyle\left[\alpha^{\alpha}_{n}+\tilde{\alpha}^{\alpha}_{-n}-\mathcal{F}^{\alpha}\;_{\beta}\left(\alpha^{\beta}_{n}-\tilde{\alpha}^{\beta}_{-n}\right)\right]|B_{x}\rangle=0,$
(2.4)
$\displaystyle\left(\alpha^{\prime}p^{\alpha}+\mathcal{F}^{\alpha}\;_{\beta}L^{\beta}\right)|B_{x}\rangle=0,$
$\displaystyle\left(\alpha^{i}_{n}-\tilde{\alpha}^{i}_{-n}\right)|B_{x}\rangle=0,$
$\displaystyle\left(x^{i}-y^{i}\right)|B_{x}\rangle=0,$ $\displaystyle
L^{i}|B_{x}\rangle=0.$

The last equation implies that the closed string cannot wind around the
compact directions which are perpendicular to the brane.

The action (2.1) is invariant under the global worldsheet supersymmetry.
Therefore, by using the replacements
$\partial_{+}X^{\mu}\left(\sigma,\tau\right)\to
i\eta\Psi^{\mu}_{+}\left(\sigma,\tau\right)$ and
$\partial_{-}X^{\mu}\left(\sigma,\tau\right)\to-\Psi^{\mu}_{-}\left(\sigma,\tau\right)$
in the bosonic Eqs. (2.2) we acquire the fermionic parts of the boundary state
equations

$\displaystyle\left[\left(\Psi^{\alpha}_{-}-\mathcal{F}^{\alpha}\;_{\beta}\Psi^{\beta}_{-}\right)-i\eta\left(\Psi^{\alpha}_{+}+\mathcal{F}^{\alpha}\;_{\beta}\Psi^{\beta}_{+}\right)\right]_{\tau=0}|B_{\Psi},\eta\rangle=0,$
(2.5)
$\displaystyle\left(\Psi^{i}_{-}+i\eta\Psi^{i}_{+}\right)_{\tau=0}|B_{\Psi},\eta\rangle=0,$

where $\eta=\pm 1$ is preserved for applying the GSO projection. In terms of
the fermions oscillators these boundary conditions take the following forms

$\displaystyle\left(\Psi^{\alpha}_{q}-i\eta\mathcal{Q}^{\alpha}\;_{\beta}\tilde{\Psi}^{\beta}_{-q}\right)|B_{\Psi},\eta\rangle=0,$
(2.6)
$\displaystyle\left(\Psi^{i}_{q}+i\eta\tilde{\Psi}^{i}_{-q}\right)|B_{\Psi},\eta\rangle=0,$

where the index $q$ is integer (half-integer) for the R-R (NS-NS) sector. The
orthogonal matrix $\mathcal{Q}$ possesses the definition
$\mathcal{Q}=\left({\bf 1}-\mathcal{F}\right)^{-1}\left({\bf
1}+\mathcal{F}\right)$.

Now we solve Eqs. (2.4) and (2.6). For the bosonic part, the solution of
equations (2.4) is given by

$\displaystyle|B_{x}\rangle$ $\displaystyle=$
$\displaystyle\frac{T_{p}}{2}\sqrt{\det\left({\bf
1}-\mathcal{F}\right)}\prod_{i}\delta\left(x^{i}-y^{i}\right)|p^{i}=0\rangle$
(2.7) $\displaystyle\times$
$\displaystyle\exp{\left(-\sum^{\infty}_{m=1}\frac{1}{m}\alpha^{\mu}_{-m}\mathcal{S}_{\mu\nu}\tilde{\alpha}^{\nu}_{-m}\right)}|0\rangle_{\alpha}|0\rangle_{\tilde{\alpha}}\prod_{\alpha}|p^{\alpha}\rangle,$

in which
$\mathcal{S}^{\mu}\;_{\nu}=\left(\mathcal{Q}^{\alpha}\;_{\beta}\;,-\delta^{i}\;_{j}\right)$,
and $T_{p}$ refers to the tension of the D$p$-brane. The factor
$\sqrt{\det\left({\bf 1}-\mathcal{F}\right)}$ is originated from the path
integral in the presence of the Gaussian boundary action [26]-[28].

Eqs. (2.6) give the boundary state of the NS-NS and R-R sectors as in the
following

$\displaystyle|B_{\Psi},\eta\rangle_{\rm NS}$ $\displaystyle=$
$\displaystyle\exp\left(i\eta\sum_{r=1/2}^{\infty}b^{\mu}_{-r}\mathcal{S}_{\mu\nu}\tilde{b}^{\nu}_{-r}\right)|0\rangle,$
(2.8) $\displaystyle|B_{\Psi},\eta\rangle_{\rm R}$ $\displaystyle=$
$\displaystyle\left[\det\left({\bf
1}-\mathcal{F}\right)\right]^{-1/2}\;\exp\left(i\eta\sum_{n=1}^{\infty}d^{\mu}_{-n}\mathcal{S}_{\mu\nu}\tilde{d}^{\nu}_{-n}\right)|B_{\Psi},\eta\rangle^{(0)}_{\rm
R}.$ (2.9)

In comparison with the bosonic equation (2.7), the factor
$\sqrt{\det\left({\bf 1}-\mathcal{F}\right)}$ has been reversed. This is
because of the presence of the Grassmannian variables in computing the path
integral. The zero-mode part of the R-R sector possesses the form

$\displaystyle|B_{\Psi},\eta\rangle^{(0)}_{\rm
R}=\left(C\Gamma^{0}\Gamma^{1}\cdots\Gamma^{p}\;\frac{1+i\eta\Gamma_{11}}{1+i\eta}U\right)_{AB}|A\rangle|\tilde{B}\rangle,$
(2.10)

where A and B represent the 32-dimensional indices for the spinors,
$\Gamma^{\mu}s$ show the Dirac matrices in the 10-dimensional spacetime,
$|A\rangle|\tilde{B}\rangle$ is the vacuum of the R-R zero modes $d^{\mu}_{0}$
and $\tilde{d}^{\mu}_{0}$, $C$ defines the charge conjugation matrix, and the
matrix $U$ has the definition

$\displaystyle
U=;\exp\left(\frac{1}{2}\mathcal{F}_{\alpha\beta}\Gamma^{\alpha}\Gamma^{\beta}\right);.$
(2.11)

The conventional notation $;\;;$ implies that we should expand the exponential
with the convention that all $\Gamma$-matrices anticommute. Thus, for each
value of $p$, we receive a finite number of terms.

The total boundary state, corresponding to the D$p$-brane, in each sector is
given by the product of the matter- and ghost-parts

$|B,\eta\rangle_{\rm NS(R)}=|B_{x}\rangle|B_{\Psi},\eta\rangle_{\rm
NS(R)}|B_{\rm{gh}}\rangle|B_{\rm{sgh}},\eta\rangle_{\rm NS(R)}.$ (2.12)

The GSO-projected boundary states in the NS-NS and R-R sectors have the forms

$|B\rangle_{\rm NS}=\frac{1}{2}\left[|B,+\rangle_{\rm NS}-|B,-\rangle_{\rm
NS}\right],\quad\;\;\;|B\rangle_{\rm R}=\frac{1}{2}\left[|B,+\rangle_{\rm
R}+|B,-\rangle_{\rm R}\right].$ (2.13)

Since the conformal ghosts and super-conformal ghosts are not influenced by
the background fields we shall apply the standard form of their boundary
states.

##  3 The D$p$-branes interaction

The overlap of the boundary states, via the closed string propagator, gives
the tree-level interaction amplitude between two parallel-dressed branes
$\mathcal{A}_{\rm NS-NS(R-R)}=_{\rm NS(R)}\langle B_{1}|D|B_{2}\rangle_{\rm
NS(R)}$. The closed string propagator is as follows

$\displaystyle D$ $\displaystyle=$
$\displaystyle\frac{\alpha^{\prime}}{4\pi}\int_{|z|\leq 1}\frac{{\rm
d}^{2}z}{|z|^{2}}z^{L_{0}}\bar{z}^{\bar{L}_{0}},$ (3.1)

where $L_{0}$ and $\tilde{L}_{0}$ are the right- and left-moving total
Virasoro generators, including the bosonic, fermionic, conformal ghosts and
super-conformal ghosts parts. The total interaction amplitude, via the
exchange of the closed superstring, is the cylinder amplitude
$\mathcal{A}_{\rm{cylinder}}=\mathcal{A}_{\rm NS-NS}+\mathcal{A}_{\rm R-R}$.
Thus, we acquire

$\displaystyle\mathcal{A}_{\rm{cylinder}}$ $\displaystyle=$
$\displaystyle\frac{2^{\left[\frac{p+1}{2}\right]-3}\;T_{p}^{2}}{\left(2\pi\right)^{9-p}}\alpha^{\prime}\left[\det{\left({\bf
1}-\mathcal{F}_{1}\right)}\det{\left({\bf
1}-\mathcal{F}_{2}\right)}\right]^{1/2}\int^{\infty}_{0}{\rm
d}t\bigg{\\{}\left(\frac{\pi}{\alpha^{\prime}t}\right)^{d_{i_{\rm n}}/2}$
(3.2) $\displaystyle\times$ $\displaystyle\exp{\left(-\frac{Y^{2}_{\rm
n}}{4\alpha^{\prime}t}\right)}\prod_{i_{\rm
c}}\Theta_{3}\left(\frac{y_{1}^{i_{\rm c}}-y_{2}^{i_{\rm c}}}{2\pi R_{i_{\rm
c}}}\bigg{|}\frac{i\alpha^{\prime}t}{\pi R_{i_{\rm
c}}^{2}}\right)\left[\eta\left(it\right)\right]^{3\left[\frac{p+1}{2}\right]-12}$
$\displaystyle\times$ $\displaystyle\left(\sum_{\\{N^{\alpha_{\rm
c}}\\}}\bigg{[}\left(2\pi\right)^{p+1}\prod_{\alpha}\delta\left(p^{\alpha}_{1}-p^{\alpha}_{2}\right)\exp{\left(-\frac{t}{\alpha^{\prime}}l^{\alpha_{\rm
c}}l^{\beta_{\rm c}}\left(\delta_{\alpha_{\rm c}\beta_{\rm
c}}+\mathcal{F}^{\alpha}_{1\;\alpha_{\rm c}}\mathcal{F}_{2\;\alpha\beta_{\rm
c}}\right)\right)}\bigg{]}\right)$ $\displaystyle\times$
$\displaystyle\bigg{(}-\Theta_{2}(0|it)^{4-[(p+1)/2]}\prod^{[\frac{p-1}{2}]}_{a=0}\Theta_{2}(\nu_{a}|it)+\Theta_{3}(0|it)^{4-[(p+1)/2]}\prod^{[\frac{p-1}{2}]}_{a=0}\Theta_{3}(\nu_{a}|it)$
$\displaystyle\times$
$\displaystyle-\Theta_{4}(0|it)^{4-[(p+1)/2]}\prod^{[\frac{p-1}{2}]}_{a=0}\Theta_{4}(\nu_{a}|it)\bigg{)}\prod^{[\frac{p-1}{2}]}_{a=0}\frac{\sin(\pi\nu_{a})}{\Theta_{1}(\nu_{a}|it)}\bigg{\\}},$

where $l^{\alpha_{\rm c}}=N^{\alpha_{\rm c}}R^{\alpha_{\rm c}}$ is the
eigenvalue of the operator $L^{\alpha_{\rm c}}$, and $d_{i_{\rm n}}$ is the
dimension of the transverse non-compact directions. The transverse compact
directions induce the $\Theta_{3}$-function in the second line. Besides, the
$\Theta_{2}$-function comes from the R-R sector, and the $\Theta_{3}$\- and
$\Theta_{4}$-functions originate from the NS-NS sector. In addition, the
quantities $\\{e^{2\pi i\nu_{a}}|a=1,2,\cdots,[(p-1)/2]\\}$ are the
eigenvalues of the orthogonal matrix $\Lambda=Q_{1}^{\rm T}Q_{2}$. They
obviously depend on the magnetic and electric fields on the branes. The
integer part of $(p\pm 1)/2$ is denoted by $[(p\pm 1)/2]$. The branes distance
along the non-compact directions is $Y_{\rm n}=\left[\sum_{i_{\rm
n}}\left(y^{i_{\rm n}}_{1}-y^{i_{\rm n}}_{2}\right)^{2}\right]^{1/2}$. Note
that in the third line, the momentum components are extracted from the second
equation of Eq. (2.4), i.e.,
$p^{\alpha}_{1(2)}=-\frac{1}{\alpha^{\prime}}\mathcal{F}^{\alpha}_{1(2)\;\beta_{\rm
c}}N^{\beta_{\rm c}}R^{\beta_{\rm c}}$.

From now on, for simplification, we restrict our setup to the case $p=3$. This
dimension also is appropriate for relating the D-branes to the real world.
However, the amplitude for two parallel D3-branes takes the feature

$\displaystyle\mathcal{A}_{\rm{cylinder}}$ $\displaystyle=$
$\displaystyle\frac{T_{3}^{2}\alpha^{\prime}}{(2\pi)^{6}}\sqrt{\det({\bf
1}-\mathcal{F}_{1})\det({\bf 1}-\mathcal{F}_{2})}\int^{\infty}_{0}{\rm
d}t\bigg{\\{}\left(\frac{\pi}{\alpha^{\prime}t}\right)^{d_{i_{\rm n}}/2}$
(3.3) $\displaystyle\times$ $\displaystyle\exp{\left(-\frac{Y^{2}_{\rm
n}}{4\alpha^{\prime}t}\right)}\prod_{i_{\rm
c}}\Theta_{3}\left(\frac{y_{1}^{i_{\rm c}}-y_{2}^{i_{\rm c}}}{2\pi R_{i_{\rm
c}}}\bigg{|}\frac{i\alpha^{\prime}t}{\pi R^{2}_{i_{\rm
c}}}\right)\left[\eta(it)\right]^{-6}$ $\displaystyle\times$
$\displaystyle\left(\sum_{\\{N^{\alpha_{\rm
c}}\\}}\left(2\pi\right)^{4}\prod_{\alpha}\delta(p^{\alpha}_{1}-p^{\alpha}_{2})\exp{\left(-\frac{t}{\alpha^{\prime}}l^{\alpha_{\rm
c}}l^{\beta_{\rm c}}\left(\delta_{\alpha_{\rm c}\beta_{\rm
c}}+\mathcal{F}^{\alpha}_{1\;\alpha_{\rm c}}\mathcal{F}_{2\;\alpha\beta_{\rm
c}}\right)\right)}\right)$ $\displaystyle\times$
$\displaystyle\left[\;\Theta^{2}_{1}\left(\frac{\nu_{0}+\nu_{1}}{2}\bigg{|}it\right)\;\Theta^{2}_{1}\left(\frac{\nu_{0}-\nu_{1}}{2}\bigg{|}it\right)\right]\prod^{1}_{a=0}\frac{\sin{\left(\pi\nu_{a}\right)}}{\Theta_{1}\left(\nu_{a}|it\right)}\bigg{\\}},$

in which we employed the following identity for the theta-functions [29],

$\displaystyle-\Theta_{2}\left(0|it\right)^{2}\Theta_{2}\left(\nu_{0}|it\right)\Theta_{2}\left(\nu_{1}|it\right)+\Theta_{3}\left(0|it\right)^{2}\Theta_{3}\left(\nu_{0}|it\right)\Theta_{3}\left(\nu_{1}|it\right)$
$\displaystyle-\Theta_{4}\left(0|it\right)^{2}\Theta_{4}\left(\nu_{0}|it\right)\Theta_{4}\left(\nu_{1}|it\right)=2\Theta^{2}_{1}(\frac{\nu_{0}+\nu_{1}}{2}|it)\Theta^{2}_{1}(\frac{\nu_{0}-\nu_{1}}{2}|it).$
(3.4)

According to the orthogonality of the matrix $\Lambda$ and the characteristics
of its eigenvalues, one of the $\nu_{a}$s is pure imaginary, e.g.
$\nu_{0}=i\upsilon_{0}$. The range of the real $\nu_{a}$s and $\upsilon_{0}$
are $\nu_{a}\in[0,1)$ and $\upsilon_{0}\in(0,\infty)$, respectively. Now
consider small distance of the branes (which is corresponding to the small
value of “$t$”). Thus, the infinite product form of $\Theta_{1}(\nu_{a}|it)$
in the denominator of Eq. (3.3) gives the factor
$\prod^{\infty}_{n=1}[1-\cosh(2\pi\upsilon_{0})]$. Therefore, the sign of the
amplitude will be proportional to $(-1)^{\infty}$, which is ambiguous. This
demonstrates that the nature of the interaction (repulsion or attraction)
obviously is ambiguous. This ambiguity indicates a new phenomenon. Precisely,
the exponential factor $\exp{\left(-Y^{2}_{\rm n}/4\alpha^{\prime}t\right)}$
for the small $t$ and also the foregoing ambiguity play a significant role for
occurrence of this new phenomenon. In fact, this new phenomenon implies the
decay of the system via the open superstring pair production.

##  4 The rate of open superstring pair creation

The Jacobi transformation $t\to 1/t^{\prime}$ will be used to calculate the
open superstring one-loop interaction amplitude from the cylinder one. This
enables us to understand the nature of the foregoing physical phenomenon. In
fact, this new phenomenon comprises the open superstring pair production,
which we will be found later. However, the annulus amplitude possesses the
feature

$\displaystyle\mathcal{A}_{\rm annulus}$ $\displaystyle=$
$\displaystyle\frac{T_{3}^{2}\pi^{d_{i_{\rm
n}}/2}}{(2\pi)^{6}\alpha^{\prime(d_{i_{\rm n}}-2)/2}}\sqrt{\det({\bf
1}-\mathcal{F}_{1})\det({\bf 1}-\mathcal{F}_{2})}\int^{\infty}_{0}{\rm
d}t^{\prime}\bigg{\\{}$ (4.1) $\displaystyle\times$
$\displaystyle\exp{\left(-t^{\prime}Y^{2}_{\rm
n}/4\alpha^{\prime}\right)}t^{\prime(d_{i_{\rm n}}-8)/2}\;\prod_{i_{\rm
c}}\Theta_{3}\left(\frac{y_{1}^{i_{\rm c}}-y_{2}^{i_{\rm c}}}{2\pi R_{i_{\rm
c}}}\bigg{|}\frac{i\alpha^{\prime}}{t^{\prime}\pi R^{2}_{i_{\rm c}}}\right)$
$\displaystyle\times$ $\displaystyle\left[\sum_{\\{N^{\alpha_{\rm
c}}\\}}\left(2\pi\right)^{4}\prod_{\alpha}\delta(p^{\alpha}_{1}-p^{\alpha}_{2})\exp{\left(-\frac{1}{t^{\prime}\alpha^{\prime}}l^{\alpha_{\rm
c}}l^{\beta_{\rm c}}\left(\delta_{\alpha_{\rm c}\beta_{\rm
c}}+\mathcal{F}^{\alpha}_{1\;\alpha_{\rm c}}\mathcal{F}_{2\;\alpha\beta_{\rm
c}}\right)\right)}\right]$ $\displaystyle\times$
$\displaystyle\frac{\sin(\pi\nu_{0})\sin(\pi\nu_{1})}{\sinh{(\pi\nu_{0}t^{\prime})}\sinh{(\pi\nu_{1}t^{\prime})}}\left[\cosh{(\pi\nu_{0}t^{\prime})}-\cosh{(\pi\nu_{1}t^{\prime})}\right]^{2}\prod^{\infty}_{m=1}\mathcal{C}_{m}(t^{\prime})\bigg{\\}},$

where $\mathcal{C}_{m}(t^{\prime})$ is given by

$\displaystyle\mathcal{C}_{m}(t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{\left[1-2\cosh{(\pi t^{\prime}(\nu_{0}+\nu_{1}))}e^{-2m\pi
t^{\prime}}+e^{-4m\pi t^{\prime}}\right]^{2}}{{\left(1-e^{-2m\pi
t^{\prime}}\right)^{4}}\left[1-2\cosh{(2\pi t^{\prime}\nu_{0})}e^{-2m\pi
t^{\prime}}+e^{-4m\pi t^{\prime}}\right]}$ $\displaystyle\times$
$\displaystyle\frac{\left[1-2\cosh{(\pi t^{\prime}(\nu_{0}-\nu_{1}))}e^{-2m\pi
t^{\prime}}+e^{-4m\pi t^{\prime}}\right]^{2}}{1-2\cosh{(2\pi
t^{\prime}\nu_{1})}e^{-2m\pi t^{\prime}}+e^{-4m\pi t^{\prime}}}.$

To acquire this feature of Eq. (4.1), we applied the infinite product form of
the $\Theta_{1}$\- and Dedekind $\eta$-functions, and also their Jacobi
transformation formulas, e.g. see Ref. [29].

Now, by applying $\nu_{0}=i\upsilon_{0}$ in the amplitude (4.1) we can obtain
its imaginary part, which exhibits the pair production of the open
superstrings, and hence the decay of the underlying system. Presence of the
factor $\sin{(\pi\upsilon_{0}t^{\prime})}$ in the denominator indicates the
simple poles in the positive direction of the $t^{\prime}$-axis, i.e.,
$t^{\prime}_{k}=k/\upsilon_{0}$ where $k$ is any positive integer number.
Hence, each pole separately leads to the creation of an open superstring pair,
and consequently deterioration of the branes system. The decay rate per unit
volume of the D3-branes is defined by $\mathcal{W}_{3,3}=-2{\rm
Im}\mathcal{A}_{\rm annulus}/V_{4}$. Thus, we obtain

$\displaystyle\mathcal{W}_{3,3}$ $\displaystyle=$
$\displaystyle\frac{2T_{3}^{2}\pi^{d_{i_{\rm
n}}/2}}{\upsilon_{0}V_{4}(2\pi)^{6}\alpha^{\prime(d_{i_{\rm
n}}-2)/2}}\sqrt{\det({\bf 1}-\mathcal{F}_{1})\det({\bf 1}-\mathcal{F}_{2})}$
(4.3) $\displaystyle\times$
$\displaystyle\sum_{k}\bigg{\\{}\left(-1\right)^{k+1}\exp{\left(-kY^{2}_{\rm
n}/4\upsilon_{0}\alpha^{\prime}\right)}\left(\frac{k}{\upsilon_{0}}\right)^{(d_{i_{\rm
n}}-8)/2}$ $\displaystyle\times$ $\displaystyle\prod_{i_{\rm
c}}\Theta_{3}\left(\frac{y_{1}^{i_{\rm c}}-y_{2}^{i_{\rm c}}}{2\pi R_{i_{\rm
c}}}\bigg{|}\frac{i\alpha^{\prime}\upsilon_{0}}{k\pi R^{2}_{i_{\rm
c}}}\right)\bigg{[}\sum_{\\{N^{\alpha_{\rm
c}}\\}}\left(2\pi\right)^{4}\prod_{\alpha}\delta(p^{\alpha}_{1}-p^{\alpha}_{2})$
$\displaystyle\times$
$\displaystyle\exp{\left(-\frac{\upsilon_{0}}{k\alpha^{\prime}}l^{\alpha_{\rm
c}}l^{\beta_{\rm c}}\left(\delta_{\alpha_{\rm c}\beta_{\rm
c}}+\mathcal{F}^{\alpha}_{1\alpha_{\rm c}}\mathcal{F}_{2\alpha\beta_{\rm
c}}\right)\right)}\bigg{]}\frac{\sinh(\pi\upsilon_{0})\sin(\pi\nu_{1})}{\sinh{(\frac{k\pi\nu_{1}}{\upsilon_{0}})}}$
$\displaystyle\times$
$\displaystyle\left[\left(-1\right)^{k}-\cosh{\left(\frac{k\pi\nu_{1}}{\upsilon_{0}}\right)}\right]^{2}\mathcal{C^{\prime}}_{k}(\upsilon_{0},\nu_{1})\bigg{\\}},$

in which the worldvolume of each D3-brane indicated by $V_{4}$, and
$\mathcal{C^{\prime}}_{k}(\upsilon_{0},\nu_{1})$ has the definition

$\displaystyle\mathcal{C^{\prime}}_{k}(\upsilon_{0},\nu_{1})$ $\displaystyle=$
$\displaystyle\prod^{\infty}_{m=1}\frac{\left[1-2\left(-1\right)^{k}\cosh(\frac{k\pi\nu_{1}}{\upsilon_{0}})e^{-2m\pi
k/\upsilon_{0}}+e^{-4m\pi k/\upsilon_{0}}\right]^{4}}{\left(1-e^{-2m\pi
k/\upsilon_{0}}\right)^{6}\left[1-2\cosh(\frac{2k\pi\nu_{1}}{\upsilon_{0}})e^{-2m\pi
k/\upsilon_{0}}+e^{-4m\pi k/\upsilon_{0}}\right]}.$ (4.4)

###  4.1 An special case

As we see, the above result is very complicated to accurately express the
influence of the parameters of the setup on the decay rate. Thus, as an
example, we shall rewrite the foregoing rate for the special matrices
$\mathcal{F}_{1}$ and $\mathcal{F}_{2}$. Therefore, let us choose the matrix
$\mathcal{F}_{1}$ $(\mathcal{F}_{2})$, with the nonzero electric field $E$
$(E^{\prime})$ and magnetic field $B$ $(B^{\prime})$ for the first (the
second) D3-brane, as in the following

$\displaystyle\mathcal{F}^{\alpha}_{1\;\beta}=\begin{pmatrix}0&-E&0&0\\\
-E&0&0&0\\\ 0&0&0&B\\\ 0&0&-B&0\\\ \end{pmatrix}.$ (4.5)

The eigenvalues $\lambda_{0}$ and $\lambda_{1}$ of the matrix $\Lambda$
satisfy the equations

$\displaystyle\lambda_{0}+\lambda^{-1}_{0}$ $\displaystyle=$
$\displaystyle\frac{2\left[1+4EE^{\prime}+E^{\prime 2}+E^{2}(1+E^{\prime
2})\right]}{(1-E^{2})(1-E^{\prime 2})},$
$\displaystyle\lambda_{1}+\lambda^{-1}_{1}$ $\displaystyle=$
$\displaystyle\frac{2\left[1+4BB^{\prime}-2B^{2}-2B^{\prime
2}(1-B^{2})\right]}{(1+B^{2})(1+B^{\prime 2})}.$ (4.6)

Thus, $\upsilon_{0}$ and $\nu_{1}$ possess the forms

$\displaystyle\tanh(\pi\upsilon_{0})$ $\displaystyle=$
$\displaystyle\frac{E+E^{\prime}}{1+EE^{\prime}},$
$\displaystyle\tan(\pi\nu_{1})$ $\displaystyle=$
$\displaystyle\frac{B^{\prime}-B}{1+BB^{\prime}}.$ (4.7)

For this configuration, we receive the decay rate as in the following

$\displaystyle\mathcal{W}_{3,3}$ $\displaystyle=$
$\displaystyle\frac{2T_{3}^{2}\pi^{d_{i_{\rm
n}}/2}\;\sinh(\pi\upsilon_{0})}{(2\pi)^{6}\alpha^{\prime(d_{i_{\rm
n}}-2)/2}}\left(1+B^{2}\right)\left(1-E^{2}\right)\sum_{k}\bigg{\\{}\frac{1}{k}\left(-1\right)^{k+1}$
(4.8) $\displaystyle\times$ $\displaystyle\exp{\left(-kY^{2}_{\rm
n}/4\upsilon_{0}\alpha^{\prime}\right)}\left(\frac{k}{\upsilon_{0}}\right)^{(d_{i_{\rm
n}}-8)/2}\prod_{i_{\rm c}}\Theta_{3}\left(\frac{y_{1}^{i_{\rm
c}}-y_{2}^{i_{\rm c}}}{2\pi R_{i_{\rm
c}}}\bigg{|}\frac{i\alpha^{\prime}\upsilon_{0}}{k\pi R^{2}_{i_{\rm
c}}}\right)$ $\displaystyle\times$
$\displaystyle\Omega\left(E,E^{\prime},B,B^{\prime},\frac{\upsilon_{0}}{k},R_{1},R_{2},R_{3}\right)\left[\left(-1\right)^{k}-1\right]^{2}\mathcal{C^{\prime}}_{k}(\upsilon_{0},0)\bigg{\\}},$
$\displaystyle\Omega$ $\displaystyle=$
$\displaystyle\frac{\left(2\pi\right)^{4}}{V_{4}}\sum^{\infty}_{n=-\infty}\sum^{\infty}_{s=-\infty}\sum^{\infty}_{m=-\infty}\bigg{\\{}\delta[(E-E^{\prime})nR_{1}/\alpha^{\prime}]$
(4.9) $\displaystyle\times$
$\displaystyle\delta[(B-B^{\prime})sR_{2}/\alpha^{\prime}]\delta[(B-B^{\prime})mR_{3}/\alpha^{\prime}]$
$\displaystyle\times$
$\displaystyle\exp{\left(-\frac{\upsilon_{0}}{k\alpha^{\prime}}\left[n^{2}R^{2}_{1}(1-E^{2})+\left(s^{2}R^{2}_{2}+m^{2}R^{2}_{3}\right)(1+B^{2})\right]\right)}\bigg{\\}}.$

The $\Omega$-function comes from the third line of Eq. (4.1). According to the
delta-functions we applied $E=E^{\prime}$ and $B=B^{\prime}$. Hence, the
second equation of Eqs. (4.7) gives $\nu_{1}=0$. Note that $-1<E<1$, and
$\mathcal{C^{\prime}}_{k}(\upsilon_{0},0)$ is calculated from Eq. (4.4) in
terms of $\nu_{1}=0$. Besides, in this special configuration, the tangential
directions $\\{x^{1},x^{2},x^{3}\\}$ are compact.

Now we compute the open superstring pair production rate. According to [30],
it is specified with the $k=1$ term. Thus, we acquire

$\displaystyle\mathcal{W}^{(k=1)}_{3,3}$ $\displaystyle=$
$\displaystyle\frac{8T_{3}^{2}\pi^{d_{i_{\rm
n}}/2}\sinh(\pi\upsilon_{0})}{(2\pi)^{6}\alpha^{\prime(d_{i_{\rm
n}}-2)/2}}\left(1+B^{2}\right)\left(1-E^{2}\right)\exp{\left(-Y^{2}_{\rm
n}/4\upsilon_{0}\alpha^{\prime}\right)}$ (4.10) $\displaystyle\times$
$\displaystyle{\upsilon_{0}}^{(8-d_{i_{\rm n}})/2}\;\prod_{i_{\rm
c}}\Theta_{3}\left(\frac{y_{1}^{i_{\rm c}}-y_{2}^{i_{\rm c}}}{2\pi R_{i_{\rm
c}}}\bigg{|}\frac{i\alpha^{\prime}\upsilon_{0}}{\pi R^{2}_{i_{\rm
c}}}\right)\Theta_{3}\left(0\bigg{|}\frac{i\upsilon_{0}R^{2}_{1}\left(1-E^{2}\right)}{\pi\alpha^{\prime}}\right)$
$\displaystyle\times$
$\displaystyle\Theta_{3}\left(0\bigg{|}\frac{i\upsilon_{0}R^{2}_{2}\left(1+B^{2}\right)}{\pi\alpha^{\prime}}\right)\Theta_{3}\left(0\bigg{|}\frac{i\upsilon_{0}R^{2}_{3}\left(1+B^{2}\right)}{\pi\alpha^{\prime}}\right)\mathcal{C^{\prime}}_{1}\left(\upsilon_{0},0\right).$

For the very small electric fields on the D3-branes, i.e. $|E|\ll 1$, Eqs.
(4.7) gives $\upsilon_{0}\approx 2|E|/\pi+\mathcal{O}(E^{3})$. This implies
that $\upsilon_{0}\ll 1$, and hence
$\mathcal{C^{\prime}}_{1}(\upsilon_{0},0)\approx 1$. Besides, for the finite
values of the phrases $R^{2}_{1}$ and $R^{2}_{2,3}(1+B^{2})$, the four
$\Theta_{3}$s in Eq. (4.10) reduce to 1. Finally, in this limit, the rate of
the open superstring pair production takes the form

$\displaystyle\mathcal{W}^{(k=1)}_{3,3}=\frac{T_{3}^{2}\pi^{d_{i_{\rm
n}}/2}}{8\pi^{5}\alpha^{\prime(d_{i_{\rm
n}}-2)/2}}\;\left(1+B^{2}\right)\left(\frac{2|E|}{\pi}\right)^{(10-d_{i_{\rm
n}})/2}\exp{\left(-\frac{\pi Y^{2}_{\rm n}}{8\alpha^{\prime}|E|}\right)}.$
(4.11)

We observe that the pair creation rate has been prominently enhanced by the
magnetic field. Though the two factors which contain the electric field are
very small, by increasing the magnetic field the rate can be adjusted to any
desirable value. In this case, the compactification radii $R_{2}$ and $R_{3}$
should be sufficiently small such that $|E|R^{2}_{2,3}(1+B^{2})\ll 1$. As we
see, the presence of the electric fields drastically plays the main role in
the production of open superstring pairs, as expected.

Let us rewrite this equation in the non-compact spacetime. Therefore, we
should replace $Y^{2}_{\rm n}\to
Y^{2}=\sum^{9}_{i=4}\left(y^{i}_{1}-y^{i}_{2}\right)^{2}$ and $d_{i_{\rm
n}}\to 6$. Accordingly, we receive

$\displaystyle\mathcal{W}^{(k=1)}_{3,3}\bigg{|}_{\rm non-
compact}=\frac{T_{3}^{2}}{8\pi^{2}\alpha^{\prime
2}}\;\left(1+B^{2}\right)\left(\frac{2|E|}{\pi}\right)^{2}\exp{\left(-\frac{\pi
Y^{2}}{8\alpha^{\prime}|E|}\right)}.$ (4.12)

The exponential factors of Eqs. (4.11) and (4.12) elaborate that by
compactifying some of the spatial directions the pair production rate
possesses an enhancement.

##  5 Conclusions

In the context of the superstring theory we introduced a boundary state,
associated with a D$p$-brane, in which some of the longitudinal and transverse
directions have been compacted on tori. The brane has been dressed with
electric and magnetic fields. The closed superstring amplitude was also
presented. Afterward, we changed the cylinder amplitude to the annulus one.

The 1-loop annulus amplitude enabled us to compute the decay rate of two
parallel D3-branes system through the pair production of the open
superstrings. Since the background fields on the branes are arbitrary and
different, the foregoing rate found a generalized form. However, for clarity
we chose special electric and magnetic fields on the branes. Thus, we observed
that for small electric fields the pair production rate is independent of the
branes distance in the compact subspace. In addition, in comparison with the
non-compact spacetime, this specific configuration demonstrated that the
compactification enhances the value of the pair creation rate. Besides, the
magnetic field also enhances the rate. Moreover, for acquiring a nonzero pair
production rate, the electric fields are clearly needed to induce a polarized
region between the branes.

## References

  * [1] J. Polchinski, Phys. Rev. Lett. 75 (1995) 4724; J. Polchinski, “String Theory”, (Cambridge University Press, Cambridge, 1998), Volumes I and II. For a review, also see J. Polchinski, hep-th/9611050.
  * [2] P. Di Vecchia, M. Frau, I. Pesando, S. Sciuto, A. Lerda and R. Russo, Nucl. Phys. B 507 (1997) 259-276.
  * [3] M. Billo, P. Di Vecchia and D. Cangemi, Phys.Lett. B 400 (1997) 63-70.
  * [4] M. Billo, P. Di Vecchia, M. Frau, A. Lerda, I. Pesando, R. Russo and S. Sciuto, Nucl.Phys. B 526 (1998) 199.
  * [5] F. Hussain, R. Iengo and C. Nunez, Nucl. Phys. B 497 (1997) 205.
  * [6] O. Bergman, M. Gaberdiel and G. Liftschytz, Nucl. Phys. B 509 (1998) 194.
  * [7] P. Di Vecchia, A. Liccardo, in: YITP Workshop on Developments in Superstring and M Theory, Kyoto, Japan, 1999, in: YITP Proceedings Series, vol. 4, 1999; P. Di Vecchia, M. Frau, A. Lerda, A. Liccardo, Nucl. Phys. B 565 (2000)397-426; P. Di Vecchia, A. Liccardo, R. Marotta, F. Pezzella, JHEP09 (2004) 050; Int. J. Mod. Phys. A 20 (2005) 5699-4795.
  * [8] C. Bachas, Phys. Lett. B 374 (1996) 37.
  * [9] C.G. Callan, I.R. Klebanov, Nucl. Phys. B 465 (1996) 473.
  * [10] S. Gukov, I.R. Klebanov, A.M. Polyakov, Phys. Lett. B 423 (1998) 64.
  * [11] H. Arfaei and D. Kamani, Phys. Lett. B 452 (1999) 54-60, arXiv:hep-th/9909167; D. Kamani, Phys. Lett. B 487 (2000) 187-191, arXiv:hep-th/0010019; E. Maghsoodi and D. Kamani, Nucl. Phys. B 922 (2017) 280-292, arXiv:1707.08383 [hep-th]; F. Safarzadeh-Maleki and D. Kamani, Phys. Rev. D 90 (2014) 107902, arXiv:1410.4948 [hep-th]; M. Saidy-Sarjoubi and D. Kamani, Phys. Rev. D 92 (2015) 046003, arXiv:1508.02084 [hep-th]; D. Kamani, Europhys. Lett. 57 (2002) 672-676, arXiv:hep-th/0112153.
  * [12] F. Safarzadeh-Maleki and D. Kamani, Phys. Rev. D 89 (2014) 026006, arXiv:1312.5489 [hep-th]; D. Kamani, Mod. Phys. Lett. A 17 (2002) 237-243, arXiv:hep-th/0107184; S. Teymourtashlou and D. Kamani, Eur. Phys. J. C 81 (2021) 761, arXiv:2108.10164 [hep-th]; D. Kamani, Eur. Phys. J. C 26 (2002) 285-291, arXiv:hep-th/0008020.
  * [13] H. Arfaei and D. Kamani, Nucl. Phys. B 561 (1999) 57-76, arXiv:hep-th/9911146; D. Kamani, Nucl. Phys. B 601 (2001) 149-168, arXiv:hep-th/0104089; N. Barghi-Janyar and D. Kamani, Int. J. Mod. Phys. A 38 (2023) 2350021, arXiv:2302.03659 [hep-th].
  * [14] D. Kamani, Ann. Phys. 354 (2015) 394-400, arXiv:1501.02453 [hep-th]; H. Arfaei and D. Kamani, Phys. Lett. B 475 (2000) 39-45, arXiv:hep-th/9909079; D. Kamani, Mod. Phys. Lett. A 15 (2000) 1655-1664, arXiv:hep-th/9910043.
  * [15] A. Hanany and E. Witten, Nucl. Phys. B 492 (1997) 152.
  * [16] T. Sato, Nucl.Phys. B 682 (2004) 117-149.
  * [17] T. Kitao, N. Ohta, J.G. ZhouPhys. Lett. B 428 (1998) 68-74.
  * [18] T. Takayanagi, JHEP 0002 (2000) 040.
  * [19] U.H. Danielsson, G. Ferretti, I.R. Klebanov, Phys. Rev. Lett. 79 (1997) 1984-1987.
  * [20] U.H. Danielsson, G. Ferretti, Nucl. Phys. Proc. Suppl. 68 (1998) 78-83.
  * [21] N. Ohta, T. Shimizu and J. Zhou, Phys. Rev. D 5 (1998) 2040-2044.
  * [22] T. Matsuo and T. Yokono, Mod. Phys. Lett. A 14 (1999) 1175-1182.
  * [23] Q. Jia, J. X. Lu, Z. Wu and X. Zhu, Nucl. Phys. B 953 (2020) 114947.
  * [24] C. Bachas and M. Porrati, Phys. Lett. B 296 (1992) 77.
  * [25] C.P. Burgess, Nucl. Phys. B 294 (1987) 427.
  * [26] M. Frau, I. Pesando, S. Sciuto, A. Lerda and R. Russo, Phys. Lett.B 400 (1997) 52.
  * [27] E. Bergshoeff, E. Sezgin, C.N. Pope and P.K. Townsend, Phys. Lett. B 188 (1987) 70.
  * [28] C.G. Callan, C. Lovelace, C.R. Nappi and S.A. Yost, Nucl. Phys. B 308 (1988) 221-284; Nucl. Phys. B 288 (1987) 525-550; Nucl. Phys. B 293 (1987).
  * [29] E.T. Whittaker and G.N. Watson, “A Course of Modern Analysis”, Cambridge University Press, doi: 10.1017/CBO9780511608759.
  * [30] A.I. Nikishov, Nucl. Phys. B 21 (1970) 346; Sov. Phys. JETP 30, 660 (1970).


Fitting a Collider in a Quantum Computer: Tackling the Challenges of Quantum
Machine Learning for Big Datasets

Miguel Caçador Peixoto1, Nuno Filipe Castro1, 2, Miguel Crispim Romão1, 3,

Maria Gabriela Jordão Oliveira1 and Inês Ochoa4

1 LIP – Laboratório de Instrumentação e Física Experimental de Partículas,
Escola de Ciências, Campus de Gualtar, Universidade do Minho, 4701-057 Braga,
Portugal

2 Departamento de Física, Escola de Ciências, Campus de Gualtar, Universidade
do Minho, 4701-057 Braga, Portugal

3 Department of Physics and Astronomy, University of Southampton, SO17 1BJ
Southampton, United Kingdom

4 LIP – Laboratório de Instrumentação e Física Experimental de Partículas, Av.
Prof. Gama Pinto, 2, 1649-003 Lisboa, Portugal

## Abstract

Current quantum systems have significant limitations affecting the processing
of large datasets with high dimensionality, typical of high energy physics. In
the present paper, feature and data prototype selection techniques were
studied to tackle this challenge. A grid search was performed and quantum
machine learning models were trained and benchmarked against classical shallow
machine learning methods, trained both in the reduced and the complete
datasets. The performance of the quantum algorithms was found to be comparable
to the classical ones, even when using large datasets. Sequential Backward
Selection and Principal Component Analysis techniques were used for feature’s
selection and while the former can produce the better quantum machine learning
models in specific cases, it is more unstable. Additionally, we show that such
variability in the results is caused by the use of discrete variables,
highlighting the suitability of Principal Component analysis transformed data
for quantum machine learning applications in the high energy physics context.

##  1 Introduction

The Standard Model of Particle Physics (SM) provides a remarkable description
of the fundamental constituents of matter and their interactions, being in
excellent agreement with the collider data accumulated so far. Nonetheless,
there are still important open questions, unaddressed by the SM, such as
gravity, dark matter, dark energy, or the matter-antimatter asymmetry in the
universe [1], motivating a comprehensive search program for new physics
phenomena beyond the SM (BSM) at the Large Hadron Collider (LHC) at CERN.

The search for BSM phenomena at colliders poses specific challenges in data
processing and analysis, given the extremely large datasets involved and the
low signal to background ratios expected. In this context, the analysis of the
collision data obtained by the LHC experiments often relies on machine
learning (ML), a field in computer science that can harness large amounts of
data to train generalizable algorithms for a variety of applications [2, 3],
such as classification tasks. These techniques have shown an outstanding
ability to find correlations in high-dimensional parameter spaces to
discriminate between potential signal and background processes. They are known
to scale with data, and usually rely on a large number of learnable parameters
to achieve their remarkable performance.

In order to train these large models, classical111 _Classical_ is used
throughout the paper as opposed to _quantum_ machine learning. machine
learning (CML) takes advantage of hardware accelerators, such as graphics
processing units (GPUs), for efficient, parallel, and fast matrix
multiplications. On the other hand, a new class of hardware is becoming
available, with the advent of noisy intermediate-scale quantum (NISQ)
computing devices. This accelerated the development of new quantum algorithms
targeted at exploiting the capacity and feasibility of this new technology for
ML applications.

Quantum machine learning (QML) is an emerging research field aiming to use
quantum circuits to tackle ML tasks. One of the motivations for using this new
technology in high energy physics (HEP) relates to the intrinsic properties of
quantum computations, namely representing the data in a Hilbert space where
the data can be in a superposition of states or in entangled states, which can
allow to explore additional information in data analysis and, eventually,
contribute to better classification of HEP events, namely in the context of
the search for BSM phenomena. Recently, this new technology has been applied
to HEP problems such as event reconstruction [4, 5, 6, 7, 8, 9, 10] and
classification [11, 12, 13, 14, 15, 16, 17, 18, 19].

Despite the promising potential of quantum computation, NISQ processors have
important limitations, such as the qubit quality (_i.e._ the accuracy with
which it is possible to execute quantum gates), the qubit lifetime and the
limited depth of quantum circuits, since for large circuits the noise
overwhelms the signal [20, 21]. This necessarily limits the complexity of the
circuits and the size of the datasets used to train them.

In this paper we perform a systematic comparison of the performance of QML and
shallow CML algorithms in HEP. The choice to focus on shallow methods rather
than state-of-the-art architectures based on deep neural networks is to
provide a fair comparison between methodologies, since neural networks are
known to require large datasets (both in terms of sample size and dimension)
to achieve good performance, something that is not feasible with current
quantum computers. By choosing CML algorithms suited for smaller datasets, we
will add to the on-going discussion regarding potential advantages of quantum
computing by comparing QML and CML in the same footing.

The use of QML algorithms in this context is studied by targeting a common
binary classification task in HEP: classifying a BSM signal against SM
background. A benchmark BSM signal leading to the $Zt$ final state is
considered, in events with multiple leptons and $b$-tagged jets, which can be
used to achieve a reasonable signal to background ratio. Variational quantum
classifiers (VQC) are trained and optimized via a grid search. The use of
reduced data is explored, considering both the number of features and the
number of events, via different strategies: ranking of features, data
transformations aiming for a richer reduced set of features, use of random
samples, and choice of representative data samples.

##  2 Quantum Machine Learning

The QML algorithms are implemented using a quantum circuit, _i.e._ a
collection of quantum gates applied to an $n$-qubit quantum state, followed by
a measurement (or multiple measurements) that represent the output of the
circuit. In order to implement a learning algorithm, the quantum circuit can
be parameterized with parameters that can be learned by confronting the
measurement to a loss function.

QML is effectively an extension of CML techniques to the Hilbert space, where
instead of representing data as vectors in a high-dimensional real space, we
encode it in state vectors of a Hilbert space. A QML algorithm, such as a
quantum neural network, can be implemented using the quantum equivalent of a
perceptron, one of the building blocks of CML. A problem arises from the
realization that the activation functions used in CML can not be expressed
using a linear operation, which is inherently required from the quantum
evolution of a state. Ideas have been proposed to imitate an activation
function in the quantum space, [22, 23] but, in the current paper, only
variational quantum classifiers [24, 25] are used for binary classification.

A VQC is a parameterized quantum circuit, a circuit type containing adjustable
gates with tunable parameters. These gates are a universal set of quantum
gates and, in the current study, rotation [$R_{X}(w)$, $R_{Y}(w)$, $R_{Z}(w)$]
and CNOT gates are used222Even if, in general, the phase shift gate $P(w)$
should be included, this gate does not change the final outcome (_i.e._ it
does not impact probabilities), so it can be discarded..

The considered VQC pipeline used has the following components:

  * •

Data Embedding: the numerical vector $X$ representing the classical
information is converted to the quantum space with the preparation of an
initial quantum state, $\ket{\psi_{X}}$, which represents a HEP event.

  * •

Quantum circuit: a unitary transformation $U(w)$, parameterized by a set of
free parameters $w$, is applied to the initial quantum state $\ket{\psi_{X}}$.
This produces the final state $\ket{\psi_{X}^{\prime}}=U(w)\ket{\psi_{X}}$.

  * •

Final prediction: a measurement of an observable is performed in one of the
qubits of the state $\ket{\psi_{X}^{\prime}}$, which will give the prediction
of the model for the task at hand. The training of a VQC aims to find the best
set of parameters $w$ to match the event labels to the prediction.

###  2.1 Data Embedding

Before passing the data through the VQC, the preparation of the initial
quantum state $\ket{\psi_{X}}$ is required. This is called data embedding, and
there are a number of proposals to perform this step [26]. The angle embedding
was tested and preliminary results found to lead to a better performance than
amplitude embedding, as reported before [17]. This embedding was therefore the
choice for the studies reported in this paper.

For an $N$-dimensional vector of classical information,
$X=(x_{1},x_{2},...,x_{N})$, the state entering the VQC will be defined via a
state preparation circuit applied to the initial ground state
$\ket{0}^{\bigotimes N}$. The information contained in $X$ is embedded as
angles: these are the values used in rotation gates applied to each qubit,
thus requiring $N$ qubits for embedding $N$ features from the original
dataset.

In the current study, the embedding is done using rotations around the
$x$-axis on the Bloch sphere, thus defining the quantum state embedded with
the classical information as:

$\ket{\psi_{X}}=\bigotimes_{i=1}^{N}R_{X}(x_{i})=\bigotimes_{i=1}^{N}\left[\cos\left(\frac{x_{i}}{2}\right)\ket{0}-i\sin\left(\frac{x_{i}}{2}\right)\ket{1}\right]\,\,,$
(1)

where $R_{X}(x)=e^{-ix\hat{\sigma}_{x}}$ and $\hat{\sigma}_{x}$ is a Pauli
operator. In this embedding each of the considered features of the original
dataset is required to be bound between $[-\pi,\pi]$.

###  2.2 Model Circuit

The model circuit is the key component of the VQC and includes the learnable
set of parameters. It is defined by a parameterized unitary circuit $U(w)$,
with $w$ being the set of tunable parameters, which will evolve a quantum
state embedded with classical information $\psi_{X}$ into the final state
$\psi_{X}^{\prime}$.

Analogously to the architecture of a classical neural network, the model
circuit is formed by layers. Each layer is composed of an assemblage of
rotation gates applied to each qubit in the system, followed by a set of CNOT
gates.

A rotation gate, $R$, is designed to be applied to one single qubit and rotate
its state. It is composed by 3 learnable parameters: $\phi,\theta,\omega$,
which enables the gate to rotate any arbitrary state to any location on the
Bloch sphere.

$R(\phi,\theta,\omega)=RZ(\omega)\,RY(\theta)\,RZ(\phi)=\begin{bmatrix}e^{-i(\phi+\omega)/2}\cos(\theta/2)&-e^{i(\phi-\omega)/2}\sin(\theta/2)\\\
e^{-i(\phi-\omega)/2}\sin(\theta/2)&e^{i(\phi+\omega)/2}\cos(\theta/2)\end{bmatrix}$
(2)

Since all the learnable parameters of the VQC are contained inside the
rotation gates, and each gate has 3 parameters, the shape of the weight vector
is $w\in\mathbb{R}^{n\times l\times 3}$, where $n$ is the number of qubits of
the current system and $l$ is the number of layers in the network. As
mentioned in the previous section, $n$ will depend on the number of features
in the data and $l$ is a hyper-parameter (HP) to be tuned.

After rotating the qubits’ state, a collection of CNOT gates will be applied
to entangle the qubits. The CNOT gate is a 2-qubit gate with no learnable
parameters. It will flip the state of the so-called target-qubit, based on the
value of the control-qubit, and it is usually represented by having two inputs
as such: CNOT(control-qubit, target-qubit). Given the number of qubits, the
CNOT arrangement is implemented as detailed in Algorithm 1.

Algorithm 1 CNOT Arrangement

$n\geq 2$, $n$ being the number of qubits.

if $n==2$ then

$CNOT(1,0)$

else

for $qubit\leftarrow 0$ to $n-1$ do

if $qubit==n-1$ then

$CNOT(qubit,0)$

else

$CNOT(qubit,qubit+1)$

end if

end for

end if

###  2.3 Final prediction

The output of the model is obtained by measuring the expectation value of the
Pauli $\hat{\sigma}_{z}$ operator in one of the qubits of the final state
$\psi_{X}^{\prime}$. An example of the implementations of a VQC is represented
in Figure 1.

Figure 1: An example circuit for the VQC architecture used. It is comprised of
2 layers and 3 features as input. The three main stages of a QML model can be
seen: embedding of the data, passing the data through the model circuit, and
the measurement of the outcome.

##  3 Classical Machine Learning Methods

Shallow CML methods are used to provide a baseline comparison to the QML
models. The specific methods chosen for the comparison are Logistic Regression
(LR) and Support Vector Machines (SVM), with these algorithms being trained
with the same data as the QML algorithms.

All the classical methods were implemented using scikit-learn [27] library
and, if not specified otherwise, the default parameters were used.

###  3.1 Logistic Regression

Logistic Regression is one of the simplest ML models and can be formulated as
one of the basic building blocks of a neural network - a single-layer
perceptron. The goal is to find the best set of weights $w$ that fit the data
$x$:

$\hat{y}(w,b,x)=\sigma(w\cdot x+b)\,,$ (3)

where $\hat{y}$ is the probability of an event to belong to class $1$, $w$ and
$b$ are learnable parameters, and $\sigma$ is the sigmoid function.

The learning process is guided by minimizing the loss function, which in our
case is the binary cross-entropy:

$L=-\mathbb{E}_{x}[y\log(\hat{y})+(1-y)\log(1-\hat{y})]\,,$ (4)

where $y$ is the binary label of whether the event is of the class signal or
not, and $\mathbb{E}_{x}$ is the expectation value over the training data,
obtaining using the event weights corresponding to each signal and background
process.

###  3.2 Support Vector Machine

An SVM classifier is trained by finding the hyperplane that best separates two
classes of data in the hyperspace of features. It does so by using support
vectors, which are the data points from the two classes closer to the
hyperplane, influencing the position and orientation of the hyperplane.

The loss function of an SVM revolves around the goal of maximizing the margin,
_i.e_ the distance between the hyperplane and the nearest data point from
either class. In other words, the goal is to find the hyperplane with the
greatest possible margin between itself and any point within the training set,
giving a greater chance of new data being classified correctly.

Just like the Logistic Regression, the base SVM classifier can only learn a
linear boundary condition. However, classification problems are rarely simple
enough for it to be separable using a hyperplane, thus usually requiring a
non-linear separation. SVM can do this by transforming the data using a non-
linear function, named kernel, after which it can be split by a hyperplane.
For this implementation, the radial-basis function (RBF) was used as kernel.
This endows the SVM with a non-linear mapping where it better separates the
two classes using a hyperplane.

##  4 Dataset

The dataset used in this work[28] is comprised of simulated events of $pp$
collisions at 13 TeV, in final states with 2 leptons, at least 1 $b$-jet, at
least 1 large-$R$ jet and large scalar sum of transverse333The transverse
plane is defined with respect to the proton colliding beams. momentum
($p_{T}$) of all reconstructed particles in the event ($H_{T}$ > 500 GeV).
Such basic selection corresponds to a topology commonly used in different
searches for BSM events at the LHC [29]. The dominant SM background for this
topology, $Zb\bar{b}$, and the BSM signal corresponding to $t\bar{t}$
production with one of the top-quarks decaying via a flavour changing neutral
current decay $t\to qZ$ ($q=c,u$) [30], were considered. Such signal was
chosen given the kinematic similitude to the background, thus providing a good
benchmark for the present study.

Both samples were generated with MadGraph5 2.6.5 [31] and Pythia 8.2 [32], and
the detector was simulated using Delphes 3 [33] with the default CMS card.
Jets were clustered using the anti-$k_{t}$ algorithm [34] with $R$-parameters
of 0.5 and 0.8 (the latter for the large-$R$ jets).

The following features were used for training of both the classical and
quantum machine learning algorithms:

  * •

$(\eta,\phi,p_{T},m,b$-tag) of the 5 leading jets, ordered by decreasing
$p_{T}$, with $b$-tag being a Boolean variable indicating if the jet is
identified as originating from a $b$-quark by the $b$-tagging algorithm
emulated by Delphes;

  * •

$(\eta,\phi,p_{T},m)$ of the leading large-$R$ jet;

  * •

$N$-subjettiness of leading large-$R$ jet, $\tau_{n}$ with $n=1,...,5$ [35].

  * •

$(\eta,\phi,p_{T})$ of the 2 leading leptons (electrons or muons);

  * •

transverse momentum ($\cancel{E}_{T}$) and $\phi$ of the missing transverse
energy;

  * •

multiplicity of jets, large-$R$ jets, electrons and muons;

  * •

$H_{T}$.

The proportion of signal and background events was kept the same as the
original simulated data during training, being $13\%$ and $87\%$ respectively.
Additionally, the Monte Carlo weights, corresponding to the theoretical
prediction for each process at target luminosity of 150 fb-1, were taken into
account in the evaluation of all the considered metrics and loss functions.

##  5 Feature Selection

As described in the previous section, a total of 47 features are available for
training. Considering the type of data embedding chosen, 47 qubits would be
needed to train a VQC using all the dataset features. Such number of qubits is
impractical given the currently available quantum computers and thus it is not
feasible to train a VQC using all the features in our dataset. For the
purposes of the current study, quantum computers with only 5 qubits were
considered and two methods for feature selection were implemented: principal
component analysis (PCA) and sequential feature selection (SFS).

A relative comparison of the best 5 features 444The area under the curve (AUC)
of the receiver operating characteristic curve (ROC) is considered as metric
for these comparisons. is shown in Table 1 while the best performance obtained
with state-of-the-art CML methods without any features or data points
restrictions can be seen in Figure 2.

Feature | AUC
---|---
$\cancel{E}_{T}$ | 0.817
Lepton1 $p_{T}$ | 0.692
Lepton2 $p_{T}$ | 0.649
large-$R$ jet $m$ | 0.609
large-$R$ jet $\tau_{1}$ | 0.576
Table 1: Top 5 features ranked by their AUC Score on the training dataset.

Figure 2: Obtained ROC curve and respective AUC score on the test dataset when
training an Boosted Decision Tree, implemented with  _xgboost_ [36] using the
full set of features and data points. The classifier has an identical
configuration as the one described in subsection 5.1.

###  5.1 Sequential Feature Selection

SFS algorithms are a widely used family of greedy search algorithms used for
automatically selecting a subset of features that is most relevant to the
problem. This is achievable by removing or adding one feature at a time based
on the classifier performance until a feature subset of the desired size, $k$,
is reached.

There are different variations of SFS algorithms but for the current paper,
the Sequential Backward Selection (SBS) algorithm was chosen. This algorithm
starts with the full set of features ($n=47$) and, at each iteration, it
generates all possible feature subsets of size $n-1$ and trains a ML model for
each one of the subsets. The performance is subsequently evaluated and the
feature that is absent from the subset of features with the highest
performance metric is removed. This process is iterated until the feature
subset contains $k$ features.

This technique was used to find subsets of 1 to 5 features. The ML model
assisting the SBS was a boosted decision tree (BDT) with a maximum number of
estimators set at $100$ and a learning rate of $1\times 10^{-5}$. The
considered loss function was a logistic regression for binary classification
and the AUC score was used as evaluation metric. The BDT was implemented using
xgboost [36] and the SBS algorithm using mlxtend [37]. The selected features
for the different values of $k$ is shown in Table 2 and the AUC scores for
each feature in Table 3.

$k$ | Selected Features
---|---
1 | $\cancel{E}_{T}$
2 |  $\cancel{E}_{T}$, Number of muons
3 |  $\cancel{E}_{T}$, Number of muons, Jet1 $b$-tag
4 |  $\cancel{E}_{T}$, Number of muons, Jet1 $b$-tag, Jet2 $p_{T}$
5 |  $\cancel{E}_{T}$, Number of muons, Jet1 $b$-tag, Jet2 $p_{T}$, large-$R$ $\tau_{3}$
Table 2: List of the features selected by the SBS algorithm for $k=1,...,5$. Feature | AUC
---|---
$\cancel{E}_{T}$ | 0.817
Number of muons | 0.534
Jet1 $b$-tag | 0.418
large-R jet $\tau_{3}$ | 0.316
Jet2 $p_{T}$ | 0.313
Table 3: Features selected by the SBS Algorithm and their respective AUC Score
on the training dataset.

###  5.2 Principal Component Analysis

The PCA transforms a highly correlated, high-dimensional dataset and into a
new one with reduced dimensionality and uncorrelated features, by rotating the
dataset in the direction of the eigenvectors of the dataset covariance matrix.
In the present paper, the PCA was performed only to remove the correlation
between the features, maintaining the same dimensionality as the original
data. The PCA transformation was learned from the training dataset and then
applied to all datasets. When training a VQC for a specific number of
features, the PCA components were ranked by AUC score and thus selected from
the highest to the lowest. This is done by introducing a priority queue,
_i.e._ if training a model using two features is desired, the 2 top-ranked PCA
components will be selected. The scikit-learn PCA implementation was used and
the obtained 5 better components are shown in Table 4.

PCA Component | AUC
---|---
Component 3 | 0.726
Component 14 | 0.606
Component 5 | 0.565
Component 41 | 0.563
Component 43 | 0.560
Table 4: Top 5 PCA components obtained with the training dataset, ranked by
their AUC.

##  6 Dataset Size Reduction

The present paper addresses the use of reduced datasets to overcome the
limitation of NISQ processors while minimizing the loss of information and
thus avoiding a performance loss of the QML algorithms in the HEP context. The
primary method used for this purpose in the current study is KMeans, where the
$k^{th}$ most representative points, _i.e._ a set of centroids, is obtained
from the original dataset. Although these centroids are the most
representative data points, they are not necessarily contained in the original
dataset and, consequently, a resampling process, allowing to choose points of
the original dataset (centrus), is required.

A study of the performance of the proposed dataset reduction method will be
done by training a logistic regression model with the original dataset and
comparing the results with those obtained when Kmeans and randomly
undersampled datasets are used.

###  6.1 KMeans Algorithm

Considering a clustering algorithm, Kmeans iteratively tries to separate data
into independent groups [38]. This separation is done using the Lloyd’s
algorithm [39], based on the minimal variability of samples within each
cluster. The KMeans algorithm requires the specification of the desired number
of clusters ($k$) a priori. The following steps were used:

  1. 1.

Initialization of the centroids: using the scikit-learn implementation, it is
possible to do it in two different ways, _random_ and _k-means++_ [40]:

     * •

Random: $k$ random samples of the original dataset are chosen.

     * •

K-means++: $k$ samples of the original dataset are chosen based on a
probabilistic approach, leading to the centroids being initialized far away
from each other.

Assuming there is enough time, the algorithm will always converge, although
the convergence to an absolute minimum is not guaranteed. The K-means++
initialization helps to address this issue. Furthermore, for both
initializations, the algorithm, by default, runs several times with different
centroid seeds, with the best result being the output.

  2. 2.

Assignment: Each data point $x_{i:}$ is addressed to a cluster
$c_{k^{\prime}}$, in such a way that the inertia is minimized:

$k^{\prime}=\mathop{argmin}_{k}\Biggl{\\{}\sum_{j=0}^{F-1}(x_{ij}-\mu_{kj})^{2}\Biggr{\\}},$
(5)

where $F$ is the dimensionality, _i.e._ the number of features, $\mu_{k}$ is
the centroid of the cluster $c_{k}$ and $j$ stands for the $(j+1)^{th}$
feature.

  3. 3.

Update of the centroids’ position: The new centroids are just the means
positions of each cluster, _i.e._

$\mu_{k:}=\frac{\sum_{i=0,x_{i:}\in c_{k}}^{n_{k}-1}x_{i:}}{n_{k}},$ (6)

with $n_{k}$ being the number of samples addressed to $c_{k}$. It should be
noted that if $n_{k^{\prime}}=0$ the centroid $\mu_{k^{\prime}}$ doesn’t
change.

  4. 4.

Iteration: Steps 2 and 3 are repeated until the maximum number of iterations
is reached or until the result converges, _i.e._ the centroids don’t change.

The KMeans algorithm was used separately for the signal and background
samples, with the corresponding weights being used.

###  6.2 Dataset Resampling

As previously mentioned, although centroids are the most representative
points, they are not necessarily contained in the original dataset. Hence, it
was chosen to consider 10 neighbors of each centroid to determine each
centrus, _i.e._ the 10 nearest points of the original dataset.

The position of each centrus was determined using the weighted mean of the
position of the neighbors,

$W=\frac{\sum_{i=0}^{9}x_{i:}\times w_{i}}{\sum_{i=0}^{9}w_{i}},$ (7)

where $W$ is the mean position, $x_{i:}$ is the $(i+1)^{th}$ nearest point and
$w_{i}$ the weight of the sample.

The sample weight of each centrus was calculated based on the number of
samples of the same label (_i.e._ signal or background) on the original
dataset:

$w_{i}=\frac{1}{n},$ (8)

with $w_{i}$ being the weight of the $(i+1)^{th}$ centrus and $n$ the number
of samples in the original dataset with the same label of this centrus.

##  7 Quantum and Classical Machine Learning Training

The training of the QML algorithms used in the current paper requires the use
of optimizers. Two different ones were considered: Adam [41] and tree-
structured Parzen estimator sampler (TPE) [42, 43].

The Adam optimizer uses an extension of stochastic gradient descent,
leveraging techniques such as adaptive moment estimation, being extensively
used in optimization problems, namely in the context of machine learning.
Nonetheless, since there is no reason to expect, _a priori_ , that it will
work equally well in the context of QML, where specific challenges are
expected, the TPE optimizer was also tested.

The TPE is a Bayesian optimization algorithm first developed for HP tuning in
the context of machine learning. In the current study, it will be used to
optimize VQC weights in a way very similar to what is typically done for HP
tuning. TPE is implemented using Optuna [44], a library focused on HP
optimization for machine learning models. TPE works by choosing a parameter
candidate that maximizes the likelihood ratio between a Gaussian Mixture Model
(GMM) fitted to the set of parameters associated with the best objective
values, with another GMM being fitted to the remaining parameter values. In
the context of HEP, TPE has also been used to explore parameter spaces of BSM
models [45].

Different machine learning methods were optimized, namely a LR, a SVM and a
VQC. The corresponding training was done for the set of HP summarized in Table
5, where the scanned values are also listed. For each set of HP, 5 models were
trained on 5 different subsets of the initial dataset (random sampling).

Variable HP | Possible values
---|---
Feature Selection | [PCA, SBS]
Number of Data points | [100, 500, 1k, 5k]
Number of Features | [1, 2, 3, 4, 5]
Number of Layers | [1, 2, 3, 4, 5]
Fixed HP | Fixed values
Max Epochs | 500
Batch Size | Size of the dataset
Learning Rate (LR) | 0.03
Table 5: List of scanned hyperparameters. The LR parameter is used only for
the VQC optimized by Adam while the number of layers is only used by the VQCs.

For both optimizers, the considered cost function used is the squared error,
with the individual Monte Carlo samples being properly weighted. During the
training of VQCs, the inference was done on the validation dataset at 5 epoch
intervals, the AUC computed and only the best-performing model, according to
the previously mentioned metric, was considered.

###  7.1 Adam Implementation Details

The training starts with the initialization of the weight vector. This is done
randomly with an order of magnitude of $10^{-2}$, which is followed by
training iterations until a maximum number of epochs is reached. At each
iteration, the model is inferred with the training dataset, the cost function
calculated and the model parameters updated via the Adam optimizer. A summary
of Adam-optimized VQC training is shown in Algorithm 2.

Algorithm 2 Adam Training

$params\leftarrow params\\_initialization()$

for $epoch\leftarrow 1$ to $max\\_epochs$ do

$loss\leftarrow cost(params)$

$params\leftarrow optimizer.step()$

if $epoch\\_number\%5==0$ or $epoch\\_number==max\\_epochs$ then

$validation\\_step()$

end if

end for

###  7.2 TPE Implementation Details

We use the Optuna implementation of the TPE sampler. Being a Bayesian
optimization algorithm, TPE works very differently to Adam, which is a
gradient descent algorithm. In TPE, for every training iteration, each
parameter is replaced by a new value acquired sampling from a Gaussian Mixture
Model of good points, which is then used to compute the loss function. At each
epoch, the algorithm computes new values for the model parameters. With the
value of the loss function of the suggested parameters, TPE will update its
internal Gaussian Mixture Models of good and bad points, which will allow it
to learn what are good suggestions as more parameter values are sampled. Since
TPE is a Bayesian algorithm, it does not need to compute derivatives of the
loss function, as Adam does, which might allow for a light workload when
running trainings on quantum computers.

###  7.3 Simulation Enviroment

Throughout this work, the PennyLane package [46] was used for hybrid quantum-
classical machine learning applications. Pennylane’s default.qubit quantum
simulator – a simple simulator for quantum circuits – was used to train and
evaluate the performance of the different QML algorithms. Later, the
performance of the trained algorithms will be measured on IBM’s quantum
computers using PennyLane integration with IBM quantum computing framework,
Qiskit [47].

##  8 Simulation Results

###  8.1 Feature Reduction

The results indicate that QML circuits trained with SBS data are generally
unstable and very susceptible to fluctuations in the randomly sampled data, as
can be seen in Figure 3. Specifically, it is evident that using PCA-originated
data produces significantly more stable results.

The performance of both optimizers, Adam and TPE, is usually saturated with
only 2 layers. This effect is most noticeable when the number of features is
greater or equal to 3. When considering only the PCA-obtained results, the two
optimizers are compatible for most of the configurations tested. Exceptions
occur when using a high number of features ($\geq 4$) and only one layer,
where TPE outperforms Adam, and when using a high number of features ($\geq
4$) and more than one layer, where the opposite happens and Adam outperforms
TPE.

The shallow ML methods trained on the same data as the VQCs are shown in
Figure 4. The AUC scores obtained in this case are more stable for both the
PCA and SBS datasets. The performance in both cases is saturated when using 2
features and the models trained with SBS data outperform the PCA-trained
models, contrary to what was observed for the QML case. It should also be
noted that the SVM outperforms LR in all cases except when only one feature is
used, which is not surprising since SVMs are more sophisticated classifiers.

For the best set of HP, VQCs trained using TPE and Adam have performed
similarly to the shallow ML methods (_c.f._ with Figure 5 and Figure 14,
respectively). It was also observed that there are no cases where QML
outperforms any of the shallow methods tested. The TPE optimizer regime
produced the best performance for QML, achieving an AUC score of $0.841\pm
0.051$, as shown in Figure 5.

Figure 3: Plot grid representing the results for both Adam and TPE-Trained
VQCs. Each data point represents the AUC score on the test dataset of a
different set of HP, as listed in Table 5. The error bar represents the
standard deviation associated with each data point since each point is the
average of 5 different random samplings from the data.

Figure 4: Plot grid representing the results for the considered shallow
methods. Each data point represents the AUC score on the test dataset of a
different set of HP, as listed in Table 5. The error bar represents the
standard deviation associated with each data point since each point is the
average of 5 different random samplings from the data.

Figure 5: ROC of the best HP set, using TPE’s QML model average AUC score as a
metric and the corresponding shallow methods ROCs for the same data. The HP
for this run are SBS for feature method, 100 data points, 1 feature and 5 VQC
layers.

The reduction algorithms studied come with an additional computational cost
compared to using the original dataset directly. In particular, the SBS
algorithm added an overhead of 1 hour for running the XGBoost algorithm and
selecting the features with more classification power. On the other hand, the
PCA algorithm took a sub-minute negligible time to complete. However, since
these algorithms only need to be run once, before the training, and given that
the grid search for the VQC, SVM and LR algorithms took over 200 CPU hours on
a dual-Intel(R) Xeon(R) Gold 6348 machine, in the end the computational cost
of the classical reduction algorithms is negligible.

####  8.1.1 VQC’s Robustness to Discrete Features

In the previous section it was noted that there was a significant variability
in the final score of QML models, especially when training with SBS data. In
fact, VQCs, being variational algorithms, are highly susceptible to small
fluctuations in the data, which can have a correspondingly significant impact
on the computed AUC. Additionally, numerical instabilities caused by
computational floating point accuracy were observed during the validation
step, leading to considerable fluctuations in the computed AUC in this regime.

To further investigate this behavior, which was not observed at the same level
on the PCA-trained circuits, we looked at the AUC distributions produced by
QML models as a function of the number of features. As shown in Figure 6, it
is clear that the instability in SBS results occurs when more than two
features are used.

Additionally, we produced visualisations of the decision regions of the models
trained using both feature selection methods. We focused on runs that used two
features, as this is where the problem originated. Figure 7 and Figure 8 show
the decision regions obtained with each model for one representative run,
illustrating the sensitivity of each boundary to variations in the data, for
SBS or PCA. The SBS features used are listed in Table 2, where the second
feature, the number of muons in the event, is a discrete variable.555In ML
literature this is called a categorical variable. However, we note that even
though it is categorical, it is still ordinal. As there are no non-ordinal or
non-binary categorical variables in our dataset, we will refer to these
variables as discrete instead of categorical for the rest of this work.

Figure 6: Distribution of the AUC values obtained for the QML model as a
function of the number of features used in training, evaluated on the test
datasest, for SBS and PCA inputs.

Figure 7: Decision regions of the 3 different architectures in a run where
large variability of results for the QML SBS-trained model was observed. This
case uses SBS data, Adam as an optimizer, 100 data points for training and 2
layers for the circuit.

Figure 8: Decision regions of the 3 different architectures in a run where
large variability of results for the QML PCA-trained model was observed. This
case uses PCA data, Adam as an optimizer, 100 data points for training and 1
layer for the circuit.

While LR and SVMs are robust in the presence of discrete variables, they may
pose a challenge for continuous learning algorithms such as VQCs. It is
therefore possible that the variability observed when using different sub-
samples of SBS data could be attributed to the use of this discrete variable.
To investigate this, we conducted the SBS feature selection once again, this
time excluding all discrete variables - yielding Table 6. The VQC circuits
where once again trained using this modified list of inputs in a limited study
of 2 features only, as illustrated in Figure 9.

Feature | AUC
---|---
$\cancel{E}_{T}$ | 0.817
large-R jet $\tau_{1}$ | 0.576
large-R jet $\tau_{3}$ | 0.316
Jet2 $p_{T}$ | 0.313
Jet1 $p_{T}$ | 0.292
Table 6: Features selected by the SBS Algorithm and their respective AUC Score
on the training dataset with all the discrete features removed.

Figure 9: Plot grid representing the results for both Adam and TPE-Trained
VQCs. Each data point represents the AUC score on the test dataset of a
different set of HP, as listed in Table 5, with the two features restriction.
The error bar represents the standard deviation associated with each data
point since each point is the average of 5 different random samplings from the
data.

Using the discrete-free SBS version to train the VQC led to significantly
better AUC scores, outperforming PCA-trained QML models with an average AUC
score of around $0.85$, although still with larger variability than that of
the PCA-trained VQCs. This is a notable departure from our previous
observations in subsection 8.1, where including discrete features in the SBS
feature selection methodology resulted in erratic performance with no instance
of outperforming PCA (except in cases where only one continuous feature was
used). Therefore, we found that excluding discrete variables during feature
selection led to better performance for VQC circuits in a limited study of 2
features, compared to when discrete variables were included. This indicates
that the choice of input features is crucial for achieving high accuracy in
quantum machine learning, and future studies should consider the impact of
discrete variables on VQC performance. The findings may inform future choices
in selecting input features for VQC circuits to optimize model performance.

###  8.2 Dataset Reduction

####  8.2.1 Implementation of KMeans

The performance of the KMeans algorithm was tested initially by training LR
models with 10 reduced datasets and selecting a different number of $k$
features ($k\in[1,2,3,4,5]$) obtained with the SBS algorithm. The KMeans
algorithm considers the sample weight and, in order to have an equal number of
signal and background centroids, it was separately applied to the signal and
background data. Since state-of-the-art quantum computing requires small
datasets, the data reduction studies were done for datasets with 100, 500,
1000, and 5000 data points and the number of features previously mentioned.

Two configurations were studied: the framework presented in subsection 6.1 was
applied to the training and test datasets; and only to the training datasets
(with test datasets obtained through random undersampling).666Throughout this
article random undersampling refers to the random selection of data points
from the original dataset. In the ML subfield of imbalanced learning, the
proper methodology is to use resampling algorithms only during training, but
not during validation or test. In this section we present results of these two
cases as a comparison, but later we will restrict to random undersampling
during validation and testing.

The mean AUC score and respective standard deviation found using KMeans for
train and test datasets are summarised in Figure 10. The results obtained
using the KMeans algorithm for the training dataset and random undersampling
in the test signal and background samples are presented in Figure 11. In order
to provide a benchmark point for comparison with the performance of the
reduced datasets, a LR model was trained on the full original dataset, with
results shown in both Figures.

Figure 10: Average AUC score and corresponding standard deviation, represented
as uncertainty bands, for different numbers of clusters as a function of the
number of features. The training and testing datasets were reduced using the
KMeans algorithm. In each case, 10 different reduced test datasets were used.

Figure 11: Average AUC score and corresponding standard deviation, represented
as uncertainty bands, for different numbers of clusters as a function of the
number of features. The training dataset was reduced using the KMeans
algorithm. In each case, 10 different randomly undersampled test datasets were
used.

It can be seen in Figure 10 and Figure 11 that using the KMeans algorithm to
reduce the training dataset results in AUC scores that are compatible with the
performance obtained using the full original dataset.

This study shows that although KMeans is a more sophisticated algorithm for
data reduction than random undersampling, in the HEP case under consideration
no significant deterioration of the performance is observed when using it,
suggesting that in this study the dataset composed of prototypes is a good
representative of the whole dataset in the small dataset regime, which is
explored in this work.

####  8.2.2 Application to QML

The QML, SVM, and LR models were trained using KMeans reduced datasets as well
as random undersampling, for different dataset sizes. In this comparison, the
HP for the VQC are the ones previously found to be the best, i.e. one feature
chosen with the SBS method and five VQC layers for the architecture. The
metric used to compare all models is the AUC score average of five different
runs.

For all cases, the test and validation sets were reduced using random
undersampling, hence, for each dataset size there are one train, five
validation and five test datasets. The choice to keep random sampling for the
test dataset, rather than KMeans reduction, is to ensure that our methodology
represents the test samples as close to the original dataset as possible,
ensuring that sophisticated resampling techniques do not significantly modify
the data.

The obtained results are shown in Figure 12. It can be seen that the
performance for the KMeans reduced dataset is compatible with the one obtained
using the dataset reduced through random undersampling, for QML and CML
models. Furthermore, the performance achieved by the simulated VQCs is
identical within the statistical uncertainties to the performances by the SVM
and LR, in agreement to what was observed in subsection 8.1.

Nonetheless, it should be emphasized that the model trained with random
undersampling needs to be trained several times for achieving these average
scores, as many times as the number of reduced datasets used. On the other
hand, the models using the KMeans reduced dataset need to be trained only
once. This can be relevant in the context of quantum computers, where access
is often subject to long queues and thus the number of accesses can be a
limiting factor. While the KMeans reduction technique brought an overall
increase in time of around 1%, this change is negligible taking into account
the reduction in number of accesses.

Figure 12: Comparison between the QML, SVM, and LR models when trained with
the TPE and the best set of HP for different dataset sizes for both random
undersampling (regular) and KMeans reduced datasets.

##  9 Real Quantum Computers Results

Until this point, only simulated quantum environments were used. In order to
test the performance in real quantum computers, and thus validate the
simulation results, the Pennylane framework was used as the integration layer
with Qiskit, which works effectively as a direct API to the quantum computers
provided by IBM.

In this study only the best performant model HP-set was used, _i.e._ the TPE-
trained VQC. This VQC was implemented and its test set was inferred on six
different quantum systems with identical architectures, all freely available.
Evaluating our model in multiple identical quantum systems allows us to get an
idea of the scale of the associated systematic uncertainty via the variability
of the observed results. Since the implemented circuits are small, no error
mitigation techniques were implemented. IBM’s transpiler optimization level
was set to 3777The level 3 of optimisation corresponds to the heaviest
optimisation inherently implemented. [47] and, for each event, the final
expectation value was computed by averaging $20k$ shots on the quantum
computer. The obtained results, shown in Figure 13, are compatible with the
simulated ones (Figure 5).

Figure 13: Final ROC curve of the best-performing model when inferred on the
test dataset in 6 different IBM systems. The average AUC scores and the
corresponding standard deviations are also shown. The colours in each subplot
stand for different runs of the same circuits in the same QC.

##  10 Conclusion

In this paper, we assessed the feasibility of applying QML to large datasets
usually found in HEP. To allow for a fair comparison between CML and QML
models, we opted to use shallow classical methods as opposed to deep methods,
which require large datasets that are not viable given the limitations of the
current quantum computers. Our results indicate that there is comparable
performance between CML and QML when tested on the same small dataset regime.

To achieve this, our study first compared feature selection techniques,
showing that while SBS can produce the best performant QML model, it generally
yielded worse and more unstable results than PCA. Additionally, we found this
was produced by using discrete variables in VQCs, highlighting the suitability
of PCA-transformed data for QML applications in the HEP context, where
discrete variables are commonly used to describe collider events.

Our grid search over different HP combinations of VQC ran in simulation
provided no evidence of quantum advantage in our study. We confirmed the
results by running the best performing configuration on real-world quantum
systems, obtaining compatible performances and therefore validating our
conclusions. We compared the performance of TPE and Adam optimizers in QML and
found that TPE achieves competitive results. Being a gradient-free optimiser,
TPE offers the advantage that it can lead to faster training with a smaller
memory usage when compared to Adam, which in principle can further facilitate
the application of QML in current quantum computers.

We then explored data reduction techniques, finding that reducing the dataset
size with the KMeans algorithm produces results that are similar to those
obtained from random undersampling. This finding is significant in that it
means that the model can achieve similar performance with fewer accesses to a
quantum computer during training, which is a considerable bottleneck in
current QML.

In conclusion, while our study found no evidence of quantum advantage in the
current state of QML within the context of large HEP datasets, the performance
of QML models was comparable to that of classical machine learning models when
restricted to small dataset regimes. Our findings suggest that using dataset
reduction techniques enables us to use large datasets more efficiently to
train VQCs, facilitating the usage of current quantum computers in large
datasets often found in HEP.

## Acknowledgments

We acknowledge the use of IBM Quantum services for this work. The views
expressed are those of the authors, and do not reflect the official policy or
position of IBM or the IBM Quantum team.

This work was supported by Fundação para a Ciência e a Tecnologia, Portugal,
through project CERN/FIS-COM/0004/2021 (“Exploring quantum machine learning as
a tool for present and future high energy physics colliders”).

We thank Declan Millar, Nuno Peres and Tiago Antão for the very useful
discussions and Ricardo Ribeiro for kindly providing access to some of the
computing systems used in this work. We also thank Henrique Carvalho for the
help in producing Figure 1.

## Code Availability

The code used in this work is publicly available in
https://github.com/mcpeixoto/QML-HEP.

## Appendix

In this appendix we show the results obtained for the Adam-trained VQCs,
namely the ROC curve for the best HP set for the QML model, shown in Figure
14. This result can be compared to the TPE-trained VQC results shown in Figure
5. As noted in subsection 8.1, the two optimizers are compatible.

Figure 14: ROC of the best HP set, using Adam’s QML model average AUC score as
a metric. The HP for this run are SBS for feature method, 100 data points, 1
feature, and 1 VQC layer. The corresponding shallow methods ROCs for the same
data have an AUC of $0.795\pm 0.105$ for SVM and $0.838\pm 0.053$ for LR.

## References

  * [1] J. Ellis,  _Outstanding questions: Physics beyond the Standard Model_ ,  Phil. Trans. Roy. Soc. Lond. A 370, 818 (2012).
  * [2] D. Guest, K. Cranmer and D. Whiteson,  _Deep learning and its application to LHC physics_ ,  Ann. Rev. Nucl. Part. Sci. 68, 161 (2018).
  * [3] M. Feickert and B. Nachman,  _A Living Review of Machine Learning for Particle Physics_ (2021),  2102.02770.
  * [4] I. Shapoval and P. Calafiura,  _Quantum associative memory in hep track pattern recognition_ ,  In _EPJ Web of Conferences_ , vol. 214, p. 01012. EDP Sciences (2019).
  * [5] C. Tüysüz, F. Carminati, B. Demirköz, D. Dobos, F. Fracas, K. Novotny, K. Potamianos, S. Vallecorsa and J.-R. Vlimant,  _Particle track reconstruction with quantum algorithms_ ,  In _EPJ Web of Conferences_ , vol. 245, p. 09013. EDP Sciences (2020).
  * [6] F. Bapst, W. Bhimji, P. Calafiura, H. Gray, W. Lavrijsen, L. Linder and A. Smith,  _A pattern recognition algorithm for quantum annealers_ ,  Computing and Software for Big Science 4(1), 1 (2020).
  * [7] L. Funcke _et al._ ,  _Studying quantum algorithms for particle track reconstruction in the LUXE experiment_ (2022),  2202.06874.
  * [8] A. Zlokapa, A. Anand, J.-R. Vlimant, J. M. Duarte, J. Job, D. Lidar and M. Spiropulu,  _Charged particle tracking with quantum annealing optimization_ ,  Quantum Machine Intelligence 3(2), 1 (2021).
  * [9] S. Das, A. J. Wildridge, S. B. Vaidya and A. Jung,  _Track clustering with a quantum annealer for primary vertex reconstruction at hadron colliders_ (2019),  1903.08879.
  * [10] A. Y. Wei, P. Naik, A. W. Harrow and J. Thaler,  _Quantum algorithms for jet clustering_ ,  Physical Review D 101(9), 094015 (2020).
  * [11] A. Mott, J. Job, J.-R. Vlimant, D. Lidar and M. Spiropulu,  _Solving a higgs optimization problem with quantum annealing for machine learning_ ,  Nature 550(7676), 375 (2017).
  * [12] A. Zlokapa, A. Mott, J. Job, J.-R. Vlimant, D. Lidar and M. Spiropulu,  _Quantum adiabatic machine learning with zooming_ ,  Bulletin of the American Physical Society 66 (2021).
  * [13] K. Terashi, M. Kaneda, T. Kishimoto, M. Saito, R. Sawada and J. Tanaka,  _Event classification with quantum machine learning in high-energy physics_ ,  Computing and Software for Big Science 5(1), 1 (2021).
  * [14] S. L. Wu, S. Sun, W. Guan, C. Zhou, J. Chan, C. L. Cheng, T. Pham, Y. Qian, A. Z. Wang, R. Zhang _et al._ ,  _Application of quantum machine learning using the quantum kernel algorithm on high energy physics analysis at the LHC_ ,  Physical Review Research 3(3), 033221 (2021).
  * [15] A. Blance and M. Spannowsky,  _Quantum machine learning for particle physics using a variational quantum classifier_ ,  Journal of High Energy Physics 2021(2), 1 (2021).
  * [16] V. S. Ngairangbam, M. Spannowsky and M. Takeuchi,  _Anomaly detection in high-energy physics using a quantum autoencoder_ ,  Physical Review D 105(9), 095004 (2022).
  * [17] A. Gianelle, P. Koppenburg, D. Lucchesi, D. Nicotra, E. Rodrigues, L. Sestini, J. de Vries and D. Zuliani,  _Quantum Machine Learning for b-jet charge identification_ ,  JHEP 08, 014 (2022),  doi:10.1007/JHEP08(2022)014,  2202.13943.
  * [18] S. Y.-C. Chen, T.-C. Wei, C. Zhang, H. Yu and S. Yoo,  _Quantum convolutional neural networks for high energy physics data analysis_ ,  Physical Review Research 4(1), 013231 (2022).
  * [19] V. Belis, S. González-Castillo, C. Reissel, S. Vallecorsa, E. F. Combarro, G. Dissertori and F. Reiter,  _Higgs analysis with quantum classifiers_ ,  In _EPJ Web of Conferences_ , vol. 251, p. 03070. EDP Sciences (2021).
  * [20] J. Preskill,  _Quantum computing in the NISQ era and beyond_ ,  Quantum 2, 79 (2018),  doi:10.22331/q-2018-08-06-79.
  * [21] G. Li, Y. Ding and Y. Xie,  _Tackling the qubit mapping problem for nisq-era quantum devices_ ,  doi:10.48550/ARXIV.1809.02573 (2018).
  * [22] S. Gupta and R. Zia,  _Quantum neural networks_ ,  Journal of Computer and System Sciences 63(3), 355 (2001).
  * [23] M. Schuld, I. Sinayskiy and F. Petruccione,  _Simulating a perceptron on a quantum computer_ ,  Physics Letters A 379(7), 660 (2015).
  * [24] E. Farhi and H. Neven,  _Classification with quantum neural networks on near term processors_ (2018),  1802.06002.
  * [25] M. Schuld, A. Bocharov, K. M. Svore and N. Wiebe,  _Circuit-centric quantum classifiers_ ,  Physical Review A 101(3) (2020),  doi:10.1103/physreva.101.032308.
  * [26] R. LaRose and B. Coyle,  _Robust data encodings for quantum classifiers_ ,  Physical Review A 102(3), 032420 (2020).
  * [27] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos _et al._ ,  _Scikit-learn: Machine learning in Python_ ,  Journal of Machine Learning Research 12, 2825 (2011).
  * [28] M. Crispim Romao, N. F. Castro and R. Pedro,  _Simulated pp collisions at 13 TeV with 2 leptons + 1 b jet final state and selected benchmark Beyond the Standard Model signals_ ,  Dataset on Zenodo,  doi:10.5281/zenodo.5126747 (2021).
  * [29] M. Crispim Romão, N. F. Castro and R. Pedro,  _Finding New Physics without learning about it: Anomaly Detection as a tool for Searches at Colliders_ ,  Eur. Phys. J. C 81(1), 27 (2021),  doi:10.1140/epjc/s10052-021-09813-2,  [Erratum: Eur.Phys.J.C 81, 1020 (2021)],  2006.05432.
  * [30] G. Durieux, F. Maltoni and C. Zhang,  _Global approach to top-quark flavor-changing interactions_ ,  Phys. Rev. D 91(7), 074017 (2015),  doi:10.1103/PhysRevD.91.074017,  1412.7166.
  * [31] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro,  _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ ,  JHEP 07, 079 (2014),  doi:10.1007/JHEP07(2014)079,  1405.0301.
  * [32] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen and P. Z. Skands,  _An introduction to PYTHIA 8.2_ ,  Comput. Phys. Commun. 191, 159 (2015),  doi:10.1016/j.cpc.2015.01.024,  1410.3012.
  * [33] M. Selvaggi,  _DELPHES 3: A modular framework for fast-simulation of generic collider experiments_ ,  J. Phys. Conf. Ser. 523, 012033 (2014),  doi:10.1088/1742-6596/523/1/012033.
  * [34] M. Cacciari, G. P. Salam and G. Soyez,  _The anti- $k_{t}$ jet clustering algorithm_,  JHEP 04, 063 (2008),  doi:10.1088/1126-6708/2008/04/063,  0802.1189.
  * [35] J. Thaler and K. Van Tilburg,  _Identifying boosted objects with n-subjettiness_ ,  Journal of High Energy Physics 2011(3), 1 (2011).
  * [36] T. Chen and C. Guestrin,  _XGBoost: A scalable tree boosting system_ ,  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pp. 785–794. ACM, New York, NY, USA,  ISBN 978-1-4503-4232-2,  doi:10.1145/2939672.2939785 (2016).
  * [37] S. Raschka,  _Mlxtend: Providing machine learning and data science utilities and extensions to python’s scientific computing stack_ ,  The Journal of Open Source Software 3(24) (2018),  doi:10.21105/joss.00638.
  * [38] D. MacKay,  _Information Theory, Inference and Learning Algorithms_ ,  Cambridge University Press (2003).
  * [39] G. A. Wilkin and H. Xiuzhen,  _A practical comparison of two k-means clustering algorithms_ ,  BMC bioinformatics 9 (2008),  doi:10.1186/1471-2105-9-S6-S19.
  * [40] A. Vouros, S. Langdell, M. Croucher and E. Vasilaki,  _An empirical comparison between stochastic and deterministic centroid initialisation for k-means variations_ ,  Machine Learning 110, 1975 (2021),  doi:10.1007/s10994-021-06021-7.
  * [41] D. P. Kingma and J. Ba,  _Adam: A Method for Stochastic Optimization_ (2014),  1412.6980.
  * [42] J. Bergstra, R. Bardenet, Y. Bengio and B. Kégl,  _Algorithms for hyper-parameter optimization_ ,  Advances in neural information processing systems 24 (2011).
  * [43] J. Bergstra, D. Yamins and D. Cox,  _Making a science of model search: Hyperparameter optimization in hundreds of dimensions for vision architectures_ ,  In _International conference on machine learning_ , pp. 115–123. PMLR (2013).
  * [44] T. Akiba, S. Sano, T. Yanase, T. Ohta and M. Koyama,  _Optuna: A next-generation hyperparameter optimization framework_ ,  In _Proceedings of the 25rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ (2019).
  * [45] F. A. de Souza, M. Crispim Romão, N. F. Castro, M. Nikjoo and W. Porod,  _Exploring Parameter Spaces with Artificial Intelligence and Machine Learning Black-Box Optimisation Algorithms_ (2022),  2206.09223.
  * [46] V. Bergholm _et al._ ,  _PennyLane: Automatic differentiation of hybrid quantum-classical computations_ (2018),  1811.04968.
  * [47] M. S. ANIS, Abby-Mitchell, H. Abraham, AduOffei, R. Agarwal, G. Agliardi, M. Aharoni, I. Y. Akhalwaya, G. Aleksandrowicz, T. Alexander, M. Amy, S. Anagolum _et al._ ,  _Qiskit: An open-source framework for quantum computing_ ,  doi:10.5281/zenodo.2573505 (2021).