# Weak integrability breaking and level spacing distribution

D. Szász-Schagrin1, B. Pozsgay2 and G. Takács3
1BME Momentum Statistical Field Theory Research Group,
Department of Theoretical Physics,
Budapest University of Technology and Economics
2MTA-ELTE “Momentum” Integrable Quantum Dynamics Research Group (ELKH),
Department of Theoretical Physics, Eötvös Loránd University, Budapest
3MTA-BME Quantum Correlations Group (ELKH),
Department of Theoretical Physics,
Budapest University of Technology and Economics

(20th March 2021)

###### Abstract

Recently it was suggested that certain perturbations of integrable spin chains
lead to a weak breaking of integrability in the sense that integrability is
preserved at the first order of the coupling. Here we examine this claim using
level spacing distribution. We find that the volume dependent cross-over
between integrable and chaotic statistics is markedly different between weak
and strong breaking of integrability, supporting the claim that perturbations
by generalised currents only break integrability at higher order. In addition,
for the massless case we find that the critical coupling as a function of the
volume $L$ scales with a $1/L^{2}$ law for weak breaking as opposed to the
previously found $1/L^{3}$ law for the strong case.

##  1 Introduction

Mean values of generalised currents in one-dimensional integrable models have
attracted considerable interest lately, mainly due to the recent theory of
Generalised Hydrodynamics which describes non-equilibrium dynamics at the
Euler scale [1, 2]. These currents express the continuity relation for
conserved charges responsible for integrability, which can be exploited for a
hydrodynamic description of the ballistic flow of the quasi-particles. This
description is useful only as far as the exact mean values of the currents in
local equilibrium conditions are known, and indeed in the thermodynamic limit
a simple exact formula was postulated in [1, 2], which was proven for
relativistic quantum field theories [1, 3], and for the spin current of the
XXZ spin chain in [4], and the classical Toda chain [5]. These exact results
were extended to finite volume in [6] with a proof that applies to systems
including Heisenberg spin chains, the Lieb-Liniger interacting Bose gas, and
integrable quantum field theories with diagonal scattering. In [7] an
algebraic construction was given for the current operators of the integrable
spin chains, which led to an alternative rigorous proof of their mean values.

Recently [8] it was discovered that these exact results are connected to the
so-called long range deformations of integrable spin chains which emerged in
the context of the AdS/CFT correspondence [9, 10, 11]. These long range spin
chains are obtained as a one-parameter deformation of ordinary short-range
spin chains (such as the XXZ model) and preserve integrability to each finite
order in the deformation parameter $g$, with the interaction range growing
order-by-order. Strictly speaking, these deformations are only defined in
infinite volume. In the work [12] it was demonstrated that these long range
deformations have a deep connection to $T\bar{T}$-deformations of integrable
quantum field theory [13, 14] (see also [15, 16]), by sharing the same
algebraic origin which allowed the proof of factorisation for the expectation
values of the operators which trigger the deformation of the spin chain.

It turns out that the existence of the exact formulae for the current
expectation values is implied by the following observation made in [8]: for
each generalised current there is a long range deformation such that the given
current operator itself is the leading perturbing operator. This implies that
the perturbation of the spin chain by the generalised current operator is
integrable to the leading order in the deformation parameter $g$, but
integrability is generally expected to be broken at higher orders since to
maintain it necessitates the inclusion of progressively longer and longer
range interaction terms at higher orders of $g$. We call this breaking of
integrability ’weak’ in distinction to ’strong’ breaking of integrability
which already happens at the first order of perturbation, and cannot be
maintained by improving the perturbing operator order-by-order in $g$.

Integrability breaking was also recently investigated in [17] within the
framework of Generalised Hydrodynamics, where it was also found that the
perturbations by the current operators do not break integrability at the
leading order: while a generic perturbation with coupling $\lambda$ is
expected to thermalise the system on a time scale of $T\sim\lambda^{-2}$, this
effect is missing for the current operators on the same time scale.

The onset of non-integrable a.k.a. chaotic behaviour in many-body systems can
be investigated via the statistics of energy spectra [18, 19, 20, 21, 22, 23,
24, 25, 26]; in some cases even analytical results can be obtained [27]. For
integrable systems the level spacing is known to follow Poissonian statistics,
while for non-integrable case random matrix theory predicts the so-called
Wigner-Dyson distribution whose precise form depends on the relevant random
matrix ensemble.

In this short paper we set out to investigate the weak breaking of
integrability by exploring the level spacing statistics for the case of the
spin-1/2 XXZ spin chain. We compare a known strong integrability breaking term
which introduces nearest neighbour interactions to the weak integrability
breaking perturbation provided by the lowest non-trivial generalised current.
We start by presenting the Hamiltonian and the perturbations in Section 2, and
then turn to the results for the level spacing statistics in Section 3, where
we examine the cross-over between the Poissonian distribution characteristic
for integrable systems to the Wigner-Dyson distribution signalling the
breaking of integrability, as a function of the integrability breaking
coupling and the volume. We present our conclusions in Section 4.

##  2 The Hamiltonian and its perturbations

###  2.1 The XXZ spin chain and the current operator

The spin-1/2 XXZ spin chain of length $L$ is defined by the Hamiltonian

$H_{XXZ}=\sum_{i=1}^{L}\left[s_{i}^{x}s_{i+1}^{x}+s_{i}^{y}s_{i+1}^{y}+\Delta
s_{i}^{z}s_{i+1}^{z}\right]$ (1)

where $s^{j}=\frac{1}{2}\sigma^{j}$ are the spin operators with $\sigma^{j}$
denoting the Pauli matrices

$\sigma^{x}=\begin{bmatrix}0&1\\\
1&0\end{bmatrix}\quad\sigma^{y}=\begin{bmatrix}0&-i\\\
i&0\end{bmatrix}\quad\sigma^{z}=\begin{bmatrix}1&0\\\ 0&-1\end{bmatrix}$ (2)

and periodic boundary conditions $s_{L+1}^{a}\equiv s_{1}^{a}$. The model has
three phases controlled by the anisotropy parameter $\Delta$: for
$-1<\Delta<1$ the spectrum is gapless, while for $\Delta>1$ and $\Delta<-1$
there is a non-zero gap. The boundary points $\Delta=\pm 1$ are special as the
$U(1)$ symmetry of the theory generated by the conserved charge

$S^{z}=\sum_{i=1}^{L}s_{i}^{z}$

is enhanced to $SU(2)$.

The XXZ chain is integrable: there exists a family of local conserved
quantities

$Q_{n}=\sum_{l=1}^{L}q_{n,l}$

where $n=1,\dots,L$ and the charge density is supported on exactly $n$
neighbouring sites. Their conservation can be expressed as a continuity
equation

$j_{n,l+1}-j_{n,l}=i\left[q_{n,l},H_{XXZ}\right]$ (3)

where $j_{n,l}$ is the conserved current corresponding to the charge $Q_{n}$.

The first two charges are $Q_{1}=S^{z}$, and $Q_{2}$ which is nothing else but
the Hamiltonian (1). These charges are related to the $U(1)$ symmetry and time
translations, and are generally present in all the models considered in this
work, including the non-integrable ones.

The first nontrivial charge corresponding to integrability of the XXZ
Hamiltonian is

$\displaystyle Q_{3}$ $\displaystyle=\sum_{l=1}^{L}q_{3,l}$ (4) $\displaystyle
q_{3,l}=s_{l-1}^{x}s_{l}^{z}s_{l+1}^{y}-s_{l-1}^{y}s_{l}^{z}s_{l+1}^{x}+\Delta\left(-s_{l-1}^{z}s_{l}^{x}s_{l+1}^{y}+s_{l-1}^{z}s_{l}^{y}s_{l+1}^{x}-s_{l-1}^{x}s_{l}^{y}s_{l+1}^{z}+s_{l-1}^{y}s_{l}^{x}s_{l+1}^{z}\right)$

with the corresponding current given by

$\displaystyle j_{3,l}=-\frac{1}{2}\bigg{[}$ $\displaystyle
2\Delta\big{(}s_{l-2}^{x}s_{l-1}^{y}s_{l}^{x}s_{l+1}^{y}+s_{l-2}^{x}s_{l-1}^{z}s_{l}^{x}s_{l+1}^{z}+s_{l-2}^{y}s_{l-1}^{x}s_{l}^{y}s_{l+1}^{x}+s_{l-2}^{y}s_{l-1}^{z}s_{l}^{y}s_{l+1}^{z}$
$\displaystyle+s_{l-2}^{z}s_{l-1}^{x}s_{l}^{z}s_{l+1}^{x}+s_{l-2}^{z}s_{l-1}^{y}s_{l}^{z}s_{l+1}^{y}-s_{l-2}^{x}s_{l-1}^{y}s_{l}^{y}s_{l+1}^{x}-s_{l-2}^{y}s_{l-1}^{x}s_{l}^{x}s_{l+1}^{y}\big{)}$
$\displaystyle-2\left(s_{l-2}^{x}s_{l-1}^{z}s_{l}^{z}s_{l+1}^{y}+s_{l-2}^{y}s_{l-1}^{z}s_{l}^{z}s_{l+1}^{y}\right)-2\Delta^{2}\left(s_{l-2}^{z}s_{l-1}^{x}s_{l}^{x}s_{l+1}^{z}+s_{l-2}^{z}s_{l-1}^{y}s_{l}^{y}s_{l+1}^{z}\right)$
$\displaystyle-\frac{1+\Delta^{2}}{4}\left(s_{l-1}^{x}s_{l}^{x}+s_{l-1}^{y}s_{l}^{y}\right)-\frac{\Delta}{2}s_{l-1}^{z}s_{l}^{z}\bigg{]}\,.$

###  2.2 Perturbations, norms and effective coupling

We are interested in perturbing the XXZ chain by the current

$H_{J}=H_{XXZ}+g_{3}J\,,$ (5)

where

$J=\sum_{l=1}^{L}j_{3,l}$ (6)

which is supposed to break integrability only at higher order. As a benchmark
to the strength of integrability breaking, we also consider a perturbation
breaking integrability by the next-to-nearest-neighbor interaction (NNNI) term

$\mathcal{O}_{NNNI}=\sum_{i=1}^{L}s_{i}^{z}s_{i+2}^{z}\,,$ (7)

which leads to the Hamiltonian

$H_{NNNI}=H_{XXZ}+g_{N}\sum_{i=1}^{L}s_{i}^{z}s_{i+2}^{z}=\sum_{i=1}^{L}\left[s_{i}^{x}s_{i+1}^{x}+s_{i}^{y}s_{i+1}^{y}+\Delta
s_{i}^{z}s_{i+1}^{z}+g_{N}s_{i}^{z}s_{i+2}^{z}\right]\,.$ (8)

To compare the strength of integrability breaking, we introduce the effective
coupling $g_{\text{eff}}$ of an operator $\mathcal{O}$ as

$g_{\text{eff}}=g_{\mathcal{O}}n_{\mathcal{O}}$ (9)

where $g_{\mathcal{O}}$ is the coupling appearing in the Hamiltonian
$H=H_{0}+g_{\mathcal{O}}\mathcal{O}$, and $n$ is defined from the norm of the
operator as the coefficient of the leading asymptotics

$\left\|\mathcal{O}\right\|_{2}=\sqrt{\text{Tr
}\mathcal{O}^{\dagger}\mathcal{O}}$ (10)

as the coefficient $n$ of the leading asymptotic term

$||\mathcal{O}||_{2}=n_{\mathcal{O}}L+...\>,$ (11)

which is linear in the volume $L$ due to the expression of $\mathcal{O}$ as
the translation invariant sum of localised terms (6,7). The norm of $J$
depends on $\Delta$, while that of $\mathcal{O}_{NNNI}$ does not; some
explicit values are shown in Table 1. As illustrated in Figure 1, these norms
change linearly with the volume, apart from some fluctuations for $J$ which
are due to the fact that $L$ is changed in steps of $2$, while the one-site
term $j_{3,l}$ is localised on $4$ sites. The coefficient $n_{\mathcal{O}}$
can be computed by fitting a linear function to the norm values as a function
of $L$, and is shown in the last column of Table 1.

| $L$ | $10$ | $12$ | $14$ | $16$ | $18$ | $20$ |  $n_{\mathcal{O}}$
---|---|---|---|---|---|---|---|---
$\left\|J\right\|_{2}$ | $\Delta=0.2$ | $0.742$ | $1.306$ | $1.232$ | $1.841$ | $1.772$ | $2.346$ |  $0.143$
$\Delta=-1.2$ | $2.376$ | $3.115$ | $3.357$ | $4.614$ | $4.726$ | $6.023$ |  $0.348$
$\Delta=1.2$ | $2.274$ | $3.115$ | $3.455$ | $4.614$ | $4.412$ | $6.023$ |  $0.340$
$\left\|\mathcal{O}_{NNNI}\right\|_{2}$ |  | $1.500$ | $2.000$ | $2.500$ | $3.000$ | $3.500$ | $4.000$ |  $0.250$
Table 1:  The norms of the perturbing operators $J$ and $\mathcal{O}_{NNNI}$
Figure 1:  Norm of $J$ and $\mathcal{O}_{NNNI}$ as a function of the system
size in the three phases.

##  3 Level spacing distribution in the XXZ chain and its perturbations

The level spacings of a given system described by a Hamiltonian $H$ with
eigenvalues $\\{\lambda_{i}\\}$ are defined as the spacings $s_{i}$ between
ordered eigenvalues:

$\left\\{s_{i}=|\lambda_{i+1}-\lambda_{i}|\bigg{|}\leavevmode\nobreak\
\lambda_{i}\leavevmode\nobreak\ \text{eigenvalue}\right\\}$ (12)

The distribution of the normalised level spacing is called the level spacing
distribution $P(s)$:

$P(s)=\left\\{s_{i}/\bar{s}|\leavevmode\nobreak\ s_{i}\leavevmode\nobreak\
\text{spacing}\right\\}$ (13)

where $\bar{s}$ is the mean spacing. For integrable systems, the level spacing
statistics is exponential:

$P(s)_{I}=e^{-s}\,,$ (14)

while for non-integrable systems it is given by the Wigner-Dyson distribution,
which takes the following form for the orthogonal Gaussian ensemble111Due to
the fact that the Hamiltonians considered here are real and symmetric, it is
the orthogonal Gaussian ensemble which is relevant here.

$P(s)_{NI}=\frac{\pi}{2}se^{-\frac{\pi}{4}s^{2}}.$ (15)

Therefore, the level spacing distribution is an explicit indicator of
integrability and its breaking. The above predictions follow from random
matrix theory. For an integrable Hamiltonian, the different levels do not
interact due to the presence of higher conserved charges, and so the
distribution of eigenvalues is a Poissonian one, leading to exponential
distribution of the normalised level spacings. Breaking integrability results
in level repulsion, and so small level spacings are suppressed. When
considering the spectrum as a function of a parameter such as volume,
integrability implies that levels approaching as a function of volume each
other simply cross, while in the non-integrable case they avoid each other due
to repulsion. In the limit of infinite matrix size, the level spacing
distribution changes suddenly from exponential to Wigner-Dyson for any non-
zero value of the integrability breaking coupling $g$.

For a spin chain of finite length, however, the Hilbert space is finite
dimensional, and so the level spacing distribution is a continuous function of
the coupling, with the transition becoming sharper for larger volumes [18,
22]. In addition, when constructing the level spacing from the full spectrum
it is found to deviate from the random matrix prediction due to the structure
dictated by quasi-particle excitations, which is due to the locality of the
Hamiltonian. This can be solved by constructing the level spacing distribution
from the middle part of the spectrum, for which we take the middle two-thirds
of the computed levels. Furthermore, to get rid of degeneracies corresponding
to trivial symmetries [28, 29] the level spacing distribution is extracted
from a sector with total momentum zero, even spatial parity and a fixed
$S_{z}$ value. In the examples below we present the results obtained for
$S_{z}=2$; similar results were obtained for $S_{z}=1$ and $3$.

###  3.1 XXZ chain

The level spacing statistics of the XXZ chain of length $L=22$ in the
different phases is shown in figure 2.

Figure 2: Level spacing statistics of the XXZ spin chain of length $L=22$ for
(left to right) $\Delta=0.2,1.6,-1.6$ with the exponential distribution fitted
(dashed red line).

###  3.2 Perturbed system

Figure 3: A typical level spacing distribution of a non-integrable system with
the Wigner-Dyson distribution fitted (red dashed line). Left belongs to
$H_{NNNI}$, $L=22$ $\Delta=0.2$, $g_{\text{eff}}=0.1$ and right belongs to
$H_{J}$, $L=22$ $\Delta=0.2$, $g_{\text{eff}}=0.42$. Figure 4: Level spacing
statistics of $H_{NNNI}$ for $\Delta=0.2$. Figure 5: Level spacing statistics
of $H_{J}$ for $\Delta=0.2$.

Results for the NNNI perturbation in the gapless phase can be seen in Figure
4. Similar results were found in the other two phases. As expected, the change
in statistics occurs faster for longer chains. The same can be said for the
$J$ perturbation which is presented in Figure 5. It can be also seen that the
current perturbation breaks integrability for much higher $g_{\text{eff}}$
than the next nearest neighbour interaction which supports the idea that it
only breaks integrability at higher orders.

###  3.3 Integrability breaking and perturbation theory

To investigate the order of integrability breaking it is tempting to try and
construct the level spacing distribution for the spectrum constructed from
perturbation theory in the coupling $g$. However it turns out that this is not
possible. The level spacing distribution computed for the spectrum in first
order perturbation theory remans exponential even for the NNNI perturbation
(7), as shown in Fig. 6.

Figure 6: Level spacing statistics of $H_{NNNI}$ as computed from the spectrum
obtained by first order perturbation theory for the case $L=22$, $\Delta=0.2$.

To understand the reason, consider how a level crossing is lifted by a
perturbation. In our model, the generic level crossing obtained as a function
of the volume $L$ (continued to real values) happens between two levels. Let
us consider a generic perturbation of the corresponding two level subsystem:

$\begin{bmatrix}a_{1}(L)+g\Delta a_{1}(L)&g\epsilon(L)\\\
g\epsilon(L)&a_{2}(L)-g\Delta a_{2}(L)\end{bmatrix}\,,$ (16)

For the unperturbed Hamiltonian $g=0$, the level crossing is located at the
volume $L_{0}$ given by $a_{1}(L_{0})=a_{2}(L_{0})$; all that happens at first
order is that the level crossing is shifted to a location given to the
location $L_{*}$ which satisfies $a_{1}(L_{*})+g\Delta
a_{1}(L_{*})=a_{1}(L_{*})-g\Delta a_{2}(L_{*})$, where $L_{*}$ can be computed
as a series in $g$:

$L_{*}(g)=L_{0}-g\frac{\Delta a_{1}(L_{0})-\Delta
a_{1}(L_{0})}{a_{1}^{\prime}(L_{0})-a_{2}^{\prime}(L_{0})}+O\left(g^{2}\right)\,.$

Therefore first order perturbation theory does not introduce a repulsion
between levels, so the exponential distribution is unchanged.

When carrying the perturbative calculation to second order, the resulting
energy levels turn out to be numerically unstable and no meaningful statistics
can be constructed. For the two-level system above, an exact calculation of
the energy levels shows that the off-diagonal term leads to level repulsion,
which is responsible for changing the level spacing distribution. However, a
perturbative evaluation of the two nearby energy levels results in

$E_{1,2}=a_{1,2}(L)+g\Delta a_{1,2}(L)\pm
g^{2}\frac{2\epsilon(L)^{2}}{a_{1}(L)-a_{2}(L)}+O\left(g^{3}\right)$ (17)

which is unstable in the vicinity of a level crossing due to the presence of
the energy difference denominator. As a result, the order of integrability
breaking cannot be deduced by considering level spacing distribution of the
spectrum obtained in perturbation theory in the coupling $g$.

###  3.4 Crossover coupling and its volume dependence

Let’s define the crossover coupling $g_{\text{cr}}$ as the value of the
effective coupling $g_{\text{eff}}$ for which the maximum of the intermediate
distribution is at $x_{0}/2$, $x_{0}=\sqrt{2/\pi}$ being the position of the
maximum of the exact Wigner-Dyson distribution. This can be obtained for all
$L$ and we expect $g_{\text{cr}}(L)$ to be a monotonically decreasing
function. More precisely, in the gapless phase it is expected to have a power-
like behaviour $g_{\text{cr}}\propto L^{-\alpha}$ [30], while in the gapped
phase finite volume corrections are expected to show exponential decay in the
volume.

####  3.4.1 Gapless phase

In the following we give results obtained by carrying out the above described
method for the NNNI and current perturbations in the gapless phase, namely for
$\Delta=0.2$. Obtaining the maximum position of the intermediate distribution
was obtained by applying a Gaussian filter of kernel radius $r$ to the raw
statistics. Determination of the position of the maximum of the level spacing
distribution is illustrated for the NNNI perturbation in figure 8, and for
current perturbation in figure 7 for a chain of length $L=22$. Even from just
a quick look at these figures the difference in the ’strength’ of
integrability breaking is already clear.

Figure 7:  Obtaining the maximum positions from the smoothed statistics for a
typical dataset. This particular case belongs to the $J$ perturbation in
$L=22$, $\Delta=0.2$. Figure 8:  Obtaining the maximum positions from the
smoothed statistics for a typical dataset. This particular case belongs to the
NNNI perturbation in $L=22$, $\Delta=0.2$.

The obtained maximum positions can be seen in figure 9 along with the parabola
$f(g)$ fitted to the resulting points around $x_{0}$. The crossover scale is
then obtained by solving $f(g_{cr})=x_{0}$.

Figure 9:  Obtaining $g_{cr}$ from the maximum positions for different
effective couplings for $J$ (left) and NNNI perturbation (right). The red
solid line is the parabola fitted to the red markers in the vicinity of
$x_{0}$ (indicated by orange dashed line).

The $g_{cr}(L)$ functions obtained this way for both the NNNI and current
perturbations can be seen in figure 10, and decay with a power of the volume
$L$. Fitting a linear function $\log g_{cr}(L)=a+b\log L$ results in the
exponents

$\displaystyle b_{J}=-1.99\pm 0.18$ (18) $\displaystyle b_{NNNI}=-3.11\pm
0.27$ (19)

The value obtained for the NNNI case is consistent with the universal exponent
$-3$ claimed for integrability breaking in a gapless chain [30], while the one
obtained for the current perturbation is in agreement with the conjecture that
the current perturbation breaks integrability at higher orders in perturbation
theory.

Figure 10:  The crossover scale as a function of the volume for the current
(blue markers) and NNNI perturbations (red markers) on normal (left) and log-
log scale (right). Dashed lines correspond to the fitted $AL^{-3}$ (red) and
$BL^{-2}$ (blue) respectively.

####  3.4.2 Gapped phase

In the gapped phase, finite size scaling is expected to be exponential in
volume. The crossover coupling as a function of the volume is shown for
$\Delta=1.6$ in Figure 11, and for $\Delta=-1.6$ in Figure 12. For the
massivel case, $\log g_{cr}(L)$ can be fitted with a linear function $c+dL$,
with the following values for the coefficients $d$ for $\Delta=1.6$:

$\displaystyle d_{J}=-0.056\pm 0.014$ $\displaystyle d_{NNNI}=-0.157\pm 0.020$

and for $\Delta=-1.6$:

$\displaystyle d_{J}=-0.063\pm 0.014$ $\displaystyle d_{NNNI}=-0.137\pm 0.043$

Again we see a marked difference between the strength of integrability
breaking for the two perturbations. The critical coupling for the current
perturbation in any given volume is an order of magnitude larger than for the
NNNI perturbation, and the coefficient $d$ describing its decay with the
volume is also significantly smaller, again in agreement with the conjecture
that the current perturbation breaks integrability at higher orders in
perturbation theory.

Figure 11:  The crossover scale as a function of the volume for the current
(blue markers) and NNNI perturbations (red markers) for $\Delta=1.6$ (gapped
phase) on normal (left) and log scale (right). Dashed lines correspond to the
fitted exponentials.

Figure 12: The crossover scale as a function of the volume for the current
(blue markers) and NNNI perturbations (red markers) for $\Delta=-1.6$ (gapped
phase) on normal (left) and log scale (right). Dashed lines correspond to the
fitted exponentials.

##  4 Conclusions

In this work we examined the crossover from integrable to chaotic behaviour
for weak breaking of integrability, defined as a perturbation of an integrable
system which preserves integrability to first order in perturbation theory. As
an example we took the XXZ spin chain perturbed by one of the currents that
appear in the continuity equation for the higher conserved charges which imply
integrability.

The tool we used was the evaluation of level spacing distribution for a finite
chain using exact diagonalisation. Since the extremal parts of the spectrum
have special properties due to locality of the Hamiltonian, they were
discarded with the middle two-thirds of levels kept. We then quantified the
cross-over between the integrable Poissonian and chaotic Wigner-Dyson
statistics in the form of a critical coupling, which was normalised to allow a
direct comparison between different perturbation. As an example for a strong
integrability breaking perturbation we considered next-to-nearest-neighbour
interaction (NNNI).

The behaviour of the critical coupling as a function of the volume was found
to be fully consistent with the suggestion in [8] that the current
perturbation only breaks integrability in higher order. For any fixed volume
$L$ the critical couplings were significantly higher for the current
perturbation compared to the NNNI case, and also their decay with the volume
was significantly slower. In particular, for the massless case while the
critical value of the NNNI perturbation follows the $1/L^{3}$ law found
previously in [30], for the current perturbation we found a $1/L^{2}$ decay.
In the massive regime the critical coupling decreases exponentially with the
volume, but again the decay for the current perturbation was significantly
slower than for the NNNI.

Since the $1/L^{3}$ law for massless models was claimed universal for strong
integrability breaking [30] (albeit without analytic support), it is tempting
to speculate that the exponent eventually depends on the order at which
integrability is broken. However, as already known and also argued for in
Subsection 3.3, perturbation theory cannot be applied to examine the cross-
over in the level spacing distribution if the level crossings resolved by the
integrability breaking perturbation are generic i.e. involve only two levels.
It is interesting to note that in quantum field theories obtained as
perturbation of conformal field theories the crossover to the Wigner-Dyson
behaviour takes place already in the perturbative regime [31]. The essential
difference with the case considered here is that the conformal spectrum the
levels are generically multiply degenerate.

To sum up, our results strongly support the observations in [8, 17] that
perturbations of spin chains by generalised currents correspond to a weak
integrability breaking which only happens at higher order in the perturbing
coupling. An interesting question left open is to clarify the dependence of
the cross-over on the order of integrability breaking.

### Acknowledgements

This work was partially supported by the National Research, Development and
Innovation Office (NKFIH) under the research grant K-16 No. 119204, and also
by the Fund TKP2020 IES (Grant No. BME-IE-NAT), under the auspices of the
Ministry for Innovation and Technology. G. T. was also supported by the the
National Research, Development and Innovation Office (NKFIH) via the Hungarian
Quantum Technology National Excellence Program, project no. 2017-1.2.1-NKP-
2017-00001.

## References

  * [1] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, “Emergent Hydrodynamics in Integrable Quantum Systems Out of Equilibrium,” Physical Review X 6 (2016) 041065, arXiv:1605.07331 [cond-mat.stat-mech].
  * [2] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, “Transport in Out-of-Equilibrium X X Z Chains: Exact Profiles of Charges and Currents,” Phys. Rev. Lett. 117 (2016) 207201, arXiv:1605.09790 [cond-mat.stat-mech].
  * [3] D.-L. VU and T. Yoshimura, “Equations of state in generalized hydrodynamics,” SciPost Physics 6 (2019) 023, arXiv:1809.03197 [cond-mat.stat-mech].
  * [4] A. Urichuk, Y. Oez, A. Klümper, and J. Sirker, “The spin Drude weight of the XXZ chain and generalized hydrodynamics,” SciPost Physics 6 (2019) 005.
  * [5] X. Cao, V. B. Bulchandani, and H. Spohn, “The GGE averaged currents of the classical Toda chain,” Journal of Physics A Mathematical General 52 (2019) 495003, arXiv:1905.04548 [cond-mat.stat-mech].
  * [6] M. Borsi, B. Pozsgay, and L. Pristyák, “Current Operators in Bethe Ansatz and Generalized Hydrodynamics: An Exact Quantum-Classical Correspondence,” Physical Review X 10 (2020) 011054, arXiv:1908.07320 [cond-mat.stat-mech].
  * [7] B. Pozsgay, “Algebraic construction of current operators in integrable spin chains,” Phys. Rev. Lett. 125 (2020) 070602, arXiv:2005.06242 [cond-mat.stat-mech].
  * [8] B. Pozsgay, “Current operators in integrable spin chains: lessons from long range deformations,” SciPost Physics 8 (2020) 016, arXiv:1910.12833 [cond-mat.stat-mech].
  * [9] N. Beisert, “The dilatation operator of N=4 super Yang Mills theory and integrability,” Physics Reports 405 (2004) 1–202.
  * [10] T. Bargheer, N. Beisert, and F. Loebbert, “LETTER: Boosting nearest-neighbour to long-range integrable spin chains,” Journal of Statistical Mechanics: Theory and Experiment 2008 (2008) L11001, arXiv:0807.5081 [hep-th].
  * [11] T. Bargheer, N. Beisert, and F. Loebbert, “Long-range deformations for integrable spin chains,” Journal of Physics A Mathematical General 42 (2009) 285205, arXiv:0902.0956 [hep-th].
  * [12] B. Pozsgay, Y. Jiang, and G. Takács, “$T{\bar{T}}$-deformation and long range spin chains,” Journal of High Energy Physics 2020 (2020) 92, arXiv:1911.11118 [hep-th].
  * [13] F. A. Smirnov and A. B. Zamolodchikov, “On space of integrable quantum field theories,” Nuclear Physics B 915 (2017) 363–383, arXiv:1608.05499 [hep-th].
  * [14] A. Cavaglià, N. Gromov, and F. Levkovich-Maslyuk, “Separation of variables and scalar products at any rank,” Journal of High Energy Physics 2019 (2019) 52, arXiv:1907.03788 [hep-th].
  * [15] A. B. Zamolodchikov, “From tricritical Ising to critical Ising by thermodynamic Bethe ansatz,” Nuclear Physics B 358 (1991) 524–546.
  * [16] G. Mussardo and P. Simon, “Bosonic-type S-matrix, vacuum instability and CDD ambiguities,” Nuclear Physics B 578 (2000) 527–551, arXiv:hep-th/9903072 [hep-th].
  * [17] J. Durnin, M. Bhaseen, and B. Doyon, “Non-Equilibrium Dynamics and Weakly Broken Integrability,” arXiv e-prints (2020) , arXiv:2004.11030 [cond-mat.stat-mech].
  * [18] D. A. Rabson, B. N. Narozhny, and A. J. Millis, “Crossover from Poisson to Wigner-Dyson level statistics in spin chains with integrability breaking,” Phys. Rev. B 69 (2004) 054403, arXiv:cond-mat/0308089 [cond-mat].
  * [19] M. Rigol and L. F. Santos, “Quantum chaos and thermalization in gapped systems,” Phys. Rev. A 82 (2010) 011604, arXiv:1003.1403 [cond-mat.stat-mech].
  * [20] L. F. Santos and M. Rigol, “Onset of quantum chaos in one-dimensional bosonic and fermionic systems and its relation to thermalization,” Phys. Rev. E 81 (2010) 036206, arXiv:0910.2985 [cond-mat.stat-mech].
  * [21] L. F. Santos and M. Rigol, “Localization and the effects of symmetries in the thermalization properties of one-dimensional quantum systems,” Phys. Rev. E 82 (2010) 031130, arXiv:1006.0729 [cond-mat.stat-mech].
  * [22] G. P. Brandino, R. M. Konik, and G. Mussardo, “Energy level distribution of perturbed conformal field theories,” Journal of Statistical Mechanics: Theory and Experiment 2010 (2010) 07013, arXiv:1004.4844 [cond-mat.stat-mech].
  * [23] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, “From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics,” Advances in Physics 65 (2016) 239–362, arXiv:1509.06411 [cond-mat.stat-mech].
  * [24] W. Beugeling, R. Moessner, and M. Haque, “Off-diagonal matrix elements of local operators in many-body quantum systems,” Phys. Rev. E 91 (2015) 012144, arXiv:1407.2043 [cond-mat.stat-mech].
  * [25] W. Beugeling, A. Bäcker, R. Moessner, and M. Haque, “Statistical properties of eigenstate amplitudes in complex quantum systems,” Phys. Rev. E 98 (2018) 022204, arXiv:1710.11433 [cond-mat.stat-mech].
  * [26] M. Haque, P. A. McClarty, and I. M. Khaymovich, “Entanglement of mid-spectrum eigenstates of chaotic many-body systems – deviation from random ensembles,” arXiv e-prints (2020) arXiv:2008.12782, arXiv:2008.12782 [cond-mat.stat-mech].
  * [27] B. Bertini, P. Kos, and T. Prosen, “Exact Spectral Form Factor in a Minimal Model of Many-Body Quantum Chaos,” Phys. Rev. Lett. 121 (2018) 264101, arXiv:1805.00931 [nlin.CD].
  * [28] K. Kudo and T. Deguchi, “Unexpected non-Wigner behavior in level-spacing distributions of next-nearest-neighbor coupled XXZ spin chains,” Phys. Rev. B 68 (2003) 052510, arXiv:cond-mat/0308099 [cond-mat.stat-mech].
  * [29] K. Kudo and T. Deguchi, “Level Statistics of XXZ Spin Chains with Discrete Symmetries: Analysis through Finite-size Effects,” Journal of the Physical Society of Japan 74 (2005) 1992–2000, arXiv:cond-mat/0409761 [cond-mat.stat-mech].
  * [30] R. Modak, S. Mukerjee, and S. Ramaswamy, “Universal power law in crossover from integrability to quantum chaos,” Phys. Rev. B 90 (2014) 075152, arXiv:1309.1865 [cond-mat.str-el].
  * [31] M. Srdinsek, T. Prosen, and S. Sotiriadis, “Signatures of Chaos in Non-integrable Models of Quantum Field Theory,” arXiv e-prints (2020) arXiv:2012.08505, arXiv:2012.08505 [cond-mat.stat-mech].


# Characterizing quantum detectors by Wigner functions

Rajveer Nehra1<EMAIL_ADDRESS>Kevin Valson Jacob2<EMAIL_ADDRESS>
1Department of Physics, University of Virginia, 382 McCormick Rd,
Charlottesville, VA 22904-4714, USA  2Hearne Institute for Theoretical Physics
and Department of Physics & Astronomy, Louisiana State University, Baton
Rouge, Louisiana 70803, USA

###### Abstract

We propose a method for characterizing a photodetector by directly
reconstructing the Wigner functions of the detector’s Positive-Operator-Value-
Measure (POVM) elements. This method extends the works of S. Wallentowitz and
Vogel [Phys. Rev. A 53, 4528 (1996)] and Banaszek and Wódkiewicz [Phys. Rev.
Lett. 76, 4344 (1996)] for quantum state tomography via weak-field homodyne
technique to characterize quantum detectors. The scheme uses displaced thermal
mixtures as probes to the detector and reconstructs the Wigner function of the
photodetector POVM elements from its outcome statistics. In order to make the
reconstruction robust to the experimental noise, we use techniques from
quadratic convex optimization.

## Introduction

Photodetection has been making consistent progress with rapidly developing
optical quantum technology Hadfield (2009); Silberhorn (2007); Wang  _et al._
(2018); Ramilli  _et al._ (2010); Mattioli  _et al._ (2015); Rosenberg  _et
al._ (2005). Not only do detectors provide us with deeper insights on the
quantum behaviour of light by allowing us to perform precise measurements,
they also are an integral part of quantum technology such as quantum
computing, quantum enhanced metrology, and quantum communication Knill  _et
al._ (2001); O’Brien (2007); Kok  _et al._ (2007); Lloyd and Braunstein
(1999); Clauser and Shimony (1978); Aspect  _et al._ (1982).

For every quantum detector, one can associate a set of measurement operators
$\\{M_{k}\\}$ called as Positive Operator Value Measures (POVMs). When such a
device measures a quantum state $\rho$, the probability of observing an
outcome ‘$k$’ is

$p(k)_{\rho}=\mathrm{Tr}[\rho M_{k}].$ (1)

Since probabilities are non-negative and sum to one, POVM elements are
positive semi-definite and satisfy the completeness property
$\sum_{k=0}^{K-1}M_{k}=\mathbb{1}$. This implies that a set of POVM elements
completely describes the measurement device. Therefore, in order to
characterize a detector, we have to determine its POVM set.

In order to identify the POVM elements of a detector, one can invert Eq. (1)
which is known as Quantum Detector Tomography (QDT) Lundeen  _et al._ (2009);
Silberhorn (2007). In optical QDT, light prepared in a set of known
tomographically complete states a.k.a probes is incident on the detector to be
characterized. The probabilities of different measurement outcomes is then
used to characterize the detector.

One such possible set of probes is composed of coherent states
$\\{|\alpha\rangle\langle\alpha|\\}$. With a coherent state $|\alpha\rangle$
as the probe, the probability of outcome $k$ is given as

$p(k)_{|\alpha\rangle}=\text{Tr}[|\alpha\rangle\langle\alpha|M_{k}]=\pi
Q_{M_{k}}(\alpha),$ (2)

where $Q_{M_{k}}(\alpha)$ is the Husimi $Q$ quasi-probability distribution
corresponding to the detector POVM element $M_{k}$. Therefore, one can simply
reconstruct the $Q$ functions for POVM elements directly from the measurement
statistics. Since $Q_{M_{k}}(\alpha)$ has complete information about the POVM
element $M_{k}$, ideally it could be used to predict the measurement outcomes
for an arbitrary quantum state as follows: consider a quantum state $\rho$
represented in Glauber–Sudarshan $P$ representation as

$\rho=\int P_{\rho}(\alpha)|\alpha\rangle\langle\alpha|d^{2}\alpha,$ (3)

where $d^{2}\alpha:=d\text{Re}(\alpha)d\text{Im}(\alpha)$. The probability of
outcome $k$ can then be obtained using the Born rule as

$p(k)_{\rho}=\text{Tr}[\rho M_{k}]=\pi\int
P_{\rho}(\alpha)Q_{M_{k}}(\alpha)d^{2}\alpha.$ (4)

Therefore, by using the $Q$ representation for detector POVM elements and $P$
representation for the input quantum state, one can, in principle determine
the outcome probabilities corresponding to detector outcomes.

But this approach suffers from an inherent shortcoming due to the divergent
nature of $P$ functions for nonclassical states of the optical field Cahill
and Glauber (1969). In addition, as discussed in Lundeen  _et al._ (2009),
experimental errors and statistical noise during the experiments may distort
$Q$ functions resulting in nonphysical POVM elements. In order to alleviate
this shortcoming, it is beneficial to determine the POVM elements in some
basis, for instance photon-number, from the measurement statistics. Several
techniques have been proposed and demonstrated to reconstruct the POVM
elements in the photon-number basis Lundeen  _et al._ (2009); Grandi  _et al._
(2017); Brida  _et al._ (2012); Piacentini  _et al._ (2015); Humphreys  _et
al._ (2015); Zhang  _et al._ (2012). One can further represent the POVMs using
Wigner quasiprobability distribution functions as discussed in Section I.

Wigner functions, originally introduced by Eugene Wigner in 1932, provide a
useful method to visualize quantum states in phase space Wigner ; Gerry and
Knight (2004). The Wigner function corresponding to a quantum state of light
has been experimentally obtained using the balanced homodyne method as well as
photon number resolving measurements Lvovsky  _et al._ (2001); Morin  _et al._
(2012); Banaszek  _et al._ (1999); Bondani  _et al._ (2010); Sridhar  _et al._
(2014); Nehra  _et al._ (2019).

It is insightful to note that there exists a symmetry between quantum states
and measurement operators. We can see this from Eq. (1) wherein, due to the
cyclicity of trace, the roles of the state and the operator can be swapped.
This is the underlying relation which we exploit in order to identify the
Wigner functions of the measurement operators (POVM elements).

By obtaining the Wigner functions of the detector, any experimental
probability can be found in terms of the Wigner functions of the state as well
as of the detector, which are well-behaved unlike $P$ functions. Thus, Eq. (4)
can be written as

$p(k)_{\rho}=\text{Tr}[\rho M_{k}]=\int
W_{\rho}(\alpha)W_{M_{k}}(\alpha)d^{2}\alpha,$ (5)

where $W_{\rho}(\alpha)$ and $W_{M_{k}}(\alpha)$ are the Wigner functions of
quantum state $\rho$ and POVM element $M_{k}$ respectively.

In this work, we propose an alternative method for QDT by directly
reconstructing the Wigner quasiprobability functions corresponding to detector
POVM elements: it alleviates the need of finding the POVMs in the photon-
number basis. Apart from the fundamental interest in obtaining Wigner
functions of a detector, our scheme is particularly beneficial to study the
decoherence of a quantum detector by observing the behaviour of the POVM
Wigner functions in certain regions of the phase space D’Auria  _et al._
(2011).

This paper is organized as follows. In section I, we detail the method for
characterizing photodetectors using displaced thermal mixtures. In section II,
we then apply this method to photon-number-resolving detectors. Section III
discusses the resources required for characterizing a phase-insensitive
detector. We discuss the number of phase space points where Wigner function
should be experimentally measured in order to have a good confidence in
reconstruction. In section IV, we use convex optimization techniques to make
our reconstruction robust to noise. Finally, we note our conclusions in
section V.

##  I Method

We use the well known result that the Wigner function operator can be
represented in Fock space as

$\hat{W}(\alpha)=\frac{2}{\pi}\sum_{n=0}^{\infty}(-1)^{n}\hat{D}(\alpha)|n\rangle\langle
n|\hat{D}^{\dagger}(\alpha),$ (6)

where $\hat{D}(\alpha)=\exp({\alpha\hat{a}^{\dagger}-\alpha^{*}\hat{a}})$ is
the displacement operator with $\alpha\in\mathbb{C}$ Royer (1977). For a
detector, our aim is to characterize the Wigner functions corresponding to its
POVM elements. Since POVMs are self-adjoint positive semi-definite operators,
one can write the Wigner function of a POVM element $M_{k}$ as

$W_{M_{k}}(\alpha)=\frac{2}{\pi}\sum_{n=0}^{\infty}(-1)^{n}\mathrm{Tr}\left[M_{k}\hat{D}(\alpha)|n\rangle\langle
n|\hat{D}^{\dagger}(\alpha)\right],$ (7)

where, for simplicity, we define

$P_{M_{k}}^{(n)}(\alpha):=\mathrm{Tr}\left[M_{k}\hat{D}(\alpha)|n\rangle\langle
n|\hat{D}^{\dagger}(\alpha)\right].$ (8)

Although the sum in Eq. (7) has infinite terms but in practice one can
truncate it to ‘$n_{0}$’ as further terms don’t significantly contribute to
the sum. Thus, we have

$\displaystyle W_{M_{k}}(\alpha)$
$\displaystyle\approx\frac{2}{\pi}\sum_{n=0}^{n_{0}}(-1)^{n}P_{M_{k}}^{(n)}(\alpha).$
(9)

From Eq. (9) we can see that finding the Wigner function corresponding to
‘$M_{k}$’ amounts to finding out all these summands.

In this paper, we will restrict ourselves to phase-insensitive detectors for
simplicity. Such detectors have the Wigner functions of their POVM elements
rotationally symmetric around the origin, and hence can be characterized on
the real line alone. However, we note that this scheme is applicable to phase-
sensitive detectors also; and for such detectors, we have to choose $\alpha$
in the complex plane.

Figure 1: Schematic of the experimental setup.

###  I.1 Proposed experimental setup

Fig. 1 shows a schematic for our proposed experiment. A laser beam is split
into two beams at the first beamsplitter. One beam is used to generate thermal
mixtures. Thermal mixtures can be generated by randomzing the phase and
amplitude of the laser beam (coherent state). To achieve that, we use a
Variable Neutral Density Filter (VNDF) along with a Rotating Ground-Glass Disk
(RGGD). These allow us to generate thermal mixtures with a variable mean
photon-number Arecchi (1965). The other beam is used as a Local Oscillator
(LO) for phase space displacements. For displacements, we interfere thermal
mixtures with the LO at a highly unbalanced beamsplitter denoted as
Displacement beamsplitter (DBS) in the experiment schematic. A local
oscillator modulator (LOMD) is used for varying the displacements, which are
required to probe of the Wigner function over the whole phase space.

In order to do this, we consider $(n_{0}+1)$ distinct thermal mixtures given
as

$\rho^{(j)}=\sum_{n=0}^{\infty}p_{n}^{(j)}|n\rangle\langle n|$ (10)

where $j$ = $0,\ldots,n_{0}$ labels the thermal states, and
$p_{n}^{(j)}=\frac{\bar{n}^{n}_{j}}{(1+\bar{n}_{j})^{n+1}}$ is the Bose-
Einstein photon-number distribution of a thermal mixture $\rho^{(j)}$ with
mean photon-number $\bar{n}_{j}$. We then displace these thermal mixtures by
amplitude $\alpha$ which is, in general, a complex number. Then, the
probability of obtaining ‘$k$’ outcome with the displaced thermal input as
input is give by

$Q_{k}^{(j)}(\alpha)\approx\sum_{n=0}^{n_{0}}p_{n}^{(j)}P_{M_{k}}^{(n)}(\alpha),$
(11)

where we choose the thermal state such that the contribution to the RHS from
the omitted terms is negligible. This is possible because the thermal state
has an exponentially decreasing photon number distribution. In matrix form, we
can write Eq. (11) as

$\begin{pmatrix}Q_{k}^{(0)}(\alpha)\\\ Q_{k}^{(1)}(\alpha)\\\ \vdots\\\
Q_{k}^{(n_{0})}(\alpha)\end{pmatrix}=\begin{pmatrix}p_{0}^{(0)}&p_{1}^{(0)}&\ldots&p_{n_{0}}^{(0)}\\\
p_{0}^{(1)}&p_{1}^{(1)}&\ldots&p_{n_{0}}^{(1)}\\\ \vdots\\\
p_{0}^{(n_{0})}&p_{1}^{(n_{0})}&\ldots&p_{n_{0}}^{(n_{0})}\end{pmatrix}\begin{pmatrix}P_{M_{k}}^{(0)}(\alpha)\\\
P_{M_{k}}^{(1)}(\alpha)\\\ \vdots\\\
P_{M_{k}}^{(n_{0})}(\alpha)\end{pmatrix}.$ (12)

We can further write Eq. (12) compactly as

$\bf{Q=P\Pi^{\alpha}_{M_{k}}},$ (13)

where $\bf{Q}$ and $\bf{\Pi^{\alpha}_{M_{k}}}$ are vectors of length
$(n_{0}+1)$, and $\bf{P}$ is the probability distribution square matrix of
dimension $(n_{0}+1)\times(n_{0}+1)$. Thus by solving Eq. (13), we can
determine $\bf{\Pi^{\alpha}_{M_{k}}}$, which allows us to calculate the
summation in Eq. (7). To solve for $\bf{\Pi^{\alpha}_{M_{k}}}$, one needs to
solve the following convex quadratic optimization problem:

Minimize $\displaystyle\\{\bf{||Q-P\Pi^{\alpha}_{M_{k}}||_{2}}\\},$
$\displaystyle\begin{split}\text{Subject
to}&\enspace\bigg{\\{}1{\bf{\geq\Pi^{\alpha}_{M_{k}}}}\geq 0,\enspace\\\
&-1\leq\sum_{n=0}^{n_{0}}(-1)^{n}P^{(n)}_{M_{k}}(\alpha)\leq
1\bigg{\\}},\end{split}$ (14)

where $||.||$ is the $l_{2}$ norm defined as
$||V||_{2}=\big{(}\sum_{i}|V_{i}|^{2}\big{)}^{1/2}$ for a vector $V$.

The optimization constraints in Eq. (14) can be understood as follows. First,
the $n^{th}$ element $P^{(n)}_{M_{k}}(\alpha)$ of $\bf{\Pi^{\alpha}_{M_{k}}}$
is essentially the probability of getting $k$-click if a displaced $n$-photon
Fock state is incident to the detector. Therefore, we have $1\geq
P^{(n)}_{M_{k}}(\alpha)\geq 0$. Second, Wigner functions are well- defined and
bounded between $[-2/\pi,2/\pi]$ for a POVM element corresponding to a phase-
insensitive detector. This is because the POVM element in such a case is a
statistical mixture of projectors. In this case we can use the second
constraint in Eq. (14).

Thus solving this optimization allows us to determine the Wigner function at a
given phase space point $\alpha$. Further, we can repeat the process with
different displacement amplitudes to reconstruct the Wigner function over the
entire phase space. In practice, the Wigner functions of various detectors are
localized around the origin and vanishes to zero for large $\alpha$, so it is
unnecessary to displace the thermal states with arbitrarily large $\alpha$.

In the following section we numerically simulate this method for a phase-
insensitive detector. We hasten to add that this method is applicable to any
type of detector.

##  II Modelling a Photon-number-resolving detector

Figure 2: Wigner functions for POVM elements corresponding to zero-, one- and
two-photon detection events. Red curves are theoretically expected Wigner
functions and blue ones the reconstructed one using the proposed method. In
the top row, A, B, and C are for a perfect PNR detector; and in the bottom row
D, E, and F are for a PNR detector with detection efficiency $\eta=0.90$

.

In this section, we reconstruct the Wigner functions of a perfect and an
imperfect photon-number-resolving (PNR) detector. In general, a POVM element
corresponding to ‘$k$’ outcome can be written in the photon-number basis as

$M_{k}=\sum_{m,n=0}^{\infty}\langle m|M_{k}|n\rangle|m\rangle\langle n|,$ (15)

where $\langle m|M_{k}|n\rangle$ are the matrix elements of the POVM operator.
One can further simplify Eq. (15) for a PNR detector with no dark counts as

$\displaystyle M_{k}=\sum_{m=k}^{m_{0}}\langle
m|M_{k}|m\rangle|m\rangle\langle m|.$ (16)

Note that Eq. (16) differs from Eq. (15) in three ways. First, the POVM is
diagonal with entries $\langle m|M_{k}|m\rangle$, which are essentially the
probabilities of detecting ‘$k$’ photons given ‘$m$’ photons are incident to
the detector. Thus for a detector with detection efficiency $\eta$, we have

$p(k|m)=\langle m|M_{k}|m\rangle={m\choose k}\eta^{k}(1-\eta)^{m-k}.$ (17)

Second, we have truncated the sum to ‘$m_{0}$’ such that it exceeds the
photon-number at which saturates the detector. Third, the sum is starting from
‘$k$’ because with no dark counts noise, one would expect ‘$k$’ clicks only if
there are $m\geq k$ photons are incident on the detector.

Eq. (16) and Eq. (17) can be interpreted as follows: The POVM elements of a
perfect PNR detector are projectors $\Pi_{m}=|m\rangle\langle m|$. However,
for an imperfect detector, its efficiency $\eta<1$. If $m$ photons impinge on
such a detector, due to its non-unity detection efficiency, $k<m$ photons
results in a detection event contributing a factor of $\eta^{k}$ to the
probability of the event; while $(m-k)$ photons remain undetected contributing
a factor of $(1-\eta)^{m-k}$ to the probability of the event. Thus, such POVMs
are statistical mixtures of projective measurements.

In numerical simulations, we considered equidistant 51 displacement amplitudes
in $\alpha\in[-3.6,3.6]$, which allowed us to probe the Wigner function
uniformly over the entire region of phase space where the Wigner function is
non-vanishing. We used 50 equally spaced thermal states of mean photon-number
in $\bar{n}\in[0,4]$. For all of our simulations in open source Python module
QuTip Johansson  _et al._ (2013), we generally limited the dimension of the
Hilbert space to 50, and the sum in Eq. (16) was truncated with $m_{0}=50$ at
which point $P(k|m_{0})$ was of the order of $10^{-10}$ for $\eta=0.90$.

In Fig. 2, we plot the Wigner functions of one, two, and three photon
detections for a perfect detector and an imperfect detector with imperfections
as modelled in Eq. (16). From Fig. 2, we notice that the extrema of the Wigner
functions of imperfect detectors are closer to the origin than those of
perfect detectors. This is due to the contribution of higher order projectors
in the Wigner functions of imperfect detectors. In particular for the
imperfect single-photon detection event, we see a reduced negativity in the
Wigner function around the origin. This is due to the contribution of the
Wigner function of the two-photon detection event which is strongly positive
around origin. Similar arguments can be made for the the reduced positivity of
the Wigner function for zero- and two-photon detection event POVMs.

In our reconstruction, we have uniformly sampled the phase space. A natural
question that now arises is whether the number of points that needs to be
probed in the phase space can be reduced. We investigate this question in the
following section.

Figure 3: Left : Wigner function corresponding to a perfect single-photon
detection POVM determined by naive summation up to 15 terms of Eq. (7) (Red),
and using an Gaussian modulated quadratic fit near the origin (Blue). Black
points represent the phase space points where the Wigner function was probed
by the proposed method here. The latter approximates well the actual Wigner
function. Right : Wigner function corresponding to an imperfect single-photon
detection POVM with $\eta=0.90$.

##  III characterizing phase-insensitive detectors with polynomial resources

Although the method outline earlier is general, it had substantial resource
requirements as we had to uniformly sample over the phase space. However, this
requirement can be drastically reduced if we have the prior knowledge that the
detector is phase-insensitive, i.e. the representations of its POVM elements
are diagonal in the Fock basis. Note that the phase sensitivity of a detector
can easily be checked by varying the phase of the LO while keeping the
amplitude fixed. In this case, unlike a PNR detector, a phase sensitive
detector outputs different measurement statistics for different phases and
fixed amplitudes of the LO.

We recall that the Wigner functions of Fock states are Gaussian modulated
Laguerre polynomials Leonhardt (1997). This allows us to write the Wigner
function the POVM elememt ‘$M_{k}$’ of a PNR detector as

$W_{M_{k}}(\alpha)=\frac{2e^{-2|\alpha|^{2}}}{\pi}\sum_{m=0}^{m_{0}}(-1)^{m}p(k|m)\
{\mathcal{L}}_{m}(4|\alpha|^{2}),$ (18)

where ${\mathcal{L}}_{m}(x)$ represents the Laguerre polynomial of $m^{th}$
degree in $|\alpha|^{2}$. As the Wigner function is a function of
$|\alpha|^{2}$, it is symmetric around the origin, and can be fully
characterized on the real line. Since the Wigner function is a Gaussian
modulated polynomial, the problem of reconstructing it is reduced to finding
out a polynomial of degree $2m_{0}$ in $\alpha$ which requires us to find the
Wigner function only at $2m_{0}+1$ points.

As an example, we considered the POVM element corresponding to a single-photon
detection event for both perfect and imperfect PNR detectors. In Fig. 3, the
red curves show the POVM determined by the naive summation up to 15 terms of
Eq. (7); and the blue curves the reconstructed Wigner functions with black
points being the phase space coordinates where the Wigner function was probed.

We see that one needs to probe the Wigner function only at three points for a
perfect detector because the Laguerre polynomial
${\mathcal{L}}_{m=1}(4|\alpha|^{2})$ is quadratic in $\alpha$, and therefore
can be fully characterized using three distinct points. Likewise, the Wigner
function for an imperfect single-photon POVM can be reconstructed using only
11 distinct points (black points in Fig. 3) if we truncate the sum in Eq. (16)
at $m_{0}=5$ where $p(k|m)$ is of the order of $10^{-6}$. In this case, we
will have to reconstruct an Gaussian modulated polynomials of degree 10
because the last term in the Eq. (16) would be a projector, $|5\rangle\langle
5|$ with Wigner function given by Gaussian modulation of
${\mathcal{L}}_{m=5}(4|\alpha|^{2})$.

Note that finding the Gaussian modulated polynomial also works for a general
detector given by Eq. (15). However, instead of reconstructing the Wigner
function on the real line, we will have to reconstruct it in the complex plane
for which appropriate polynomial interpolation schemes have to be used Gasca
and Sauer (2000).

##  IV Robustness against experimental noise

In this section, we discuss the robustness of this method against experimental
noise. In general, inverting Eq. (13) is ill-conditioned as seen by the large
ratio of the largest and smallest singular values of the matrix $\bf{P}$. This
makes the reconstructed POVM elements extremely sensitive to small
fluctuations in the measurement statistics, and can lead to nonphysical POVMs.

However, the effects of ill-conditioning can be remarkably suppressed by
adding a regularization to the optimization problem. Several types of
regularization techniques are discussed in detail in Lundeen  _et al._ (2009),
and for this work we use Tikhonov regularization Il'in (1967). Using this
technique, inverting Eq. (13) can be mathematically formulated as the
following optimization problem:

Minimize
$\displaystyle\\{{||\bf{Q-P\Pi^{\alpha}_{M_{k}}}||_{2}+\gamma||\bf{\Pi^{\alpha}_{M_{k}}||_{2}}}\\},$
$\displaystyle\begin{split}\text{Subject
to}&\enspace\bigg{\\{}1\geq{\bf{\Pi^{\alpha}_{M_{k}}}}\geq 0,\enspace\\\
&-1\leq\sum_{n=0}^{n_{0}}(-1)^{n}P^{(n)}_{M_{k}}(\alpha)\leq
1\bigg{\\}},\end{split}$ (19)

where $\gamma$ is the regularization parameter. Solving this problem
translates to a convex quadratic optimization which can be efficiently solved
using a semi-definite problem solver, for instance, the Python package CVXOPT.

Figure 4: Blue: Reconstructed Wigner functions using regularization for
zero-,one- and two-photon detection event of a detector with $\eta=1$. Red
curves are theoretically expected Wigner functions. Gray areas are error
$(1\sigma)$ obtained using $N=40$ iterations. Bottom right: Dashed-diamond
curve illustrates the robustness of the reconstruction against regularization
parameter $\gamma$ and black solid line is without regularization, i.e,
$\gamma=0$.

In order to simulate the presence of noise in our reconstruction, we introduce
noise in the LO’s amplitude $|\alpha|$. We model this noise as a Gaussian
distribution of mean zero and standard deviation $\sigma=1\%|\alpha|^{2}$.
This is the typical noise level present in currently available stabilized
lasers. Therefore, the displacement amplitudes are $(\alpha_{1}+\delta
d_{1},\alpha_{2}+\delta d_{2}\dots,\alpha_{max}+\delta d_{max})$, where each
$\delta d_{i}$ is a random variable sampled from the Gaussian distribution. To
further reduce the effects of the fluctuations, we average the Wigner
functions obtained over $N=40$ iterations of the optimization. As a result, we
get

$\overline{W}_{{M_{k}}}(\alpha)=\frac{\sum_{j=1}^{N}W^{j}_{{M_{k}}}(\alpha+\delta{\alpha}_{j})}{N}.$
(20)

Having obtained $\overline{W}_{{M_{k}}}(\alpha)$, we then we utilize robust
nonlinear regression methods to further suppress the fluctuations. We recall
that for a phase insensitive detector, the POVMs are Gaussian modulated
polynomials of degree $2m_{0}$ in $\alpha$, where $m_{0}$ is the saturation
limit given in Eq. (18). Therefore, once we have experimentally probed the
Wigner function at $2m_{0}+1$ distinct points of the phase space, we could
simply fit a Gaussian modulated polynomial of degree $2m_{0}$ in $\alpha$ to
reconstruct the Wigner function over the entire phase space. Keeping that in
mind, we set an optimization problem as:

Minimize:
$\displaystyle\left\\{\frac{1}{2}\sum_{i=1}L\left[\left(e^{-2|\alpha_{i}|^{2}}\text{Poly}(2m_{0},\alpha_{i})-\overline{W}_{{M_{k}}}(\alpha_{i})\right)^{2}\right]\right\\},$

where $L$ is defined as

$L(y)=2(\sqrt{1+y}-1),$ (21)

and $\text{Poly}(2m_{0},\alpha_{i})$ is a polynomial of degree $2m_{0}$. Note
that this approach of finding the Gaussian modulated polynomial has an
advantage of not being biased unlike the simple least-square fitting method
which tends to significantly bias in order to avoid high residuals in the data
Triggs  _et al._ (2000).

We further evaluate the quality of reconstruction method by using the relative
error defined with $l_{2}$ norm as:

$\triangle:=\frac{||W^{theory}_{M_{k}}(\alpha)-W^{reconstruted}_{M_{k}}(\alpha)||_{2}}{||W^{theory}_{M_{k}}(\alpha)||_{2}}$
(22)

The result of our reconstruction is shown in Fig. 4. Since the fluctuations
grow with increasing local oscillator amplitude, the reconstruction of the
Wigner function around the origin of phase space is the least disturbed, but
with higher displacements the fluctuations grow stronger as seen in Fig. 4.
Therefore, it may be beneficial to probe the Wigner function around the origin
densely, and sparsely at the higher displacements, in particular $|\alpha|>1$.
Note that probing near the origin doesn’t undermine the quality of
reconstruction as long as we probe the Wigner function at $2m_{0}+1$ distinct
points because we need only $2m_{0}+1$ distinct points to reconstruct a
polynomial of degree $2m_{0}$ as seen in Fig. 3.

In fact, we can further exploit the rotational symmetry of the POVMs
corresponding to phase insensitive detector, which means the Wigner function
at $\alpha$ has the same value at $-\alpha$. This allows us to only probe the
Wigner function at $m_{0}+1$ distinct points to fully characterize a quantum
detector that saturates at the photon-number $m_{0}$. However, in this work we
numerically probe the phase space at equidistant displacement amplitudes.

We now investigate how sensitive our reconstruction is to the choice of
$\gamma$. To evaluate that, we calculate the relative error defined in Eq.
(22) for several values of $\gamma\in$ $[10^{-4},0.012]$. The result is
illustrated on the bottom right in Fig. 4 for the POVM element corresponding
to $n=1$ and $\eta=0.90$. We can clearly see that even if we vary $\gamma$ by
an order of magnitude (from $10^{-3}$ to $10^{-2}$), the relative error only
changes by less than one percent. This shows that there is sufficient freedom
in the choice of $\gamma$.

##  V Conclusions

We have developed a method for characterizing photodetectors by experimentally
obtaining the Wigner functions corresponding to the POVMs describing the
detector measurements. The proposed experimental scheme is simple and easily
accessible, in particular, for phase insensitive detector. Augmented with
convex quadratic optimization and robust nonlinear fitting techniques, we
demonstrated its robustness to the experimental fluctuations.

Future work on this method may involve an account for mode mismatch between
the local oscillator and the optical mode of thermal mixtures. This direction
of research is motivated by the fact that unlike in the balanced homodyne
technique, mode mismatch cannot simply be treated as losses in this method.

## Acknowledgements

KVJ acknowledges the support from Louisiana State University System Board of
Regents via an Economic Development Assistantship. RN was supported by NSF
grant PHY-1708023. The authors are grateful to Prof. Olivier Pfister and
Carlos Andres Gonzalez Arciniegas for the feedback on the manuscript. RN is
indebted to Prof. Jonathan P. Dowling for financially supporting the visit to
LSU. The authors thank Miller Eaton, Austin P. Lund, Arshdeep Sekhon, and
Vikesh Siddhu for helpful discussions.

The authors contributed equally to this work.

## References

  * Hadfield (2009) R. H. Hadfield, Nature Photonics 3, 696 EP (2009).
  * Silberhorn (2007) C. Silberhorn, Contemporary Physics 48, 143 (2007).
  * Wang  _et al._ (2018) Y. Wang, Z. Wang, Q. Yu, X. Xie, T. Posavitz, M. Jacob-Mitos, A. Ramaswamy, E. J. Norberg, G. A. Fish, and A. Beling, IEEE Journal of Selected Topics in Quantum Electronics 24, 1 (2018).
  * Ramilli  _et al._ (2010) M. Ramilli, A. Allevi, V. Chmill, M. Bondani, M. Caccia, and A. Andreoni, J. Opt. Soc. Am. B 27, 852 (2010).
  * Mattioli  _et al._ (2015) F. Mattioli, Z. Zhou, A. Gaggero, R. Gaudio, S. Jahanmirinejad, D. Sahin, F. Marsili, R. Leoni, and A. Fiore, Superconductor Science and Technology 28, 104001 (2015).
  * Rosenberg  _et al._ (2005) D. Rosenberg, A. E. Lita, A. J. Miller, and S. W. Nam, Phys. Rev. A 71, 061803 (2005).
  * Knill  _et al._ (2001) E. Knill, R. Laflamme, and G. J. Milburn, Nature 409, 46 (2001).
  * O’Brien (2007) J. L. O’Brien, Science 318, 1567 (2007).
  * Kok  _et al._ (2007) P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowling, and G. J. Milburn, Rev. Mod. Phys. 79, 135 (2007).
  * Lloyd and Braunstein (1999) S. Lloyd and S. L. Braunstein, Phys. Rev. Lett. 82, 1784 (1999).
  * Clauser and Shimony (1978) J. F. Clauser and A. Shimony, Rep. Prog. Phys. 41, 1881 (1978).
  * Aspect  _et al._ (1982) A. Aspect, P. Grangier, and G. Roger, Phys. Rev. Lett. 49, 91 (1982).
  * Lundeen  _et al._ (2009) J. S. Lundeen, A. Feito, H. Coldenstrodt-Ronge, K. L. Pregnell, C. Silberhorn, T. C. Ralph, J. Eisert, M. B. Plenio, and I. A. Walmsley, Nat. Phys. 5, 27 (2009).
  * Cahill and Glauber (1969) K. E. Cahill and R. J. Glauber, Phys. Rev. 177, 1882 (1969).
  * Grandi  _et al._ (2017) S. Grandi, A. Zavatta, M. Bellini, and M. G. A. Paris, New J. Phys. 19, 053015 (2017).
  * Brida  _et al._ (2012) G. Brida, L. Ciavarella, I. P. Degiovanni, M. Genovese, L. Lolli, M. G. Mingolla, F. Piacentini, M. Rajteri, E. Taralli, and M. G. A. Paris, New J. Phys. 14, 085001 (2012).
  * Piacentini  _et al._ (2015) F. Piacentini, M. P. Levi, A. Avella, M. López, S. Kück, S. V. Polyakov, I. P. Degiovanni, G. Brida, and M. Genovese, Opt. Lett. 40, 1548 (2015).
  * Humphreys  _et al._ (2015) P. C. Humphreys, B. J. Metcalf, T. Gerrits, T. Hiemstra, A. E. Lita, J. Nunn, S. W. Nam, A. Datta, W. S. Kolthammer, and I. A. Walmsley, New J. Phys. 17, 103044 (2015).
  * Zhang  _et al._ (2012) L. Zhang, A. Datta, H. B. Coldenstrodt-Ronge, X.-M. Jin, J. Eisert, M. B. Plenio, and I. A. Walmsley, New J. Phys. 14, 115005 (2012).
  * (20) E. P. Wigner,  _Quantum Theory of Angular Momentum_ , edited by L. C. Biedenharn and H. V. Dam (Academic Press, (New York, $1965$)).
  * Gerry and Knight (2004) C. Gerry and P. Knight,  _Introductory Quantum Optics_ (Cambridge University Press, 2004).
  * Lvovsky  _et al._ (2001) A. I. Lvovsky, H. Hansen, T. Aichele, O. Benson, J. Mlynek, and S. Schiller, Phys. Rev. Lett. 87, 050402 (2001).
  * Morin  _et al._ (2012) O. Morin, V. D’Auria, C. Fabre, and J. Laurat, Opt. Lett. 37, 3738 (2012).
  * Banaszek  _et al._ (1999) K. Banaszek, C. Radzewicz, K. Wódkiewicz, and J. S. Krasiński, Phys. Rev. A 60, 674 (1999).
  * Bondani  _et al._ (2010) M. Bondani, A. Allevi, and A. Andreoni, J. Opt. Soc. Am. B 27, 333 (2010).
  * Sridhar  _et al._ (2014) N. Sridhar, R. Shahrokhshahi, A. J. Miller, B. Calkins, T. Gerrits, A. Lita, S. W. Nam, and O. Pfister, J. Opt. Soc. Am. B 31, B34 (2014).
  * Nehra  _et al._ (2019) R. Nehra, A. Win, M. Eaton, N. Sridhar, R. Shahrokhshahi, T. Gerrits, A. Lita, S. W. Nam, and O. Pfister, (2019), arXiv:1906.02093 .
  * D’Auria  _et al._ (2011) V. D’Auria, N. Lee, T. Amri, C. Fabre, and J. Laurat, Phys. Rev. Lett. 107, 050504 (2011).
  * Royer (1977) A. Royer, Phys. Rev. A 15, 449 (1977).
  * Arecchi (1965) F. T. Arecchi, Phys. Rev. Lett. 15, 912 (1965).
  * Johansson  _et al._ (2013) J. R. Johansson, P. D. Nation, and F. Nori, Computer Physics Communications 184, 1234 (2013).
  * Leonhardt (1997) U. Leonhardt,  _Measuring the Quantum State of Light_ (Cambridge University Press, Cambridge, U.K., 1997).
  * Gasca and Sauer (2000) M. Gasca and T. Sauer, Advances in Computational Mathematics 12, 377 (2000).
  * Il'in (1967) V. A. Il'in, Russian Mathematical Surveys 22, 142 (1967).
  * Triggs  _et al._ (2000) B. Triggs, P. F. McLauchlan, R. I. Hartley, and A. W. Fitzgibbon, in  _Vision Algorithms: Theory and Practice_, edited by B. Triggs, A. Zisserman, and R. Szeliski (Springer Berlin Heidelberg, Berlin, Heidelberg, 2000) pp. 298–372.


NLO inclusive $J/\psi$ photoproduction at large $P_{T}$ at HERA and the EIC

C. Flore1$\star$, J.-P. Lansberg1, H.-S. Shao2 and Y. Yedelkina1

1 Université Paris-Saclay, CNRS, IJCLab, 91405 Orsay, France

2 Laboratoire de Physique Théorique et Hautes Energies, UMR 7589, Sorbonne
Université et CNRS, 4 place Jussieu, 75252 Paris Cedex 05, France

*<EMAIL_ADDRESS>
Proceedings for the XXVIII International Workshop
on Deep-Inelastic Scattering and Related Subjects,
Stony Brook University, New York, USA, 12-16 April 2021
10.21468/SciPostPhysProc.?

## Abstract

We study inclusive $J/\psi$ photoproduction at NLO at large $P_{T}$ at HERA
and the EIC. Our computation includes NLO QCD leading-$P_{T}$ corrections, QED
contributions via an off-shell photon as well as those from $J/\psi$+charm
channels. For the latter, we employ the variable-flavour-number scheme. Our
results are found to agree with the latest HERA data by H1 and provide, for
the first time, a reliable estimate of the EIC reach for such a measurement.
Finally, we demonstrate the observability of $J/\psi$+charm production and the
sensitivy to probe the non-perturbative charm content of the proton at high
$x$, also known as intrinsic charm, at the EIC.

##  1 Introduction

Inclusive $J/\psi$ photoproduction, when an almost on-shell photon hits and
breaks a proton producing the $J/\psi$, is a useful tool for understanding the
quarkonium-production mechanism and to learn more about the gluon content of
the proton.

In Ref. [1], we presented a comprehensive analysis on the NLO $J/\psi$
photoproduction in the large-$P_{T}$ region ($P_{T}\gg M_{J/\psi}$) which was
seldom studied at HERA [2, 3, 4, 5, 6, 7, 8], and which could be
experimentally studied in more details at the future US Electron-Ion Collider
(EIC) [9]. The analysis is based on the Colour-Singlet Model (CSM) [10], which
is the leading-$v$ contribution of Non-Relativistic QCD (NRQCD) [11].

At variance with previous works on inclusive $J/\psi$ photoproduction, where
only full QCD contributions were considered, we also included two, so far
overlooked, contributions: the associated $J/\psi$+charm production, and the
pure QED contribution, where the $J/\psi$ is produced by an off-shell photon,
that happens to be relevant at large $P_{T}$. The full list of the considered
partonic subprocesses is the following:

  1. (a)

$\gamma+g\to\mbox{$J/\psi$}+g(+g)$;

  2. (b)

$\gamma+\\{q,\bar{q}\\}\to\mbox{$J/\psi$}+\\{q,\bar{q}\\}(+g)\,$ ($q=u,d,s$);

  3. (c)

$\gamma+g\to\mbox{$J/\psi$}+c+\bar{c}$ and
$\gamma+\\{c,\bar{c}\\}\to\mbox{$J/\psi$}+\\{c,\bar{c}\\}$.

$\gamma+g\to\mbox{$J/\psi$}+g$ represents the
$\mathcal{O}(\alpha\alpha_{s}^{2})$ LO QCD contribution, while
$\gamma+q\to\mbox{$J/\psi$}+q$ is the $\mathcal{O}(\alpha^{3})$ LO QED
subprocess. The $\mathcal{O}(\alpha\alpha_{s}^{3})$ contributions are given by
$\gamma\\{g,q,\bar{q}\\}\to\mbox{$J/\psi$}+\\{g,q,\bar{q}\\}$ $+g$. The sum of
$\mathcal{O}(\alpha\alpha_{s}^{2})$ and $\mathcal{O}(\alpha\alpha_{s}^{3})$
terms is the NLO QCD contribution. Finally, we note that the partonic
subprocesses (c) can be calculated respectively in the 3 and 4 Flavour Schemes
(3FS and 4FS). A proper treatment of such channels is required, and is given
by the LO Variable Flavour Number Scheme (LO VFNS) illustrated in Section 2.4
of Ref. [1], to which we guide the reader for more details.

##  2 $J/\psi$ photoproduction at finite $P_{T}$ and the
NLO${}^{\star}\,$approximation

We first start by recalling some elements of kinematics. We define
$s_{ep}=(P_{e}+P_{p})^{2}=4E_{e}E_{p}$ ($E_{e(p)}$ is the electron (proton)
beam energy) and $s_{\gamma p}=W_{\gamma p}^{2}=(P_{\gamma}+P_{p})^{2}$.
Introducing $x_{\gamma}$ as $P_{\gamma}=x_{\gamma}P_{e}$, it turns out that
$s_{\gamma p}=x_{\gamma}s_{ep}$. Another important variable, called
elasticity, is defined as $z=(P_{\mathcal{Q}}\cdot P_{p})/(P_{\gamma}\cdot
P_{p})$. $z$ reduces to the photon energy taken by the $J/\psi$ in the proton
rest frame, and can also be expressed as $z=(2\,E_{p}\,m_{T})/(W^{2}_{\gamma
p}\,e^{y})$ in terms of the $J/\psi$ rapidity $y$ (with $y$ and $E_{p}$
defined in the same frame) and the quarkonium transverse mass,
$m_{T}=\sqrt{m_{\mathcal{Q}}^{2}+P_{\mathcal{Q}T}^{2}}$.

In our evaluation of the NLO corrections to $J/\psi+g$, we resorted to the
NLO⋆ approximation [12, 13], which allows one to consider the leading-$P_{T}$
contributions by imposing a democratic, lower cut on the invariant mass of
every pair of massless partons, denoted by $s_{ij}^{\rm min}$. Such an
approximation has already been tested in the case of hadroproduction (see e.g.
Ref. [13]), and it is automated within HELAC-Onia [14, 15], that we used for
our computations. As for any leading-$P_{T}$ topology associated to real NLO
emission $s_{ij}$ for any $i,j$ pair grows with $P_{T}$, the result becomes
insensitive to $s_{ij}^{\rm min}$, and only the subleading $P_{T}$
contributions exhibit a $\log(s_{ij}^{\rm min})$ dependence which is thus
$P_{T}$-power suppressed.

For testing the validity of the NLO⋆ approximation, we checked our approximate
calculation against a full NLO calculation by Butenschön and Kniehl [16]. It
turned out that a suitable interval for our analysis is $\sqrt{s_{ij}^{\rm
min}}/m_{c}\in\left[1:3\right]$ ($m_{c}=1.5$ GeV), and that $\sqrt{s_{ij}^{\rm
min}}=2m_{c}$ remarkably reproduces the complete NLO calculation. The results
are shown in Fig. 2 of Ref. [1]. In view of such a result, we revisited the
latest HERA data from H1 Collaboration and then extended our computation to
the case of the future EIC.

##  3 Results

In the following, we present our results for HERA and the future EIC. In every
plot, we show the following CS contributions: LO (blue) and
NLO${}^{\star}\,$QCD (red), LO QED (light blue), LO VFNS (green) and the
overall sum (orange). A feed-down contribution
$\psi^{\prime}\to\mbox{$J/\psi$}$ of $20\%$ is also taken into account111We
guide the reader to Section 3 of Ref. [1] for a detailed discussion about the
considered feed-downs..

For our predictions, we used the CT14NLO proton PDFs set [17]. The
corresponding theoretical uncertainty is automatically evaluated by HELAC-
Onia, as well as the factorisation- and renormalisation-scale uncertainties,
evaluated from an independent variation in the interval
$\mu_{F},\mu_{R}\in[1/2:2]\,\mu_{0}$, with $\mu_{0}=m_{T}$. The mass
uncertainty is obtained by varying $m_{c}$ in the range $m_{c}=1.5\pm 0.1$ GeV
for all but the LO QED channel. For this latter channel, the invariant mass of
the photon, $2m_{c}$, should coincide with $M_{J/\psi}$ and we chose
$m_{c}=1.55$ GeV. We have also considered
$\langle\mathcal{O}_{J/\psi}^{{}^{3\\!}S_{1}^{[1]}}\rangle=1.45$ GeV3.

We first start by revisiting HERA data. Our results are computed at H1 Run 2
kinematics [7]: $\sqrt{s_{ep}}=319$ GeV, $Q^{2}<2.5$ GeV2,
$P_{T}>1\,\text{GeV}$, $0.3<z<0.9$, $60\,\text{GeV}<W_{\gamma
p}<240\,\text{GeV}$. From Fig. 1, we note the following: 1) the LO QCD
contribution well describes the bulk of the data at low $P_{T}$; 2) the LO QED
contribution is smaller in size, but its spectrum is harder with respect to
the LO QCD one; 3) the LO VFNS $\mbox{$J/\psi$}+$charm contribution is not
negligible and matters at large $P_{T}$; 4) the QCD NLO${}^{\star}\,$is close
to the data points, and the overall sum nearly agrees with them. Such an
agreement is enhanced when subtracting the expected $b\to\mbox{$J/\psi$}$
feed-down from the data222See Appendix A of Ref. [1] for a detailed discussion
about the feed-downs from $b$ quarks.. This being said, we argue that the CSM
up to order $\alpha\alpha_{s}^{3}$ is able to reproduce HERA data. For this
reason, we restrict our EIC predictions to the CSM.

Figure 1: Comparison between H1 data [7] (black: inclusive yield; grey:
estimated prompt yield) and various CS contributions: LO (blue),
NLO${}^{\star}\,$(red), LO VNFS $J/\psi+c$ (green), LO QED (light blue) and
their combination (orange). The solid bands indicate the mass uncertainty
while the patterns display the scale uncertainty. Figure taken from Ref. [1].

Moving now to the EIC, we remark that this future collider will be able to run
at different $\sqrt{s_{ep}}$ with very large luminosities. We consider two
energy configurations: $E_{e}=5\,(18)$ GeV and $E_{p}=100\,(275)$ GeV,
resulting in $\sqrt{s_{ep}}=45\,(140)$ GeV. Also here, we apply some
kinematical cuts: $P_{T}>1$ GeV, $0.05<z<0.9$, $Q^{2}<1$ GeV2; we also
consider two different $W_{\gamma p}$ regions: $[10:40]$ GeV and $[20:80]$ GeV
for $\sqrt{s_{ep}}=45$ and $140$ GeV respectively.

Figure 2: Predictions for the future EIC at $\sqrt{s_{ep}}=45$ GeV (left) and
$\sqrt{s_{ep}}=140$ GeV (right). The calculation is performed adopting the
same $\mu_{F}$, $\mu_{R}$ and PDFs as Fig. 1 with the same meaning for bands.
Figure taken from Ref. [1].

By looking at Fig. 2, we can note some differences between the two kinematical
configurations. At $\sqrt{s_{ep}}=45$ GeV (Fig. 2, left panel), as $P_{T}$
increases, one enters the valence region, and the QED contribution becomes the
dominant one at the largest measurable $P_{T}\simeq 11$ GeV at
$\mathcal{L}=100$ fb-1. Moreover, $\gamma+q$ fusion contributes more than
$30\%$ for $P_{T}>8$ GeV, and the $\mbox{$J/\psi$}\,+$ unidentified charm
contribution is comparable to the $\gamma+g(q)$ fusion subprocesses. Hence,
both of the so far overlooked contributions will be relevant at the EIC.
Looking then at $\sqrt{s_{ep}}=140$ GeV (right panel in Fig. 2), one can note
that the yield is measurable up to $P_{T}\sim 18$ GeV. Also at this energy,
the QED contribution is the leading one at the largest reachable $P_{T}$,
while $\gamma+g$ fusion results the dominant contribution up to $P_{T}\sim 15$
GeV. More generally, it turns out that the production of $\mbox{$J/\psi$}+2$
hard partons (i.e.$\,\mbox{$J/\psi$}+\\{gg,qg,c\bar{c}\\}$) is dominant for
$P_{T}\sim 8-15$ GeV. This could lead to the observation of
$\mbox{$J/\psi$}+2$ jets with moderate $P_{T}$, with the leading jet1
recoiling on the $\mbox{$J/\psi$}+{\rm jet}_{2}$ pair.

###  3.1 $\mbox{$J/\psi$}\,+$ charm production

Another interesting aspect that could be studied at the future EIC is the
associated production of a $J/\psi$ and a charmed particle. In addition, one
may also check to which extent the $\mbox{$J/\psi$}+c$ yield could help
detecting a valence-like, non-perturbative charm content in the proton, to
which we refer to as Intrinsic Charm (IC) [18]. By considering the same LO
VFNS computation used so far, but now including a $10\%$ charm detection
efficiency $\varepsilon_{c}$ as follows:

$d\sigma^{\rm VFNS}=d\sigma^{\rm
3FS}\left[1-\left(1-\varepsilon_{c}\right)^{2}\right]+\left(d\sigma^{\rm
4FS}-d\sigma^{\rm CT}\right)\varepsilon_{c}$ (1)

(where $d\sigma^{\rm CT}$ is a proper counterterm to properly merge 3FS and
4FS contributions), we studied the $\mbox{$J/\psi$}+c$ yield at the two
different EIC energy configurations. To do so, we employed the CT14nnloIC PDFs
set [19], that encodes some eigensets with different IC effects: a “sea-like”
one (in green) and a “valence-like” one, also called “BHPS” (in red). The
central eigenset, to which we refer to as “no IC”, is depicted in blue. Note
that these IC effects are only affecting the 4FS contribution in Eq. (1).

By looking at Fig. 3, left panel, we see that at $\sqrt{s_{ep}}=45$ GeV the
$\mbox{$J/\psi$}+c$ yield is limited to low $P_{T}$ even with the largest
integrated luminosity. Nonetheless, the yield is clearly observable if
$\varepsilon_{c}=0.1$ with $\mathcal{O}(500,50,5)$ events for
$\mathcal{L}=(100,10,1)$ fb-1. On the other hand, at $\sqrt{s_{ep}}=140$ GeV
(Fig. 3, right panel), the $P_{T}$ range is up to $10$ GeV, and we expect
$\mathcal{O}(10^{3})$ events at $\mathcal{L}=100$ fb-1. Such events could be
observed by measuring a charmed jet. Finally, note that as the valence region
at high $x$ is not probed, no clear IC effect is visible at
$\sqrt{s_{ep}}=140$ GeV, while at $\sqrt{s_{ep}}=45$ GeV, we observe a
measurable effect, where the BHPS valence-like peak is visible.

Figure 3: Predictions for $\mbox{$J/\psi$}+$charm production at the future EIC
at $\sqrt{s_{ep}}=45$ GeV (left) and $\sqrt{s_{ep}}=140$ GeV (right) computed
with charm PDFs with “no IC” (blue), sea-like (green) and “BHPS” valence-like
(red) ICs. The lower panel shows the ratio to the “no IC” curves and its
relative uncertainty. The kinematical cuts are the same as in Fig. 2. A $10\%$
charm detection efficiency is considered when drawing the horizontal
observability lines. Figure taken from Ref. [1].

##  4 Conclusion

We have analysed inclusive, large-$P_{T}$ $J/\psi$ photoproduction at $ep$
colliders. By including new, so far overlooked partonic subprocesses and feed-
down contributions, we have shown that the CSM at
$\mathcal{O}(\alpha\alpha_{s}^{3})$ is able to describe the latest HERA data
measured by the H1 Collaboration. We then provided predictions for the future
EIC, where the new contributions, the $\mathcal{O}(\alpha^{3})$ LO QED and the
$\mbox{$J/\psi$}+$charm, can have a non-negligible role at this future
collider. Furthermore, we showed that EIC can represent an ideal playground to
constrain the proton charm content via the associated production of a $J/\psi$
and an identified charmed particle. Such studies can be surely extended to
other planned future colliders, such as the LHeC and FCC-eh [20, 21].

## Acknowledgements

We thank Y. Feng, M.A. Ozcelik, J.W. Qiu, I. Schienbein, H. Spiesberger for
useful discussions.

#### Funding information

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under grant agreement No. 824093 in order to
contribute to the EU Virtual Access NLOAccess. This work was also partly
supported by the French CNRS via the IN2P3 project GLUE@NLO, via the Franco-
Chinese LIA FCPPL (Quarkonium4AFTER), by the Paris-Saclay U. via the P2I
Department and by the P2IO Labex via the Gluodynamics project, the French ANR
under the grant ANR-20-CE31-0015 (PrecisOnium), and the CNRS IEA (GlueGraph).
L.Y. is supported by the EU Erasmus+ Paris-Saclay U.-Ukraine program.

## References

  * [1] C. Flore, J. P. Lansberg, H. S. Shao and Y. Yedelkina, “Large-$P_{T}$ inclusive photoproduction of $J/\psi$ in electron-proton collisions at HERA and the EIC,” Phys. Lett. B 811 (2020), 135926,arXiv:2009.08264 [hep-ph].
  * [2] H1 Collaboration, S. Aid et al., ‘‘Elastic and inelastic photoproduction of $J/\psi$ mesons at HERA,’’ Nucl. Phys. B 472 (1996) 3--31, arXiv:hep-ex/9603005.
  * [3] ZEUS Collaboration, J. Breitweg et al., ‘‘Measurement of inelastic $J/\psi$ photoproduction at HERA,’’ Z. Phys. C 76 (1997) 599--612, arXiv:hep-ex/9708010.
  * [4] ZEUS Collaboration, S. Chekanov et al., ‘‘Measurements of inelastic $J/\psi$ and $\psi^{\prime}$ photoproduction at HERA,’’ Eur. Phys. J. C 27 (2003) 173--188, arXiv:hep-ex/0211011.
  * [5] H1 Collaboration, C. Adloff et al., ‘‘Inelastic photoproduction of $J/\psi$ mesons at HERA,’’ Eur. Phys. J. C25 (2002) 25--39, arXiv:hep-ex/0205064 [hep-ex].
  * [6] ZEUS Collaboration, S. Chekanov et al., ‘‘Measurement of J/psi helicity distributions in inelastic photoproduction at HERA,’’ JHEP 12 (2009) 007, arXiv:0906.1424 [hep-ex].
  * [7] H1 Collaboration, F. Aaron et al., ‘‘Inelastic Production of $J/\psi$ Mesons in Photoproduction and Deep Inelastic Scattering at HERA,’’ Eur. Phys. J. C 68 (2010) 401--420, arXiv:1002.0234 [hep-ex].
  * [8] ZEUS Collaboration, H. Abramowicz et al., ‘‘Measurement of inelastic $J/\psi$ and $\psi^{\prime}$ photoproduction at HERA,’’ JHEP 02 (2013) 071, arXiv:1211.6946 [hep-ex].
  * [9] A. Accardi et al., ‘‘Electron Ion Collider: The Next QCD Frontier: Understanding the glue that binds us all,’’ Eur. Phys. J. A 52 no. 9, (2016) 268, arXiv:1212.1701 [nucl-ex].
  * [10] E. L. Berger and D. L. Jones, ‘‘Inelastic Photoproduction of J/psi and Upsilon by Gluons,’’ Phys. Rev. D23 (1981) 1521--1530.
  * [11] G. T. Bodwin, E. Braaten, and G. P. Lepage, ‘‘Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium,’’ Phys. Rev. D51 (1995) 1125--1171, arXiv:hep-ph/9407339 [hep-ph]. [Erratum: Phys. Rev.D55,5853(1997)].
  * [12] P. Artoisenet, J. M. Campbell, J. Lansberg, F. Maltoni, and F. Tramontano, ‘‘$\Upsilon$ Production at Fermilab Tevatron and LHC Energies,’’ Phys. Rev. Lett. 101 (2008) 152001, arXiv:0806.3282 [hep-ph].
  * [13] J. Lansberg, ‘‘On the mechanisms of heavy-quarkonium hadroproduction,’’ Eur. Phys. J. C 61 (2009) 693--703, arXiv:0811.4005 [hep-ph].
  * [14] H.-S. Shao, ‘‘HELAC-Onia: An automatic matrix element generator for heavy quarkonium physics,’’ Comput. Phys. Commun. 184 (2013) 2562--2570, arXiv:1212.5293 [hep-ph].
  * [15] H.-S. Shao, ‘‘HELAC-Onia 2.0: an upgraded matrix-element and event generator for heavy quarkonium physics,’’ Comput. Phys. Commun. 198 (2016) 238--259, arXiv:1507.03435 [hep-ph].
  * [16] M. Butenschoen and B. A. Kniehl, ‘‘World data of $J/\psi$ production consolidate NRQCD factorization at NLO,’’ Phys. Rev. D 84 (2011) 051501, arXiv:1105.0820 [hep-ph].
  * [17] S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin, C. Schmidt, D. Stump, and C. Yuan, ‘‘New parton distribution functions from a global analysis of quantum chromodynamics,’’ Phys. Rev. D 93 no. 3, (2016) 033006, arXiv:1506.07443 [hep-ph].
  * [18] S. J. Brodsky, P. Hoyer, C. Peterson, and N. Sakai, ‘‘The Intrinsic Charm of the Proton,’’ Phys. Lett. 93B (1980) 451--455.
  * [19] T.-J. Hou, S. Dulat, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, C. Schmidt, J. Winter, K. Xie, and C. P. Yuan, ‘‘CT14 Intrinsic Charm Parton Distribution Functions from CTEQ-TEA Global Analysis,’’ JHEP 02 (2018) 059, arXiv:1707.00657 [hep-ph].
  * [20] LHeC Study Group Collaboration, J. Abelleira Fernandez et al., ‘‘A Large Hadron Electron Collider at CERN: Report on the Physics and Design Concepts for Machine and Detector,’’ J. Phys. G 39 (2012) 075001, arXiv:1206.2913 [physics.acc-ph].
  * [21] LHeC, FCC-he Study Group Collaboration, P. Agostini et al., ‘‘The Large Hadron-Electron Collider at the HL-LHC,’’ arXiv:2007.14491 [hep-ex].