# Reverse Carleson measures for Hardy spaces in the unit ball

Evgueni Doubtsov  Department of Mathematics and Computer Science, St.
Petersburg State University, Line 14th (Vasilyevsky Island), 29, St.
Petersburg 199178, Russia<EMAIL_ADDRESS>
###### Abstract.

Let $H^{p}=H^{p}(B_{d})$ denote the Hardy space in the open unit ball $B_{d}$
of $\mathbb{C}^{d}$, $d\geq 1$. We characterize the reverse Carleson measures
for $H^{p}$, $1<p<\infty$, that is, we describe all finite positive Borel
measures $\mu$, defined on the closed ball $\overline{B}_{d}$, such that

$\|f\|_{H^{p}}\leq c\|f\|_{L^{p}(\overline{B}_{d},\mu)}$

for all $f\in H^{p}(B_{d})\cap C(\overline{B}_{d})$ and a universal constant
$c>0$. Given a non-inner holomorphic function $b:B_{d}\to B_{1}$, we obtain
properties of the reverse Carleson measures for the de Branges–Rovnyak space
$\mathcal{H}(b)$.

This research was supported by the Russian Science Foundation (grant No.
19-11-00058), https://rscf.ru/project/19-11-00058/

Key words: Hardy spaces, reverse Carleson measures, de Branges–Rovnyak spaces

##  1\. Introduction

Let $B_{d}$ denote the open unit ball of $\mathbb{C}^{d}$, $d\geq 1$, and let
$\sigma=\sigma_{d}$ denote the normalized Lebesgue measure on the unit sphere
$\partial B_{d}$. For the unit disk $B_{1}$, the measure $\sigma_{1}$ and the
unit circle $\partial B_{1}$, we also use the symbols $\mathbb{D}$, $m$ and
$\partial\mathbb{D}$, respectively.

Let $\mathcal{H}ol(B_{d})$ denote the space of holomorphic functions in the
ball $B_{d}$. For $0<p<\infty$, the Hardy space $H^{p}=H^{p}(B_{d})$ consists
of functions $f\in\mathcal{H}ol(B_{d})$ such that

$\|f\|_{H^{p}}^{p}=\sup_{0<r<1}\int_{\partial
B_{d}}|f(r\zeta)|^{p}\,d\sigma(\zeta)<\infty.$

As usual, the Hardy space $H^{p}(B_{d})$, $p>0$, is identified with the space
$H^{p}(\partial B_{d})$, which consists of the corresponding boundary values.
Further details about the spaces ${H^{p}}(B_{d})$ and related objects are
contained, for example, in monographs [8] and [10].

###  1.1. Carleson measures

Let $M_{+}(B_{d})$ denote the space of all finite positive Borel measures
defined on $B_{d}$. Firstly, assume that $d=1$ and $\mu\in M_{+}(\mathbb{D})$.
The classical Carleson Theorem [3] states that $H^{p}(\mathbb{D})$ is embedded
into $L^{p}(\mathbb{D},\mu)$, that is, there exists a constant $C>0$ such that

(1.1) $\|f\|_{L^{p}(\mathbb{D},\mu)}\leq C\|f\|_{H^{p}}\quad\textrm{for all}\
f\in{H^{p}}$

if and only if $\mu$ satisfies the following Carleson condition: there exists
a constant $C>0$ such that

(1.2) $\mu(S_{I})\leq Cm(I)\quad\textrm{for all arcs}\
I\subset\partial\mathbb{D},$

where

$S_{I}=\left\\{z\in\overline{\mathbb{D}}:\ 1-m(I)\leq|z|\leq 1,\
\frac{z}{|z|}\in I\right\\}.$

For an arbitrary dimension $d\geq 1$, consider the sets

$Q=Q(\zeta,\delta)=\\{\xi\in\partial
B_{d}:|1-\langle\zeta,\xi\rangle|\leq\delta\\},\quad\zeta\in\partial B_{d},\
\delta>0.$

Observe that $Q$ is a closed ball with respect to the nonisotropic metric

$\rho(\zeta,\xi)=|1-\langle\zeta,\xi\rangle|^{\frac{1}{2}},\quad\zeta,\xi\in\partial
B_{d}.$

For $r>0$ and $Q=Q(\zeta,\delta)$, put $rQ=Q(\zeta,r\delta)$ by definition.

By the theorem on Carleson measures for the Hardy space in the unit ball (see,
for example, [10, Ch. 5]), the estimate

(1.3) $\|f\|_{L^{p}(\mathbb{D},\mu)}\leq C\|f\|_{H^{p}},\quad f\in{H^{p}},$

is equivalent to the following condition on the measure $\mu$:

(1.4) $\mu(S_{Q})\leq C\sigma(Q)\quad\textrm{for all nonisotropic balls}\
Q\subset\partial B_{d},$

where

$S_{Q}=\left\\{z\in B_{d}:\ 1-\sigma(Q)\leq|z|\leq 1,\ \frac{z}{|z|}\in
Q\right\\}$

is a standard Carleson window.

###  1.2. Reverse Carleson measures

Remark that Luecking [7] studied Carleson measures for Bergman spaces in the
unit disk together with those measures for which the corresponding inverse
inequality holds. In the case of several complex variables, we mention the
recent related results by Green and Wagner [4] about the dominating sets for
Bergman spaces on general domains in $\mathbb{C}^{d}$.

In the first part of the present paper, we investigate the reverse Carleson
inequality

(1.5) $\|f\|_{H^{p}}\leq C\|f\|_{L^{p}(\overline{B}_{d},\mu)},\quad f\in
C(\overline{B}_{d})\cap H^{p}(B_{d}),\ 1<p<\infty,$

for $\mu\in M_{+}(\overline{B}_{d})$.

Remark that the above problem is formulated under assumption that the measure
$\mu$ is defined on the closed ball $\overline{B}_{d}$. In fact, for $\mu\in
M_{+}(B_{d})$ this problem degenerates and has no solutions; see Theorem 1
stated below.

Observe that the theorem about the Carleson measures for ${H^{p}}(B_{d})$ also
has a formulation under the assumption $\mu\in M_{+}(\overline{B}_{d})$.
Indeed, in this case, condition (1.3) makes sense for the functions $f$ in the
set $C(\overline{B}_{d})\cap H^{p}(B_{d})$, which is dense in $H^{p}(B_{d})$.
Condition (1.4) implies that the restriction $\mu_{|\partial B_{d}}$ is
absolutely continuous with respect to Lebesgue measure $\sigma$, and the
Radon–Nikodym derivative ${d\mu_{|\partial B_{d}}}/{d\sigma}$ is bounded.

For $d=1$ and a measure $\mu\in M_{+}(\overline{\mathbb{D}})$ with property
(1.2), Lefèvre, Li et al. [6] proved that property (1.5) is equivalent to the
following condition:

$\mu(S_{I})\geq Cm(I)\quad\textrm{for all arcs}\ I\subset\partial\mathbb{D}.$

Further, this result was proved for all $\mu\in M_{+}(\overline{\mathbb{D}})$
in [5], where other equivalent properties have been also obtained.

The main result of the present paper about the reverse Carleson measures for
the Hardy space $H^{p}(B_{d})$ with $d\geq 1$ is Theorem 1 formulated in the
next section.

###  1.3. Reverse Carleson measures for Hardy spaces in the unit ball

Let $1<p<\infty$. For $w\in B_{d}$, consider the Cauchy kernels

$k_{w}(z)=\frac{1}{(1-\langle z,w\rangle)^{d}},\quad z\in B_{d},$

and the corresponding normalized functions

$K_{w}=\frac{k_{w}}{\|k_{w}\|_{{H^{p}}}}.$

Proposition 1.4.10 from [8] guarantees that

$\|k_{w}\|_{p}=(1-|w|^{2})^{-\frac{d}{q}},\quad\textrm{where}\
\frac{1}{p}+\frac{1}{q}=1.$

For $d=1$, the following result was obtained in [5].

######  Theorem 1.

Let $1<p<\infty$ and $\mu\in M_{+}(\overline{B}_{d})$, $d\geq 1$. Then the
following properties are equivalent.

  1. (i)

There exists a constant $C>0$ such that

$\int_{\overline{B}_{d}}|f|^{p}\,d\mu\geq C\|f\|_{H^{p}}^{p}\quad\textrm{for
all}\ f\in C(\overline{B}_{d})\cap H^{p}.$

  2. (ii)

There exists a constant $C>0$ such that

$\int_{\overline{B}_{d}}|K_{w}|^{p}\,d\mu\geq C\quad\textrm{for all}\ w\in
B_{d}.$

  3. (iii)

There exists a constant $C>0$ such that

$\mu(Q)\geq C\sigma(Q)\quad\textrm{for all nonisotropic balls}\
Q\subset\partial B_{d}.$

It is mentioned above that the problem on inverse Carleson measures for
${H^{p}}(B_{d})$, $1<p<\infty$, has no solutions for the measures $\mu$
defined on the open unit ball $B_{d}$. Indeed, this fact follows immediately
from property (iii) in Theorem 1. Further comments and the proof of Theorem 1
are given in Section 2.

###  1.4. Reverse Carleson measures for de Branges–Rovnyak spaces

The second type of spaces studied in this paper is that of de Branges–Rovnyak
spaces in the unit ball.

Let $b:B_{d}\to\mathbb{D}$ be a holomorphic function. Direct inspection shows
that the function

$\mathcal{K}^{b}(z,w)=\frac{1-b(z)\overline{b(w)}}{(1-\langle
z,w\rangle)^{d}},\quad z,w\in B_{d},$

has the properties of a reproducing kernel. The reproducing kernel Hilbert
space $\mathcal{H}(b)\subset H^{2}$ with kernel $\mathcal{K}^{b}(z,w)$ is
called a de Branges–Rovnyak space. Further details about these spaces in the
unit disk are given, for example, in monograph [9].

On the one hand, if $\|b\|_{\infty}<1$, then the space $\mathcal{H}(b)$
coincides with $H^{2}$, up to an equivalent norm. On the other hand, the space
$\mathcal{H}(b)$ has an explicit form for an inner function $b$. Recall that a
holomorphic function $I:B_{d}\to\mathbb{D}$ is called inner if $|I(\zeta)|=1$
for $\sigma_{d}$-almost all $\zeta\in\partial B_{d}$. Here and in what
follows, $I(\zeta)$ denotes the limit $\lim_{r\to 1-}I(r\zeta)$, as usual. It
is known that the corresponding limit exists $\sigma_{d}$-almost everywhere.
Direct computations of the reproducing kernel for the space $H^{2}\ominus
IH^{2}$ show that $\mathcal{H}(b)=H^{2}\ominus IH^{2}$. If $I$ is an inner
function in the unit disk $\mathbb{D}$, then $H^{2}(\mathbb{D})\ominus
IH^{2}(\mathbb{D})$ is a classical model space.

A number of results on the reverse Carleson measures for the spaces
$\mathcal{H}(b)$ in the unit disk are obtained in [2]. In that paper, for
certain functions $b$, descriptions of the measures under consideration are
obtained and it is shown, in particular, that it suffices to check the reverse
Carleson inequality on the corresponding reproducing kernels. For arbitrary
spaces $\mathcal{H}(b)$ in the unit disk, it is known that such a reproducing
kernel thesis does not hold (see [5]).

In the final Section 3 of the present paper, an analog of property (ii) from
Theorem 1 is used to deduce properties of the reverse Carleson measures for
the de Branges–Rovnyak spaces $\mathcal{H}(b)$ in the unit ball. The results
in this direction are obtained for non-inner functions $b$ and are contained
in Theorem 2 and Corollary 1, see Section 3.

### Notation

As usual, $C>0$ denotes a constant whose value can change from line to line.

##  2\. Reverse Carleson measures for Hardy spaces

In this section, Theorem 1 is proved and one more description of the measures
under investigation is given in terms of the standard Carleson windows
$S_{Q}$, see Proposition 1 below.

###  2.1. Auxiliary results

For $0<h\leq 1$ and a nonisotropic ball $Q$, define the Carleson $h$-window
$S_{Q,h}$ by the following equality:

(2.1) $S_{Q,h}=\left\\{z\in B_{d}:\ 1-h\leq|z|\leq 1,\ \frac{z}{|z|}\in
Q\right\\}.$

Let $\nu=\nu_{d}$ denote Lebesgue measures on the complex Euclidean space
$\mathbb{C}^{d}$.

######  Lemma 1.

Let $Q\subset\partial B_{d}$ be a nonisotropic ball. For $0<h\leq 1$, put

(2.2)
$\Phi_{h}(z)=\frac{1}{h}\int_{S_{Q,h}}\frac{(1-|w|^{2})^{pd-d}}{|1-\langle
z,w\rangle|^{pd}}\,d\nu(w),\quad z\in\overline{B}_{d}.$

Then

  1. (a)

$\Phi_{h}(z)\leq C$ for all $z\in\overline{B}_{d}$ and $0<h\leq 1;$

  2. (b)

$\lim_{h\to 0}\Phi_{h}(z)=0$ for all $z\in\overline{B}_{d}\setminus Q$.

###### Proof.

(a) It is clear that

$\Phi_{h}(z)\leq\frac{1}{h}\int_{1-h\leq|w|\leq
1}\frac{(1-|w|^{2})^{pd-d}}{|1-\langle z,w\rangle|^{pd}}\,d\nu(w),\quad z\in
B_{d}.$

Next, integrating in polar coordinates, we have

$\begin{split}\Phi_{h}(z)&\leq
C(d)\frac{1}{h}\int_{1-h}^{1}(1-r^{2})^{pd-d}\int_{\partial
B_{d}}\frac{d\sigma(\xi)}{|1-r\langle z,\xi\rangle|^{pd}}\,dr\\\ &\leq
C(d)\frac{1}{h}\int_{1-h}^{1}(1-r^{2})^{pd-d}(1-r|z|)^{d-pd}\,dr\\\ &\leq
C(d)\frac{1}{h}\int_{1-h}^{1}\,dr\leq C(d)\end{split}$

by [8, Proposition 1.4.10].

(b) Assume that ${z}\notin Q$. Then there exists a $\delta=\delta(z)>0$ such
that $|1-\langle z,w\rangle|\geq\delta$ for all $w\in S_{Q,h}$, $0<h\leq 1$.
Hence,

$\begin{split}0\leq\Phi_{h}(z)&\leq\frac{1}{h}\int_{S_{Q,h}}\frac{(1-|w|^{2})^{pd-d}}{\delta^{pd}}d\nu(w)\\\
&\leq C\frac{1}{h}\frac{\sigma(Q)h^{pd-d+1}}{\delta^{pd}}\to
0\quad\textrm{as}\ h\to 0+.\end{split}$

The proof of the lemma is finished. ∎

###  2.2. Proof of Theorem 1

$(\mathrm{i})\Rightarrow(\mathrm{ii})$ This implication is trivial, since
$\|K_{w}\|_{H^{p}}=1$ for $w\in B_{d}$.

$(\mathrm{ii})\Rightarrow(\mathrm{iii})$ Fix a nonisotropic ball
$Q\subset\partial B_{d}$. Integrating the estimate from condition (ii) over
the set

$\widetilde{S}_{Q,h}=\left\\{z\in B_{d}:\ 1-h\leq|z|<1,\ \frac{z}{|z|}\in
Q\right\\},\quad 0<h\leq 1,$

with respect to the measure $\nu$, we have

$\displaystyle C_{\mathrm{ii}}|Q|h$
$\displaystyle\leq\int_{\widetilde{S}_{Q,h}}\int_{\overline{B}_{d}}|K_{w}|^{p}\,d\mu\,d\nu(w)$
$\displaystyle=\int_{\overline{B}_{d}}\int_{\widetilde{S}_{Q,h}}\frac{(1-|w|^{2})^{\frac{pd}{q}}}{|1-\langle
z,w\rangle|^{pd}}\,d\nu(w)\,d\mu(z).$

The last expression does not change after replacing the set
${\widetilde{S}_{Q,h}}$ by ${S}_{Q,h}$ in the inner integral. Therefore, the
above estimate has the following form in terms of the function $\Phi_{h}$
defined by identity (2.2):

(2.3) $\int_{\overline{B}_{d}}\Phi_{h}(z)\,d\mu(z)\geq
C_{\mathrm{ii}}\sigma(Q).$

On the one hand, properties (a) and (b) from Lemma 1 guarantee that

(2.4) $\lim_{h\to 0+}\int_{\overline{B}_{d}\setminus
Q}\Phi_{h}(z)\,d\mu(z)=0.$

On the other hand, by property (b) from Lemma 1, we have

(2.5) $\int_{Q}\Phi_{h}(z)\,d\mu(z)\leq C\mu(Q)\quad\textrm{for all}\ h>0.$

Combining properties (2.3)–(2.5), we obtain

$C\mu(Q)\geq C_{\mathrm{ii}}\sigma(Q),$

as required.

$(\mathrm{iii})\Rightarrow(\mathrm{i})$ Let

$\mu_{|\partial B_{d}}=g\sigma+\mu^{s}_{|\partial B_{d}},$

be the Lebesgue decomposition, where $\mu^{s}_{|\partial B_{d}}\bot\sigma$ and
$g\in L^{1}(\sigma)$ denotes the corresponding Radon–Nikodym derivative.
Theorems 5.3.1 and 5.3.2 from monograph [8] guarantee that

$g(\zeta)=\lim_{\delta\to 0+}\frac{\mu_{|\partial
B_{d}}(Q(\zeta,\delta))}{\sigma(Q(\zeta,\delta))}$

for $\sigma$-almost all $\zeta\in\partial B_{d}$. Thus, condition (iii)
implies that $g\geq C>0$ $\sigma$-almost everywhere. Therefore,

$\int_{\overline{B}_{d}}|f|^{p}\,d\mu\geq\int_{\partial
B_{d}}|f|^{p}\,d\mu\geq C\int_{\partial
B_{d}}|f|^{p}\,d\sigma=C\|f\|_{H^{p}}^{p}$

for all $f\in H^{p}\cap C(\overline{B}_{d})$. The proof of Theorem 1 is
finished.

###  2.3. Comments and remarks

######  Remark 1.

Equivalence of properties (i) and (ii) in Theorem 1 means that to prove the
reverse Carleson inequality, it suffices to verify this inequality on the
reproducing kernels $k_{w}$, $w\in B_{d}$. In other words, the reproducing
kernel thesis holds for the problem on reverse Carleson measures for
${H^{p}}$, $1<p<\infty$.

######  Remark 2.

On the one hand, property (iii) guarantees that Lebesgue measure $\sigma$,
defined on the sphere $\partial B_{d}$, is absolutely continuous with respect
to the restriction $\mu_{|\partial B_{d}}$. On the other hand, as it is shown
in the proof of Theorem 1, the Radon–Nikodym derivative $d\mu^{a}_{|\partial
B_{d}}/d\sigma$ is separated from zero by a positive constant. Here
$\mu^{a}_{|\partial B_{d}}$ denotes the absolutely continuous part (with
respect to $\sigma$) of $\mu_{|\partial B_{d}}$.

The following assertion shows that condition (iii) from Theorem 1 is
equivalent to its analog in terms of the sets $S_{Q}$ used to characterize the
classical Carleson measures.

######  Proposition 1.

Equivalent properties from Theorem 1 hold if and only if there exists a
constant $C>0$ such that

(2.6) $\mu(S_{Q})\geq C\sigma(Q)\quad\textrm{for all nonisotropic balls}\
Q\subset\partial B_{d}.$

###### Proof.

Clearly, property (iii) from Theorem 1 implies (2.6). To prove the reverse
implication, assume that property (2.6) holds. Fix a $\delta>0$ and a
nonisotropic ball $Q=Q(\zeta,\delta)\subset\partial B_{d}$ for some point
$\zeta\in\partial B_{d}$. For $0<h<\delta$, choose a maximal family of
pairwise disjoint nonisotropic balls $Q_{j}=Q_{j}(\zeta_{j},h)\subset 2Q$,
$j=1,\dots,J$. Observe that $Q\subset\cup_{1\leq j\leq J}2Q_{j}$. Since the
$h$-window $S_{Q_{j},h}$ coincides with the set $S_{Q_{j}}$ for $j=1,\dots,J$,
condition (iii) guarantees that

(2.7)
$\begin{split}C_{0}(d)\mu(S_{Q,h})=\mu(S_{2Q,h})&\geq\mu\sum_{j=1}^{J}S_{Q_{j},h}=\sum_{j=1}^{J}\mu(S_{Q_{j},h})\\\
&\overset{\mathrm{(iii)}}{\geq}C\sum_{j=1}^{J}\sigma(Q_{j})=C(d)\sum_{j=1}^{J}\sigma(2Q_{j})\\\
&\geq C(d)\sigma(Q).\end{split}$

If $Q\subset G\subset\overline{B}_{d}$ and the set $G$ is relatively open in
$\overline{B}_{d}$, then $G\supset S_{Q,h}$ for a sufficiently small parameter
$h>0$. Thus, applying (2.7) and the outer regularity of the measure $\mu\in
M_{+}(\overline{B}_{d})$, we have

$\begin{split}\mu(Q)&=\inf\\{\mu(G):Q\subset G,\ G\ \textrm{is open in}\
\overline{B}_{d}\\}\\\ &\geq\inf_{h>0}\mu(S_{Q,h})\geq
C(d)\sigma(Q),\end{split}$

as required. ∎

##  3\. Reverse Carleson measures for de Branges–Rovnyak spaces

To define a reverse Carleson measure $\mu$ for $H^{p}(B_{d})$, we compute the
norm $\|f\|_{L^{p}(\overline{B}_{d},\mu)}$ only for $f\in H^{p}(B_{d})\cap
C(\overline{B}_{d})$, see (1.5). For the spaces $\mathcal{H}(b)$, we apply a
different approach based on the following notion.

######  Definition 1.

Let $\mu\in M_{+}(\overline{B}_{d})$. A function $f\in\mathcal{H}ol(B_{d})$ is
called $\mu$-admissible if $f$ has radial limits $\mu$-almost everywhere on
the sphere $\partial B_{d}$. Given a holomorphic function
$b:B_{d}\to\mathbb{D}$, the set of all $\mu$-admissible functions in the space
$\mathcal{H}(b)$ is denoted by $\mathcal{H}_{\mu}(b)$.

Now, assume that $f\in\mathcal{H}_{\mu}(b)$. On the support of the measure
$\mu_{\left|\partial B_{d}\right.}$, the function $f$ is defined via its
radial limits. As a result, the quantity $\|f\|_{L^{2}(\overline{B}_{d},\mu)}$
takes value in the infinite interval $[0,+\infty]$.

######  Definition 2.

Let $b:B_{d}\to\mathbb{D}$ be a holomorphic function. A measure $\mu\in
M_{+}(\overline{B}_{d})$ is called a reverse Carleson measure for the space
$\mathcal{H}(b)$ if

$\|f\|_{\mathcal{H}(b)}\leq C||f\|_{L^{2}(\overline{B}_{d},\mu)}$

for all $f\in\mathcal{H}_{\mu}(b)$.

The main result of the present section is the following theorem on properties
of the reverse Carleson measures. For $d=1$, this result was obtained in [2].

######  Theorem 2.

Let $\mu\in M_{+}(\overline{B}_{d})$ and $b:B_{d}\to\mathbb{D}$ be a
$\mu$-admissible holomorphic function. Put

$g=\frac{d\mu_{\left|\partial B_{d}\right.}}{dm}.$

Assume that $\mu$ is a reverse Carleson measure for the space
$\mathcal{H}(b)$. Then

(3.1) $1-|b(\zeta)|^{2}\leq C(1-|b(\zeta)|^{2})^{2}g(\zeta)$

for $\sigma_{d}$-almost all points $\zeta\in\partial B_{d}$ and a universal
constant $C>0$.

On the one hand, if $|b(\zeta)|=1$, then condition (3.1) imposes no
restriction at the point $\zeta$. In particular, this condition degenerates
completely for an inner function $b$. On the other hand, condition (3.1)
naturally simplifies for $|b(\zeta)|<1$ and provides the following explicit
information.

######  Corollary 1.

Let $\mu\in M_{+}(\overline{B}_{d})$ and $b:B_{d}\to\mathbb{D}$ be a
$\mu$-admissible non-inner holomorphic function. Assume that $\mu$ is a
reverse Carleson measure for the space $\mathcal{H}(b)$. Then

$\frac{d\mu_{\left|\partial B_{d}\right.}}{d\sigma}\neq\mathbf{0}$

and

(3.2) $\int_{\\{\zeta\in\partial
B_{d}:\,|b(\zeta)|<1\\}}\frac{1}{1-|b|}\,d\sigma<\infty.$

###### Proof.

Firstly, by Theorem 2, inequality (3.1) holds, therefore, $g\neq\mathbf{0}$
because $b$ is not an inner function. Secondly, property (3.1) implies that

$\frac{1}{1-|b(\xi)|}\leq Cg(\xi)$

for $\sigma$-almost all $\xi\in\\{\zeta\in\partial B_{d}:\,|b(\zeta)|<1\\}$.
Since $g\in L^{1}(\partial B_{d},\sigma)$, we obtain the required property
(3.2). ∎

Remark that property (3.2) imposes a restriction on those functions $b$ for
which the set of the reverse Carleson measures for $\mathcal{H}(b)$ is not
empty.

######  Corollary 2.

There exists a holomorphic function $b:B_{d}\to\mathbb{D}$ with the following
property: if $\mu\in M_{+}(\overline{B}_{d})$ and $b$ is $\mu$-admissible,
then $\mu$ is not a reverse Carleson measure for the space $\mathcal{H}(b)$.

###### Proof.

Choose a positive function $\psi\in C(\partial B_{d})$ such that
$\psi(\zeta_{0})=1$, $0<\psi(\zeta)<1$ for $\zeta\in\partial
B_{d}\setminus\\{\zeta_{0}\\}$ and

$\int_{\partial B_{d}}\frac{1}{1-\psi}\,d\sigma=\infty.$

By [1, Corollary 2], there exists a holomorphic function
$b:B_{d}\to\mathbb{D}$ such that $|b|=\psi$ $\sigma$-almost everywhere on the
sphere $\partial B_{d}$. Therefore, if $\mu\in M_{+}(\overline{B}_{d})$ and
$b$ is $\mu$-admissible, then condition (3.2) does not hold. Hence, $\mu$ is
not a reverse Carleson measure for $\mathcal{H}(b)$ by Corollary 1. ∎

Further, Corollary 1 provides negative information on the sampling sequences
for $\mathcal{H}(b)$. Recall that $\\{w_{j}\\}_{j=1}^{\infty}\subset B_{d}$ is
called a sampling sequence for $\mathcal{H}(b)$ provided that there exist
constants $C_{2}\geq C_{1}>0$ such that

$C_{1}\|f\|^{2}_{\mathcal{H}(b)}\leq\sum_{j=1}^{\infty}\frac{|f(w_{j})|^{2}}{\
\|\mathcal{K}_{w_{j}}\|^{2}_{\mathcal{H}(b)}}\leq
C_{2}\|f\|^{2}_{\mathcal{H}(b)}$

for all $f\in\mathcal{H}(b)$. For $d=1$, the following assertion is proved in
[2].

######  Corollary 3.

Let $b:B_{d}\to\mathbb{D}$ be a non-inner holomorphic function. Then there is
no sampling sequence for the space $\mathcal{H}(b)$.

###### Proof.

Assume that $\\{w_{j}\\}_{j=1}^{\infty}\subset B_{d}$ is a sampling sequence
for $\mathcal{H}(b)$. Then

$\mu=\sum_{j=1}^{\infty}\frac{\delta_{w_{j}}}{\|\mathcal{K}_{w_{j}}\|^{2}_{\mathcal{H}(b)}}$

is a reverse Carleson measure for $\mathcal{H}(b)$. This contradicts Corollary
1 because $d\mu_{|\partial B_{d}}/d\sigma=\mathbf{0}$. ∎

Now, we proceed to the proof of Theorem 2.

###  3.1. Auxiliary results

Given a holomorphic function $F:B_{d}\to\mathbb{D}$, recall that the boundary
function $F(\xi)$ is defined for $\sigma$-almost all $\xi\in\partial B_{d}$ in
terms of radial limits. The following assertion refines Lemma 1 for $p=2$.

######  Lemma 2.

Let $Q\subset\partial B_{d}$ be a nonisotropic ball and $F:B_{d}\to\mathbb{D}$
be a holomorphic function. Then the limit

$\lim_{r\to 1-}\int_{Q}|F(r\xi)|\frac{(1-r^{2})^{d}}{|1-\langle
z,r\xi\rangle|^{2d}}d\sigma(\xi)$

  1. (a)

is equal to zero for $z\in\overline{B}_{d}\setminus Q$,

  2. (b)

is equal $|F(z)|$ for $\sigma$-almost all $z$ in the interior of $Q$.

###### Proof.

(a) Since $F$ is a bounded function, it suffices to apply Lemma 1(b) for
$p=2$.

(b) Let $\zeta\in\partial B_{d}$. Observe that

$\lim_{r\to 1-}\int_{\partial
B_{d}\setminus\widetilde{Q}(\zeta,\tau)}\frac{(1-r^{2})^{d}}{|1-\langle\zeta,r\xi\rangle|^{2d}}d\sigma(\xi)=0,$

where $\widetilde{Q}(\zeta,\tau)\subset\partial B_{d}$ is a nonisotropic ball.
Thus, to prove property (b), it suffices to verify that

(3.3) $\lim_{r\to 1-}\int_{\partial
B_{d}}|F(r\xi)|\frac{(1-r^{2})^{d}}{|1-\langle\zeta,r\xi\rangle|^{2d}}d\sigma(\xi)=|F(\zeta)|$

for $\sigma$-almost all points $\zeta\in\partial B_{d}$.

Now, to prove property (3.3), assume that $\zeta\in\partial B_{d}$ is a
Lebesgue point for the function $F(\xi)$, i.e.,

(3.4) $\lim_{\tau\to
0+}\int_{\widetilde{Q}(\zeta,\tau)}|F(\xi)-F(\zeta)|\,d\sigma(\xi)=0.$

Consider the following Poisson integral:

$u(w)=\int_{\partial B_{d}}|F(\xi)-F(\zeta)|\frac{(1-|w|^{2})^{d}}{|1-\langle
w,\xi\rangle|^{2d}}d\sigma(\xi),\quad w\in B_{d}.$

Remark that $\lim_{r\to 1-}u(r\zeta)=0$ by (3.4) and Theorem 5.4.8 from [8].

The function $|F(w)-F(\zeta)|$, $w\in B_{d}$, is plurisubharmonic, hence,

$\begin{split}0&\leq\limsup_{r\to 1-}\int_{\partial
B_{d}}|F(r\xi)-F(\zeta)|\frac{(1-r^{2})^{d}}{|1-\langle\zeta,r\xi\rangle|^{2d}}\,d\sigma(\xi)\\\
&\leq\limsup_{r\to 1-}\int_{\partial
B_{d}}u(r\xi)\frac{(1-r^{2})^{d}}{|1-\langle
r\zeta,\xi\rangle|^{2d}}\,d\sigma(\xi)=\limsup_{r\to
1-}u(r^{2}\zeta)=0.\end{split}$

One has $\left||F(r\xi)|-|F(\zeta)|\right|\leq|F(r\xi)-F(\zeta)|$, thus, the
above estimate implies the required property (3.3) at the point
$\zeta\in\partial B_{d}$ under consideration. To finish the argument, recall
that property (3.4) holds for $\sigma$-almost all $\zeta\in\partial B_{d}$. ∎

###  3.2. Proof of Theorem 2

The function $\mathcal{K}^{b}_{w}$, $w\in B_{d}$, is $\mu$-admissible, since
$b$ is $\mu$-admissible by the hypothesis of the theorem. Hence, by the
definition of a reverse Carleson measure, we have

$\|\mathcal{K}^{b}_{w}\|^{2}_{\mathcal{H}(b)}\leq
C||\mathcal{K}^{b}_{w}\|^{2}_{L^{2}(\mu)}$

for all $w\in B_{d}$. Rewriting this inequality in explicit terms, we obtain

(3.5) $\frac{1-|b(w)|^{2}}{(1-|w|^{2})^{d}}\leq
C\int_{\overline{B}_{d}}\left|\frac{1-b(z)\overline{b(w)}}{(1-\langle
z,w\rangle)^{d}}\right|^{2}d\mu(z).$

Let

$\mu_{\left|\partial B_{d}\right.}=g\sigma_{d}+\eta$

be the Lebesgue decomposition, where $\eta\bot\mu_{\left|\partial
B_{d}\right.}$ and $g$ denotes the Radon–Nikodym derivative
$\frac{d\mu_{\left|\partial B_{d}\right.}}{d\sigma_{d}}$. Put

$\rho:=\mu_{\left|B_{d}\right.}=\mu-g\sigma-\eta.$

In these terms, estimate (3.5) has the following form:

(3.6)
$\begin{split}C(1-|b(w)|^{2})&\leq\int_{B_{d}}\frac{(1-|w|^{2})^{d}}{|1-\langle
z,w\rangle|^{2d}}|1-b(z)\overline{b(w)}|^{2}\,d\rho(z)\\\
&\quad+\int_{\partial B_{d}}\frac{(1-|w|^{2})^{d}}{|1-\langle
z,w\rangle|^{2d}}|1-b(z)\overline{b(w)}|^{2}\,d\eta(z)\\\
&\quad+\int_{\partial B_{d}}\frac{(1-|w|^{2})^{d}}{|1-\langle
z,w\rangle|^{2d}}|1-b(z)\overline{b(w)}|^{2}\,g(z)\,d\sigma_{d}(z)\\\
&\quad:=J_{1}+J_{2}+J_{3}.\end{split}$

Now, choose a point $\zeta\in\partial B_{d}$ and consider a non-isotropic ball
$Q=(\zeta,\tau)\subset\partial B_{d}$, $\tau>0$. Recall that the Carleson
$h$-window $S_{Q,h}$ for $0<h\leq 1$ is defined by the following equality:

$S_{Q,h}=\left\\{z\in B_{d}:\ 1-h\leq|z|\leq 1,\ \frac{z}{|z|}\in Q\right\\}.$

We integrate both parts of inequality (3.6) over the set $S_{Q,h}$ with
respect to the volume Lebesgue measure $\nu$, divide by $h$ and study the
behavior of the resulting expressions as $h\to 0+$.

Firstly, for the integral of the left-hand side of inequality (3.6), after
integration in polar coordinates we have

(3.7) $\begin{split}C\int_{S_{Q,h}}&(1-|b(w)|^{2})\,d\nu(w)\\\
&=\frac{C(d)}{h}\int_{1-h}^{1}\int_{Q}(1-|b(r\xi)|^{2})\,d\sigma(\xi)\,r^{2d-1}dr\\\
&\to C(d)\int_{Q}(1-|b(\xi)|^{2})\,d\sigma(\xi)\end{split}$

as $h\to 0$.

Secondly, applying Fubini’s theorem, we obtain

$\begin{split}&\frac{C(d)}{h}\int_{S_{Q,h}}J_{1}(w)\,d\nu(w)\\\
&=\int_{B_{d}}\frac{1}{h}\int_{1-h}^{1}\int_{Q}\frac{(1-r^{2})^{d}}{|1-\langle
z,r\xi\rangle|^{2d}}|1-b(z)\overline{b(r\xi)}|^{2}\,d\sigma(\xi)\,r^{2d-1}dr\,d\rho(z).\end{split}$

Remark that

(3.8) $\begin{split}\int_{Q}\frac{(1-r^{2})^{d}}{|1-\langle
z,r\xi\rangle|^{2d}}&|1-b(z)\overline{b(r\xi)}|^{2}\,d\sigma(\xi)\\\ &\leq
C\int_{Q}\frac{(1-r^{2})^{d}}{|1-\langle
z,r\xi\rangle|^{2d}}\,d\sigma(\xi)\leq C\end{split}$

for all $z\in\overline{B}_{d}$, and

$\frac{1}{h}\int_{1-h}^{1}\int_{Q}\frac{(1-r^{2})^{d}}{|1-\langle
z,r\xi\rangle|^{2d}}\,d\sigma(\xi)\,r^{2d-1}dr\overset{h\to
0}{\longrightarrow}0\quad\textrm{for all}\ z\in B_{d}$

by Lemma 2(i) with $F\equiv 1$. Therefore,

(3.9) $\frac{1}{h}\int_{S_{Q,h}}J_{1}(w)\,d\nu(w)\overset{h\to
0}{\longrightarrow}0.$

Next,

$\frac{1}{h}\int_{1-h}^{1}\int_{Q}\frac{(1-r^{2})^{d}}{|1-\langle
z,r\xi\rangle|^{2d}}\,d\sigma(\xi)\,r^{2d-1}dr\overset{h\to
0}{\longrightarrow}0\quad\textrm{for all}\ z\in\partial B_{d}\setminus Q$

by Lemma 2(i) with $F\equiv 1$. Therefore, applying Fubini’s theorem and
property (3.8), we obtain

(3.10) $\begin{split}&\frac{C(d)}{h}\int_{S_{Q,h}}J_{2}(w)\,d\nu(w)\\\
&=\int_{\partial
B_{d}}\frac{1}{h}\int_{1-h}^{1}\int_{Q}\frac{(1-r^{2})^{d}}{|1-\langle
z,r\xi\rangle|^{2d}}|1-b(z)\overline{b(r\xi)}|^{2}\,d\sigma(\xi)\,r^{2d-1}dr\,d\eta(z)\\\
&\leq\int_{Q}d\eta+\varepsilon(h),\end{split}$

where $\varepsilon(h)\to 0$ as $h\to 0$.

Finally, after application of Fubini’s theorem, we have

$\begin{split}&\frac{C(d)}{h}\int_{S_{Q,h}}J_{3}(w)\,d\nu(w)\\\
&=\int_{\partial
B_{d}}\frac{1}{h}\int_{1-h}^{1}\int_{Q}\frac{(1-r^{2})^{d}}{|1-\langle
z,r\xi\rangle|^{2d}}|1-b(r\xi)\overline{b(z)}|^{2}d\sigma(\xi)r^{2d-1}drg(\zeta)d\sigma(z).\end{split}$

Applying Lemma 2 with $F(w)=(1-b(w)\overline{b(\zeta)})^{2}$, we obtain

$C\geq\int_{Q}\frac{(1-r^{2})^{d}}{|1-\langle
z,r\xi\rangle|^{2d}}|1-b(r\xi)\overline{b(z)}|^{2}\,d\sigma(\xi)\overset{r\to
1-}{\longrightarrow}(1-|b(z)|^{2})^{2}\chi_{Q}(z)$

for $\sigma$-almost all $z\in\partial B_{d}$. Therefore,

(3.11)
$\frac{1}{h}\int_{S(Q,h)}J_{3}(w)\,d\nu(w)\to\int_{Q}(1-|b|^{2})^{2}g\,d\sigma$

as $h\to 0$.

Recall that $Q=Q(\zeta,\tau)$. So, combining properties (3.7) and
(3.9)–(3.11), we obtain

$\int_{Q(\zeta,\tau)}(1-|b|^{2})\,d\sigma\leq
C\int_{Q(\zeta,\tau)}(1-|b|^{2})^{2}g\,d\sigma+C\int_{Q(\zeta,\tau)}d\eta$

for all $\zeta\in\partial B_{d}$.

Now, we divide both parts of the above estimate by $\sigma(Q(\zeta,\tau))$ and
let $\tau$ tend to zero. Since $\eta$ is a singular measure, Theorem 5.3.2
from monograph [8] guarantees that

$\frac{1}{\sigma(Q(\zeta,\tau))}\int_{Q(\zeta,\tau)}d\eta\overset{\tau\to
0}{\longrightarrow}0$

for $\sigma$-almost all $\zeta\in\partial B_{d}$. Therefore, the required
inequality (3.1) holds for $\sigma$-almost all $\zeta\in\partial B_{d}$ by the
measure differentiation theorem on the sphere $\partial B_{d}$ (see [8,
Theorem 5.3.1]). The proof of Theorem 2 is finished.

## References

  * [1] A. B. Aleksandrov, _The existence of inner functions in the ball_ , Mat. Sb. (N.S.) 118(160) (1982), no. 2, 147–163 (Russian); English transl.: Math. USSR–Sb. 46 (1983), no. 2, 143–159.
  * [2] A. Blandignères, E. Fricain, F. Gaunard, A. Hartmann, and W. T. Ross, _Direct and reverse Carleson measures for $\mathcal{H}(b)$ spaces_, Indiana Univ. Math. J. 64 (2015), no. 4, 1027–1057.
  * [3] L. Carleson, _An interpolation problem for bounded analytic functions_ , Amer. J. Math. 80 (1958), 921–930.
  * [4] A. W. Green and N. A. Wagner, _Dominating sets in Bergman spaces on strongly pseudoconvex domains_ , Constr. Approx. (2023), https://doi.org/10.1007/s00365–023–09639–z.
  * [5] A. Hartmann, X. Massaneda, A. Nicolau, and J. Ortega-Cerdà, _Reverse Carleson measures in Hardy spaces_ , Collect. Math. 65 (2014), no. 3, 357–365.
  * [6] P. Lefèvre, D. Li, H. Queffélec, and L. Rodríguez-Piazza, _Some revisited results about composition operators on Hardy spaces_ , Rev. Mat. Iberoam. 28 (2012), no. 1, 57–76.
  * [7] D. H. Luecking, _Forward and reverse Carleson inequalities for functions in Bergman spaces and their derivatives_ , Amer. J. Math. 107 (1985), no. 1, 85–111.
  * [8] W. Rudin, _Function theory in the unit ball of ${\bf C}^{n}$_, Grundlehren der Mathematischen Wissenschaften, vol. 241, Springer-Verlag, New York-Berlin, 1980.
  * [9] D. Sarason, _Sub-Hardy Hilbert spaces in the unit disk_ , University of Arkansas Lecture Notes in the Mathematical Sciences, vol. 10, John Wiley & Sons, Inc., New York, 1994, A Wiley-Interscience Publication.
  * [10] K. Zhu, _Spaces of holomorphic functions in the unit ball_ , Graduate Texts in Mathematics, vol. 226, Springer-Verlag, New York, 2005.


# Double parton distributions with flavor interference from lattice QCD

Daniel Reitinger b Christian Zimmermann c Markus Diehl a and Andreas Schäfer
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Institute for Theoretical Physics,
University of Regensburg, 93040 Regensburg, GermanyAix Marseille Univ,
Université de Toulon, CNRS, CPT, Marseille, FranceDeutsches Elektronen-
Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany

###### Abstract

We study double parton distributions with flavor interference in the nucleon
and compare them with previous results for the flavor diagonal case. We
investigate both unpolarized and polarized partons, and we compare our lattice
results with those obtained from the simple description of the proton in terms
of an SU(6) symmetric three-quark wave function.

††preprint:  DESY-24-014

##  1 Introduction

The Standard Model of particle physics has been extremely successful in
describing experimental data at high energies. As part of this, the
fundamental QCD Lagrangian has been established beyond any reasonable doubt.
However, many aspects of the resulting hadron properties are still only poorly
understood, e.g. the entanglement of two partons in the wave functions of a
proton. Work on this front is not only motivated by the goal to better
understand the strong interaction, but also by the fact that QCD effects often
limit the sensitivity of searches for physics beyond the Standard Model.

An intriguing phenomenon in hadron collisions is double parton scattering
(DPS), a mechanism in which two partons in each hadron take part in a hard
scattering subprocess. Building on pioneering work from the 1970s and 1980s
Landshoff:1978fq ; Kirschner:1979im ; Politzer:1980me ; Paver:1982yp ;
Shelest:1982dg ; Mekhfi:1983az ; Sjostrand:1986ep , substantial progress has
been made during the last decade in an effort to develop a systematic
description of DPS in QCD Blok:2010ge ; Gaunt:2011xd ; Ryskin:2011kk ;
Blok:2011bu ; Diehl:2011yj ; Manohar:2012jr ; Manohar:2012pe ; Ryskin:2012qx ;
Gaunt:2012dd ; Blok:2013bpa ; Diehl:2017kgu ; Cabouat:2019gtm ;
Cabouat:2020ssr . Experimental investigations of DPS started in the 1980s
AxialFieldSpectrometer:1986dfj  and were followed by a wealth of studies at
the Tevatron and the LHC, see for instance CDF:1993sbj ; CDF:1997lmq ;
D0:2009apj ; D0:2015dyx  and LHCb:2012aiv ; ATLAS:2013aph ; CMS:2013huw ;
ATLAS:2016rnd ; LHCb:2015wvu ; LHCb:2020jse ; CMS:2022pio .

A crucial input for computing DPS are double parton distributions (DPDs),
which describe the joint distribution of two partons inside a hadron. They
quantify several types of two-particle correlations in the proton wave
function and are not well known. It is natural to explore to which extent
lattice QCD calculations can provide guidance in this context. In two previous
publications Bali:2020mij ; Bali:2021gel , we presented lattice computations
of two-current correlation functions that can be related with the Mellin
moments of DPDs Diehl:2011yj . This generalizes the well-known relation
between single-current matrix elements and the Mellin moments of single-parton
distributions, which has been extensively studied in the literature
Lin:2017snn . More recently, it has been proposed in Zhang:2023wea ;
Jaarsma:2023woo  to study the full functional dependence of DPDs on the
lattice in the LaMET approach Ji:2013dva ; Ji:2020ect . It will be very
interesting to see to which extent this can be done in practice.

One should bear in mind that neither of these lattice approaches is sensitive
to parton momentum fractions smaller than, say, $10^{-2}$, which are
responsible for much of the phase space where DPS is observed. (The same holds
for quark models, which have been used extensively to compute DPDs
Chang:2012nw ; Rinaldi:2013vpa ; Broniowski:2013xba ; Rinaldi:2014ddl ;
Broniowski:2016trx ; Kasemets:2016nio ; Rinaldi:2016jvu ; Rinaldi:2016mlk ;
Rinaldi:2018zng ; Courtoy:2019cxq ; Broniowski:2019rmu ; Broniowski:2020jwk .)
However, partons with larger momentum fractions at low resolution scales are
the “seeds” of evolution and thus have an imprint on partons with smaller
momentum fractions at high scales. At high scales, large parton momentum
fractions are probed in the production of heavy particles and thus of interest
in searches for new physics. Moreover, larger parton momentum fractions become
relevant when the products of the two hard scatters in DPS have a large
rapidity separation. This is an interesting kinematic region, where DPS is
often appreciable compared to single hard scattering. Finally, DPDs at larger
momentum fractions are of interest in their own right from the point of view
of exploring hadron structure.

With this in mind, we complement in the present paper our previous study
Bali:2021gel  of DPDs in the nucleon. Specifically, we investigate flavor
interference DPDs, which are characterized by different flavors for the quark
(or antiquark) initiating a specific hard scattering in the amplitude and in
its complex conjugate. Such distributions contribute for instance to the
double Drell-Yan process. They were introduced in Diehl:2011yj , were it was
also pointed out that they do not mix with gluons under evolution, such that
at small momentum fractions one can expect them to be small compared with
flavor diagonal DPDs. At moderate or large momentum fractions (which are
relevant for Mellin moments) there is however no argument that flavor
interference should be suppressed. We will investigate whether this is the
case by comparing the corresponding two-current matrix elements computed on
the lattice. We will furthermore compare our lattice results (both for the
flavor diagonal and the interference case) with the predictions obtained from
an $SU(6)$ symmetric three-quark wave function of the proton.

Our paper is organized as follows. In section 2, we review the different
quantities relevant to our study and explain how they are related to each
other. Details of the lattice setup we use are given in 3. Sections 4 and 5
contain the results of our calculations. In section 4, the Mellin moments for
different combinations of flavor and polarization are presented and compared
with the $SU(6)$ predictions. In section 5 we explore to which extent DPDs can
be factorized in terms of single-parton distributions. We summarize our
findings in section 6.

##  2 Theory background

###  2.1 Definitions and properties

In the following, we review some basic definitions and properties of double
parton distributions in the context of our lattice simulation. For more
details the reader is referred to Bali:2020mij ; Bali:2021gel . DPDs describe
the joint probability of finding two quarks with given polarization in a
hadron. In this work we focus on the proton. We average over its polarization
$\lambda$, which is indicated by the notation
$\sum^{\prime}_{\lambda}=\frac{1}{2}\sum_{\lambda}$. The definition of DPDs is
given by

$\displaystyle
F_{a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{y})=2p^{+}\int\mathrm{d}y^{-}\int\frac{\mathrm{d}z^{-}_{1}}{2\pi}\,\frac{\mathrm{d}z^{-}_{2}}{2\pi}\,$
$\displaystyle e^{i\mskip 1.5mu(x_{1}z_{1}^{-}+x_{2}z_{2}^{-})\mskip
1.5mup^{+}}$
$\displaystyle\times\sideset{}{{}^{\prime}}{\sum}_{\lambda}\bra{p,\lambda}\mathcal{O}_{a_{1}}(y,z_{1})\,\mathcal{O}_{a_{2}}(0,z_{2})\ket{p,\lambda}\,,$
(1)

where we use light-cone coordinates $v^{\pm}:=(v^{0}\pm v^{3})/\sqrt{2}$ and
$\boldsymbol{v}:=(v^{1},v^{2})$ for a given four-vector $v^{\mu}$. The light-
cone operators are defined as

$\displaystyle\mathcal{O}_{a}(y,z)$
$\displaystyle=\bar{q}\bigl{(}y-{\textstyle\frac{1}{2}}z\bigr{)}\,\Gamma_{a}\,q\bigl{(}y+{\textstyle\frac{1}{2}}z\bigr{)}\Big{|}_{z^{+}=y^{+}=0,\,\boldsymbol{z}=\boldsymbol{0}}\,,$
(2)

where $a$ specifies the quark flavor and polarization, which is determined by
the spin projections

$\displaystyle\Gamma_{q}$ $\displaystyle={\textstyle\frac{1}{2}}\gamma^{+}\,,$
$\displaystyle\Gamma_{\Delta q}$
$\displaystyle={\textstyle\frac{1}{2}}\gamma^{+}\gamma_{5}\,,$
$\displaystyle\Gamma_{\delta
q}^{j}={\textstyle\frac{1}{2}}i\sigma^{j+}\gamma_{5}\quad(j=1,2)\,.$ (3)

$q$ refers to an unpolarized quark, $\Delta q$ to a longitudinally polarized
quark and $\delta q$ to a transversely polarized quark. The DPDs (2.1) can be
decomposed in terms of rotationally invariant functions:

$\displaystyle F_{q_{1}q_{2}}(x_{1},x_{2},\boldsymbol{y})$
$\displaystyle=f_{q_{1}q_{2}}(x_{1},x_{2},y^{2})\,,$ $\displaystyle F_{\Delta
q_{1}\Delta q_{2}}(x_{1},x_{2},\boldsymbol{y})$ $\displaystyle=f_{\Delta
q_{1}\Delta q_{2}}(x_{1},x_{2},y^{2})\,,$ $\displaystyle F_{\delta
q_{1}q_{2}}^{j_{1}}(x_{1},x_{2},\boldsymbol{y})$
$\displaystyle=\epsilon^{j_{1}k}\boldsymbol{y}^{k}\,mf_{\delta
q_{1}q_{2}}(x_{1},x_{2},y^{2})\,,$ $\displaystyle F_{q_{1}\delta
q_{2}}^{j_{2}}(x_{1},x_{2},\boldsymbol{y})$
$\displaystyle=\epsilon^{j_{2}k}\boldsymbol{y}^{k}\,mf_{q_{1}\delta
q_{2}}(x_{1},x_{2},y^{2})\,,$ $\displaystyle F_{\delta q_{1}\delta
q_{2}}^{j_{1}j_{2}}(x_{1},x_{2},\boldsymbol{y})$
$\displaystyle=\delta^{j_{1}j_{2}}f_{\delta q_{1}\delta
q_{2}}(x_{1},x_{2},y^{2})$
$\displaystyle\quad+\bigl{(}2\boldsymbol{y}^{j_{1}}\boldsymbol{y}^{j_{2}}-\delta^{j_{1}j_{2}}\boldsymbol{y}^{2}\bigr{)}\mskip
1.5mum^{2}f^{\mskip 1.5mut}_{\delta q_{1}\delta q_{2}}(x_{1},x_{2},y^{2})\,,$
(4)

where $m$ denotes the proton mass and $\epsilon^{jk}$ is the antisymmetric
tensor in two dimensions ($\epsilon^{12}\leavevmode\nobreak\
=\leavevmode\nobreak\ 1$). The definitions given above can be extended for the
case of flavor-changing operators, i.e. (2) is modified so that the quark
field has a different quark flavor than the conjugate quark field, $q\neq
q^{\prime}$:

$\displaystyle\mathcal{O}_{a}(y,z)$
$\displaystyle=\bar{q}\bigl{(}y-{\textstyle\frac{1}{2}}z\bigr{)}\,\Gamma_{a}\,q^{\prime}\bigl{(}y+{\textstyle\frac{1}{2}}z\bigr{)}\Big{|}_{z^{+}=y^{+}=0,\,\boldsymbol{z}=\boldsymbol{0}}\,,$
(5)

with $a=(qq^{\prime}),\Delta(qq^{\prime}),\delta(qq^{\prime})$. The
corresponding functions $F_{a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{y})$ given by
inserting the operators (5) in (2.1) are called flavor interference
distributions or flavor interference DPDs. For instance, an interference
contribution to the unpolarized channel is given by $F_{(ud)(du)}$, which
corresponds to the operator combination
$\mathcal{O}_{(ud)}\mathcal{O}_{(du)}$. Notice that in contrast to flavor
diagonal DPDs, these interference DPDs cannot be interpreted as parton density
distributions. Hence, there is no positivity constraint for them. Flavor
interference DPDs are relevant in the description of double parton scattering,
where they represent flavor interference contributions in the cross section.
An example is shown in figure 1.

Figure 1: Example graphs for flavor interference in the double Drell-Yan
process.

In the context of lattice calculations, it is useful to introduce so-called
skewed DPDs, which are a generalization of ordinary DPDs in the sense that one
introduces a difference between the momentum fractions of the quark in the
wave function and in its complex conjugate. This difference is quantified by
the skewness parameter $\zeta$:

$\displaystyle F_{a_{1}a_{2}}(x_{1},x_{2},\zeta,\boldsymbol{y})$
$\displaystyle=2p^{+}\int\mathrm{d}y^{-}e^{-i\zeta
y^{-}p^{+}}\int\frac{\mathrm{d}z^{-}_{1}}{2\pi}\,\frac{\mathrm{d}z^{-}_{2}}{2\pi}\,e^{i\mskip
1.5mu(x_{1}z_{1}^{-}+x_{2}z_{2}^{-})\mskip 1.5mup^{+}}$
$\displaystyle\qquad\times\sideset{}{{}^{\prime}}{\sum}_{\lambda}\bra{p,\lambda}\mathcal{O}_{a_{1}}(y,z_{1})\,\mathcal{O}_{a_{2}}(0,z_{2})\ket{p,\lambda}\,.$
(6)

This skewed DPD has already been used in Bali:2020mij ; Bali:2021gel , where
also the region of support w.r.t. the parameters $x_{1}$, $x_{2}$ and $\zeta$
has been discussed in detail. Let us recall that this support region is given
by:

$\displaystyle|x_{i}\pm{\textstyle\frac{1}{2}}\zeta|\leq
1\,,\qquad|x_{1}|+|x_{2}|\leq 1\,,\qquad|\zeta|\leq 1\,.$ (7)

In figure 2 we give a graphical representation of skewed DPDs, where we
indicate the longitudinal momentum fraction of the quarks in the wave function
and its complex conjugate.

Figure 2: Skewed DPDs and their parton momentum fractions for the flavor
diagonal case (left) and for flavor interference (right).

Moreover, we define Mellin moments w.r.t. $x_{1}$ and $x_{2}$:

$\displaystyle I_{a_{1}a_{2}}(\zeta,y^{2})$
$\displaystyle=\int_{-1}^{1}\mathrm{d}x_{1}\int_{-1}^{1}\mathrm{d}x_{2}\;f_{a_{1}a_{2}}(x_{1},x_{2},\zeta,y^{2})\,.$
(8)

The symmetry properties of flavor diagonal DPDs are discussed in Bali:2021gel
. In the following we briefly show how these properties generalize to the
flavor interference case.

For unpolarized distributions, one finds

$\displaystyle f_{(ud)(du)}(x_{1},x_{2},\zeta,y^{2})$
$\displaystyle=f_{(du)(ud)}(x_{1},x_{2},-\zeta,y^{2})\,,$
$\displaystyle\bigl{[}f_{(ud)(du)}(x_{1},x_{2},\zeta,y^{2})\bigr{]}^{*}$
$\displaystyle=f_{(du)(ud)}(x_{1},x_{2},-\zeta,y^{2})\,,$ $\displaystyle
f_{(ud)(du)}(x_{1},x_{2},\zeta,y^{2})$
$\displaystyle=f_{(du)(ud)}(x_{2},x_{1},-\zeta,y^{2})\,,$ (9)

where the first relation follows from $PT$ invariance, the second one from
taking the hermitian conjugate of the definition (2.1), and the third one from
interchanging the two operators. From the first two relations, it follows that
$f_{(ud)(du)}$ is real valued. Equations analogous to (2.1) relate
$f_{\Delta(ud)\mskip 1.5mu\Delta(du)}$ with $f_{\Delta(du)\mskip
1.5mu\Delta(ud)}$, $f_{\delta(ud)\mskip 1.5mu\delta(du)}$ with
$f_{\delta(du)\mskip 1.5mu\delta(ud)}$, and $f^{t}_{\delta(ud)\mskip
1.5mu\delta(du)}$ with $f^{t}_{\delta(du)\mskip 1.5mu\delta(ud)}$. For a
single transverse polarization, we have instead

$\displaystyle f_{(ud)\mskip 1.5mu\delta(du)}(x_{1},x_{2},\zeta,y^{2})$
$\displaystyle=f_{(du)\mskip 1.5mu\delta(ud)}(x_{1},x_{2},-\zeta,y^{2})\,,$
$\displaystyle\bigl{[}f_{(ud)\mskip
1.5mu\delta(du)}(x_{1},x_{2},\zeta,y^{2})\bigr{]}^{*}$
$\displaystyle=f_{(du)\mskip 1.5mu\delta(ud)}(x_{1},x_{2},-\zeta,y^{2})\,,$
$\displaystyle f_{(ud)\mskip 1.5mu\delta(du)}(x_{1},x_{2},\zeta,y^{2})$
$\displaystyle={}-f_{\delta(du)(ud)}(x_{2},x_{1},-\zeta,y^{2})\,.$ (10)

Also in this case, we find that the distributions are real valued. Parity
invariance implies that $F_{(ud)\mskip 1.5mu\Delta(du)}$, $F_{(du)\mskip
1.5mu\Delta(ud)}$, $F_{\Delta(ud)(du)}$, and $F_{\Delta(du)(ud)}$ are zero.

###  2.2 Euclidean matrix elements

Information about DPDs can be obtained from first principles on the lattice
through Euclidean two-current matrix elements. This has been worked out in
detail for the pion Bali:2020mij  and the nucleon Bali:2021gel . In the
following, we recall the important definitions and relations to DPDs and
extend them to the case of flavor interference. The Euclidean two-current
matrix element of the nucleon is defined as:

$\displaystyle
M^{\mu_{1}\cdots\mu_{2}\cdots}_{q_{1}q_{2}q_{3}q_{4},i_{1}i_{2}}(p,y)$
$\displaystyle:=\sideset{}{{}^{\prime}}{\sum}_{\lambda}\bra{p,\lambda}J^{\mu_{1}\cdots}_{q_{1}q_{2},i_{1}}(y)\,J^{\mu_{2}\cdots}_{q_{3}q_{4},i_{2}}(0)\ket{p,\lambda}\,,$
(11)

where we take the average of the proton spin. The currents
$J^{\mu\dots}_{qq^{\prime},i}$ are local quark bilinear operators. In this
work we focus on three types of currents, which are defined as:

$\displaystyle J_{qq^{\prime},V}^{\mu}(y)$ $\displaystyle=\bar{q}(y)\mskip
1.5mu\gamma^{\mu}\mskip 1.5muq^{\prime}(y)\,,$ $\displaystyle
J_{qq^{\prime},A}^{\mu}(y)$ $\displaystyle=\bar{q}(y)\mskip
1.5mu\gamma^{\mu}\gamma_{5}\,q^{\prime}(y)\,,$ $\displaystyle
J_{qq^{\prime},T}^{\mu\nu}(y)$ $\displaystyle=\bar{q}(y)\mskip
1.5mu\sigma^{\mu\nu}\mskip 1.5muq^{\prime}(y)\,.$ (12)

In order to make contact with the DPDs, we decompose the matrix elements in
terms of Lorentz invariant functions:

$\displaystyle
M^{\\{\mu\nu\\}}_{q_{1}q_{2}q_{3}q_{4},VV}-\tfrac{1}{4}g^{\mu\nu}g_{\alpha\beta}M^{\alpha\beta}_{q_{1}q_{2}q_{3}q_{4},VV}$
$\displaystyle=u_{VV,A}^{\mu\nu}\,A_{(q_{1}q_{2})(q_{3}q_{4})}+u_{VV,B}^{\mu\nu}\,m^{2}\mskip
1.5muB_{(q_{1}q_{2})(q_{3}q_{4})}$
$\displaystyle\quad+u_{VV,C}^{\mu\nu}\,m^{4}\mskip
1.5muC_{(q_{1}q_{2})(q_{3}q_{4})}\,,$ $\displaystyle
M^{\mu\nu\rho}_{q_{1}q_{2}q_{3}q_{4},TV}+\tfrac{2}{3}g_{\vphantom{q_{1}T}}^{\rho[\mu}M^{\nu]\alpha\beta}_{q_{1}q_{2}q_{3}q_{4},TV}g_{\alpha\beta}$
$\displaystyle=u_{TV,A}^{\mu\nu\rho}\,m\mskip
1.5muA_{\delta(q_{1}q_{2})(q_{3}q_{4})}+u_{TV,B}^{\mu\nu\rho}\,m^{3}\mskip
1.5muB_{\delta(q_{1}q_{2})(q_{3}q_{4})}\,,$
$\displaystyle\tfrac{1}{2}\,\bigl{[}M^{\mu\nu\rho\sigma}_{q_{1}q_{2}q_{3}q_{4},TT}+M^{\rho\sigma\mu\nu}_{q_{1}q_{2}q_{3}q_{4},TT}\bigr{]}$
$\displaystyle=\tilde{u}_{TT,A}^{\mu\nu\rho\sigma}\,A_{\delta(q_{1}q_{2})\delta(q_{3}q_{4})}+\tilde{u}_{TT,B}^{\mu\nu\rho\sigma}\,m^{2}\mskip
1.5muB_{\delta(q_{1}q_{2})\delta(q_{3}q_{4})}$
$\displaystyle\quad+\tilde{u}_{TT,C}^{\mu\nu\rho\sigma}\,m^{2}\mskip
1.5muC_{\delta(q_{1}q_{2})\delta(q_{3}q_{4})}+\tilde{u}_{TT,D}^{\mu\nu\rho\sigma}\,m^{4}\mskip
1.5muD_{\delta(q_{1}q_{2})\delta(q_{3}q_{4})}$
$\displaystyle\quad+u_{TT,E}^{\mu\nu\rho\sigma}\,m^{2}\mskip
1.5mu\widetilde{E}_{\delta(q_{1}q_{2})\delta(q_{3}q_{4})}\,,$ (13)

where the basis tensors $u$ and $\tilde{u}$, which depend on the Lorentz
vectors $y$ and $p$, have been defined in Bali:2021gel , equation (2.26). The
quantities $A_{(q_{1}q_{2})(q_{3}q_{4})}$, etc. are Lorentz scalar functions
depending only on $y^{2}$ and $py$. At leading twist, we only need to consider
the functions $A_{(q_{1}q_{2})(q_{3}q_{4})}$,
$A_{\Delta(q_{1}q_{2})\Delta(q_{3}q_{4})}$,
$A_{\delta(q_{1}q_{2})(q_{3}q_{4})}$, $A_{(q_{1}q_{2})\delta(q_{3}q_{4})}$,
$A_{\delta(q_{1}q_{2})\delta(q_{3}q_{4})}$, and
$B_{\delta(q_{1}q_{2})\delta(q_{3}q_{4})}$. These so-called twist-two
functions are directly related to the Mellin moments (8):

$\displaystyle I_{a_{1}a_{2}}(\zeta,y^{2})$
$\displaystyle=\int_{-\infty}^{\infty}\mathrm{d}(py)\,e^{-i\zeta
py}\,A_{a_{1}a_{2}}(py,y^{2})\,,$ (14) $\displaystyle
I^{t}_{a_{1}a_{2}}(\zeta,y^{2})$
$\displaystyle=\int_{-\infty}^{\infty}\mathrm{d}(py)\,e^{-i\zeta
py}\,B_{a_{1}a_{2}}(py,y^{2})\,,$ (15)

where (15) is only defined for transverse polarization of $a_{1}$ and $a_{2}$.
In this work, we shall restrict ourselves to $\vec{\mskip 0.5mup}\mskip
1.5mu=\vec{\mskip 0.5mu0}\mskip 1.5mu$ and, therefore, $py=0$. In this case,
the twist-two functions correspond to the first moment in $\zeta$ of the DPD:

$\displaystyle
A_{a_{1}a_{2}}(py=0,y^{2})=\frac{1}{2\pi}\int_{-1}^{1}\mathrm{d}\zeta\,I_{a_{1}a_{2}}(\zeta,y^{2})\,.$
(16)

In Bali:2021gel  we found similar patterns for the twist-two functions and the
reconstructed DPDs themselves regarding their dependence on the distance $y$
and on the quark polarization.

###  2.3 DPDs in the $SU(6)$ quark model

In the following, we consider a simple $SU(6)$-symmetric quark model and
derive its predictions for matrix elements of two-quark currents. The spin-
flavor part of the $SU(6)$-symmetric proton wave function $\ket{p^{\mskip
1.5mu\uparrow}}$ is given by:

$\displaystyle\ket{p^{\mskip 1.5mu\uparrow}}$
$\displaystyle=\frac{1}{3\sqrt{2}}\left[\ket{u^{\mskip 1.5mu\uparrow}u^{\mskip
1.5mu\downarrow}d^{\mskip 1.5mu\uparrow}}+\ket{u^{\mskip
1.5mu\downarrow}u^{\mskip 1.5mu\uparrow}d^{\mskip
1.5mu\uparrow}}-2\ket{u^{\mskip 1.5mu\uparrow}u^{\mskip
1.5mu\uparrow}d^{\mskip 1.5mu\downarrow}}+\ket{u^{\mskip
1.5mu\uparrow}d^{\mskip 1.5mu\uparrow}u^{\mskip
1.5mu\downarrow}}+\ket{u^{\mskip 1.5mu\downarrow}d^{\mskip
1.5mu\uparrow}u^{\mskip 1.5mu\uparrow}}-2\ket{u^{\mskip
1.5mu\uparrow}d^{\mskip 1.5mu\downarrow}u^{\mskip 1.5mu\uparrow}}+\right.$
$\displaystyle\quad\left.+\ket{d^{\mskip 1.5mu\uparrow}u^{\mskip
1.5mu\uparrow}u^{\mskip 1.5mu\downarrow}}+\ket{d^{\mskip
1.5mu\uparrow}u^{\mskip 1.5mu\downarrow}u^{\mskip
1.5mu\uparrow}}-2\ket{d^{\mskip 1.5mu\downarrow}u^{\mskip
1.5mu\uparrow}u^{\mskip 1.5mu\uparrow}}\right]\,,$ (17)

where $\uparrow$ ($\downarrow$) indicates polarization along the positive
(negative) $z$ axis. Moreover, the quark operators are expressed as

$\displaystyle\mathcal{O}_{(ud)}$
$\displaystyle=\frac{1}{2}\left[(\bar{u}^{\mskip
1.5mu\uparrow}\gamma^{+}d^{\mskip 1.5mu\uparrow})+(\bar{u}^{\mskip
1.5mu\downarrow}\gamma^{+}d^{\mskip 1.5mu\downarrow})\right]\,,$
$\displaystyle\mathcal{O}_{\Delta(ud)}$
$\displaystyle=\frac{1}{2}\left[(\bar{u}^{\mskip
1.5mu\uparrow}\gamma^{+}d^{\mskip 1.5mu\uparrow})-(\bar{u}^{\mskip
1.5mu\downarrow}\gamma^{+}d^{\mskip 1.5mu\downarrow})\right]\,,$ (18)

and likewise for the other flavor combinations. Considering matrix elements of
the form $\bra{p^{\mskip
1.5mu\uparrow}}(\bar{q}_{1}\gamma^{+}q_{2})(\bar{q}_{3}\gamma^{+}q_{4})\ket{p^{\mskip
1.5mu\uparrow}}$ for the proton state (17), we obtain:

$\displaystyle\bra{p^{\mskip 1.5mu\uparrow}}(\bar{u}^{\mskip
1.5mu\uparrow}\gamma^{+}u^{\mskip 1.5mu\uparrow})(\bar{d}^{\mskip
1.5mu\uparrow}\gamma^{+}d^{\mskip 1.5mu\uparrow})\ket{p^{\mskip
1.5mu\uparrow}}$ $\displaystyle=a\,,$ $\displaystyle\quad\bra{p^{\mskip
1.5mu\uparrow}}(\bar{u}^{\mskip 1.5mu\uparrow}\gamma^{+}u^{\mskip
1.5mu\uparrow})(\bar{d}^{\mskip 1.5mu\downarrow}\gamma^{+}d^{\mskip
1.5mu\downarrow})\ket{p^{\mskip 1.5mu\uparrow}}$ $\displaystyle=4a\,,$
$\displaystyle\bra{p^{\mskip 1.5mu\uparrow}}(\bar{u}^{\mskip
1.5mu\downarrow}\gamma^{+}u^{\mskip 1.5mu\downarrow})(\bar{d}^{\mskip
1.5mu\uparrow}\gamma^{+}d^{\mskip 1.5mu\uparrow})\ket{p^{\mskip
1.5mu\uparrow}}$ $\displaystyle=a\,,$ $\displaystyle\quad\bra{p^{\mskip
1.5mu\uparrow}}(\bar{u}^{\mskip 1.5mu\downarrow}\gamma^{+}u^{\mskip
1.5mu\downarrow})(\bar{d}^{\mskip 1.5mu\downarrow}\gamma^{+}d^{\mskip
1.5mu\downarrow})\ket{p^{\mskip 1.5mu\uparrow}}$ $\displaystyle=0\,,$
$\displaystyle\bra{p^{\mskip 1.5mu\uparrow}}(\bar{u}^{\mskip
1.5mu\uparrow}\gamma^{+}u^{\mskip 1.5mu\uparrow})(\bar{u}^{\mskip
1.5mu\uparrow}\gamma^{+}u^{\mskip 1.5mu\uparrow})\ket{p^{\mskip
1.5mu\uparrow}}$ $\displaystyle=4a\,,$ $\displaystyle\quad\bra{p^{\mskip
1.5mu\uparrow}}(\bar{u}^{\mskip 1.5mu\uparrow}\gamma^{+}u^{\mskip
1.5mu\uparrow})(\bar{u}^{\mskip 1.5mu\downarrow}\gamma^{+}u^{\mskip
1.5mu\downarrow})\ket{p^{\mskip 1.5mu\uparrow}}$ $\displaystyle=a\,,$
$\displaystyle\bra{p^{\mskip 1.5mu\uparrow}}(\bar{u}^{\mskip
1.5mu\downarrow}\gamma^{+}u^{\mskip 1.5mu\downarrow})(\bar{u}^{\mskip
1.5mu\uparrow}\gamma^{+}u^{\mskip 1.5mu\uparrow})\ket{p^{\mskip
1.5mu\uparrow}}$ $\displaystyle=a\,,$ $\displaystyle\quad\bra{p^{\mskip
1.5mu\uparrow}}(\bar{u}^{\mskip 1.5mu\downarrow}\gamma^{+}u^{\mskip
1.5mu\downarrow})(\bar{u}^{\mskip 1.5mu\downarrow}\gamma^{+}u^{\mskip
1.5mu\downarrow})\ket{p^{\mskip 1.5mu\uparrow}}$ $\displaystyle=0\,,$
$\displaystyle\bra{p^{\mskip 1.5mu\uparrow}}(\bar{d}^{\mskip
1.5mu\uparrow}\gamma^{+}u^{\mskip 1.5mu\uparrow})(\bar{u}^{\mskip
1.5mu\uparrow}\gamma^{+}d^{\mskip 1.5mu\uparrow})\ket{p^{\mskip
1.5mu\uparrow}}$ $\displaystyle=a\,,$ $\displaystyle\quad\bra{p^{\mskip
1.5mu\uparrow}}(\bar{d}^{\mskip 1.5mu\uparrow}\gamma^{+}u^{\mskip
1.5mu\uparrow})(\bar{u}^{\mskip 1.5mu\downarrow}\gamma^{+}d^{\mskip
1.5mu\downarrow})\ket{p^{\mskip 1.5mu\uparrow}}$ $\displaystyle=-2a\,,$
$\displaystyle\bra{p^{\mskip 1.5mu\uparrow}}(\bar{d}^{\mskip
1.5mu\downarrow}\gamma^{+}u^{\mskip 1.5mu\downarrow})(\bar{u}^{\mskip
1.5mu\uparrow}\gamma^{+}d^{\mskip 1.5mu\uparrow})\ket{p^{\mskip
1.5mu\uparrow}}$ $\displaystyle=-2a\,,$ $\displaystyle\quad\bra{p^{\mskip
1.5mu\uparrow}}(\bar{d}^{\mskip 1.5mu\downarrow}\gamma^{+}u^{\mskip
1.5mu\downarrow})(\bar{u}^{\mskip 1.5mu\downarrow}\gamma^{+}d^{\mskip
1.5mu\downarrow})\ket{p^{\mskip 1.5mu\uparrow}}$ $\displaystyle=0\,,$ (19)

where it is understood that the four field operators are taken at different
positions as specified by (2.1) and (2). Here the factor $a$ depends on the
orbital part of the proton wave function, which we leave unspecified. The DPDs
are expressed in terms of proton matrix elements of the operators (18), i.e.
$f_{a_{1}a_{2}}\propto\bra{p^{\mskip
1.5mu\uparrow}}\mathcal{O}_{a_{1}}\mathcal{O}_{a_{2}}\ket{p^{\mskip
1.5mu\uparrow}}$. Using the results given in (19), we find:

$\displaystyle f_{ud}$ $\displaystyle=+6\tilde{a}\,,$ $\displaystyle\quad
f_{uu}$ $\displaystyle=+6\tilde{a}\,,$ $\displaystyle\quad f_{(du)(ud)}$
$\displaystyle=-3\tilde{a}\,,$ $\displaystyle f_{\Delta u\Delta d}$
$\displaystyle=-4\tilde{a}\,,$ $\displaystyle\quad f_{\Delta u\Delta u}$
$\displaystyle=+2\tilde{a}\,,$ $\displaystyle\quad f_{\Delta(du)\Delta(ud)}$
$\displaystyle=+5\tilde{a}\,.$ (20)

with an overall factor $\tilde{a}$ that depends again on the orbital part of
the wave function. Since this orbital part is isotropic, there is no
difference regarding the direction of the quark polarization. Hence, the
values obtained for $f_{\Delta u\Delta d}$ are the same for $f_{\delta u\delta
d}$, and likewise for all other flavor combinations. Predictions that are
independent of the factor $\tilde{a}$ are obtained for _ratios_ of DPDs, such
as:

$\displaystyle\frac{f_{(du)(ud)}}{f_{ud}}=-\frac{1}{2}\,,\qquad\quad\frac{f_{(du)(ud)}}{f_{uu}}=-\frac{1}{2}\,,\qquad\quad\frac{f_{ud}}{f_{uu}}=+1\,,$
(21)

and

$\displaystyle\frac{f_{\Delta(du)\Delta(ud)}}{f_{ud}}=+\frac{5}{6}\,,\qquad\quad\frac{f_{\Delta
u\Delta d}}{f_{ud}}=-\frac{2}{3}\,,\qquad\quad\frac{f_{\Delta u\Delta
u}}{f_{uu}}=+\frac{1}{3}\,.$ (22)

Analogous predictions hold for the ratios of the twist-two functions
$A_{a_{1}a_{2}}$ and can be directly checked against lattice results.

##  3 Lattice calculation

In the following, we give an overview of the calculation of two-current matrix
elements on the lattice and briefly review the techniques that are used. A
detailed explanation is given in Bali:2021gel . This shall now be extended to
flavor-changing operators.
Let us first recall that the two-current matrix element (11) of the nucleon at
$y^{0}=0$ is related to the nucleon four-point function
$C_{4\mathrm{pt}}(\vec{\mskip 0.5muy}\mskip 1.5mu,t,\tau)$ by the following
formula:

$\displaystyle\left.M_{ij}(p,y)\vphantom{\sum}\right|_{y^{0}=0}=C^{ij,\vec{\mskip
0.5mup}\mskip 1.5mu}_{\mathrm{4pt}}(\vec{\mskip 0.5muy}\mskip
1.5mu):=2V\sqrt{m^{2}+\vec{\mskip 0.5mup}\mskip
1.5mu^{2}}\left.\frac{C^{ij,\vec{\mskip 0.5mup}\mskip
1.5mu}_{\mathrm{4pt}}(\vec{\mskip 0.5muy}\mskip 1.5mu,t,\tau)}{C^{\vec{\mskip
0.5mup}\mskip 1.5mu}_{\mathrm{2pt}}(t)}\right|_{0\ll\tau\ll t}\,,$ (23)

where $V$ is the spatial lattice volume and the four-point function
$C^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{\mathrm{4pt}}(\vec{\mskip
0.5muy}\mskip 1.5mu,t,\tau)$ is given by:

$\displaystyle C^{ij,\vec{\mskip 0.5mup}\mskip
1.5mu}_{\mathrm{4pt}}(\vec{\mskip 0.5muy}\mskip 1.5mu,t,\tau)$
$\displaystyle:=a^{6}\sum_{\vec{\mskip 0.5muz}\mskip
1.5mu^{\prime},\vec{\mskip 0.5muz}\mskip 1.5mu}e^{-i\vec{\mskip 0.5mup}\mskip
1.5mu(\vec{\mskip 0.5muz}\mskip 1.5mu^{\prime}-\vec{\mskip 0.5muz}\mskip
1.5mu)}\ \left\langle\operatorname{tr}\left\\{P_{+}\mathcal{P}(\vec{\mskip
0.5muz}\mskip 1.5mu^{\prime},t)\ J_{i}(\vec{\mskip 0.5muy}\mskip 1.5mu,\tau)\
J_{j}(\vec{\mskip 0.5mu0}\mskip 1.5mu,\tau)\
\overline{\mathcal{P}}(\vec{\mskip 0.5muz}\mskip
1.5mu,0)\right\\}\right\rangle\,.$ (24)

Here $P_{+}=\left(\mathbb{1}+\gamma_{4}\right)/2$ projects onto positive
parity, and $\overline{\mathcal{P}}(\vec{\mskip 0.5mux}\mskip 1.5mu,t)$ and
$\mathcal{P}(\vec{\mskip 0.5mux}\mskip 1.5mu,t)$ are the nucleon
interpolators, for which we choose:

$\displaystyle\overline{\mathcal{P}}(\vec{\mskip 0.5mux}\mskip 1.5mu,t)$
$\displaystyle:=\left.\epsilon_{abc}\ \left[\bar{u}_{a}(x)\ C\gamma_{5}\
\bar{d}_{b}^{\,T}(x)\right]\bar{u}_{c}(x)\vphantom{\sum}\right|_{x^{4}=t}\,,$
$\displaystyle\mathcal{P}(\vec{\mskip 0.5mux}\mskip 1.5mu,t)$
$\displaystyle:=\left.\epsilon_{abc}\ u_{a}(x)\left[u_{b}^{T}(x)\ C\gamma_{5}\
d_{c}(x)\right]\vphantom{\sum}\right|_{x^{4}=t}\,,$ (25)

where $C$ is the charge conjugation matrix in Dirac space. The two-point
function $C^{\vec{\mskip 0.5mup}\mskip 1.5mu}_{\mathrm{2pt}}(t)$ appearing in
(23) is defined as:

$\displaystyle C^{\vec{\mskip 0.5mup}\mskip 1.5mu}_{\mathrm{2pt}}(t)$
$\displaystyle:=a^{6}\sum_{\vec{\mskip 0.5muz}\mskip
1.5mu^{\prime},\vec{\mskip 0.5muz}\mskip 1.5mu}e^{-i\vec{\mskip 0.5mup}\mskip
1.5mu(\vec{\mskip 0.5muz}\mskip 1.5mu^{\prime}-\vec{\mskip 0.5muz}\mskip
1.5mu)}\ \left\langle\operatorname{tr}\left\\{P_{+}\mathcal{P}(\vec{\mskip
0.5muz}\mskip 1.5mu^{\prime},t)\ \overline{\mathcal{P}}(\vec{\mskip
0.5muz}\mskip 1.5mu,0)\right\\}\right\rangle\,.$ (26)

#### Wick contractions:

The four-point function (24) decomposes into a definite set of Wick
contractions w.r.t. the fermion fields. There are five types of contractions,
which we call $C_{1}$, $C_{2}$, $S_{1}$, $S_{2}$ and $D$. These contractions
are represented by the graphs in figure 3. Notice that, depending on the quark
flavor of the operators, there are several contributions for each contraction
type. For $C_{1}$-type graphs, we denote this by the flavor indices of the
operator insertions $J_{q_{1}q_{2},i}$ and $J_{q_{3}q_{4},j}$, i.e.
$C_{1,q_{1}q_{2}q_{3}q_{4}}$. For the proton, only the contractions
$C_{1,uudd}$, $C_{1,uuuu}$ and $C_{1,duud}$ contribute (together with
$C_{1,dduu}$ and $C_{1,uddu}$, which are obtained by exchanging the two
currents). In our calculation, we consider only proton matrix elements in the
iso-symmetric limit, where the $u$-quark and the $d$-quark have the same mass.
In that case, there are two independent contributions with $C_{2}$\- or
$S_{1}$-topology, namely $C_{2,u}$, $C_{2,d}$, $S_{1,u}$ and $S_{1,d}$. The
indicated flavor refers to the quark line that connects one of the currents
with the proton source.

Figure 3: Depiction of the five types of Wick contractions that contribute to
a nucleon four-point function. In case of $C_{1}$, $C_{2}$ and $S_{1}$ the
explicit contraction depends on all involved flavors. Provided that the quark
masses are taken to be equal, $C_{2}$ depends only on the flavor of the quark
line that is connected to one of the insertions (red) and the source. For the
disconnected diagrams $S_{1}$, $S_{2}$ and $D$ we also indicate the
disconnected parts $G_{3\mathrm{pt}}$ and $G_{2\mathrm{pt}}$ (blue), as well
as the loops $L_{1}$ and $L_{2}$ (orange).

The explicit set of contractions contributing to the four-point function
depends on the quark flavors of the considered operators. In the following we
list the decomposition of all possible matrix elements of two light-quark
operators in a proton:

$\displaystyle\left.M_{uudd,ij}(p,y)\right|_{y^{0}=0}$
$\displaystyle=C^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{1,uudd}(\vec{\mskip
0.5muy}\mskip 1.5mu)+S^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{1,u}(\vec{\mskip
0.5muy}\mskip 1.5mu)+S^{ji,\vec{\mskip 0.5mup}\mskip 1.5mu}_{1,d}(-\vec{\mskip
0.5muy}\mskip 1.5mu)+D^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}(\vec{\mskip
0.5muy}\mskip 1.5mu)\,,$
$\displaystyle\left.M_{uuuu,ij}(p,y)\right|_{y^{0}=0}$
$\displaystyle=C^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{1,uuuu}(\vec{\mskip
0.5muy}\mskip 1.5mu)+C^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{2,u}(\vec{\mskip
0.5muy}\mskip 1.5mu)+C^{ji,\vec{\mskip 0.5mup}\mskip 1.5mu}_{2,u}(-\vec{\mskip
0.5muy}\mskip 1.5mu)+S^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{1,u}(\vec{\mskip
0.5muy}\mskip 1.5mu)+S^{ji,\vec{\mskip 0.5mup}\mskip 1.5mu}_{1,u}(-\vec{\mskip
0.5muy}\mskip 1.5mu)$ $\displaystyle\quad+S_{2}^{ij,\vec{\mskip 0.5mup}\mskip
1.5mu}(\vec{\mskip 0.5muy}\mskip 1.5mu)+D^{ij,\vec{\mskip 0.5mup}\mskip
1.5mu}(\vec{\mskip 0.5muy}\mskip 1.5mu)\,,$
$\displaystyle\left.M_{dddd,ij}(p,y)\right|_{y^{0}=0}$
$\displaystyle=C^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{2,d}(\vec{\mskip
0.5muy}\mskip 1.5mu)+C^{ji,\vec{\mskip 0.5mup}\mskip 1.5mu}_{2,d}(-\vec{\mskip
0.5muy}\mskip 1.5mu)+S^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{1,d}(\vec{\mskip
0.5muy}\mskip 1.5mu)+S^{ji,\vec{\mskip 0.5mup}\mskip 1.5mu}_{1,d}(-\vec{\mskip
0.5muy}\mskip 1.5mu)$ $\displaystyle\quad+S_{2}^{ij,\vec{\mskip 0.5mup}\mskip
1.5mu}(\vec{\mskip 0.5muy}\mskip 1.5mu)+D^{ij,\vec{\mskip 0.5mup}\mskip
1.5mu}(\vec{\mskip 0.5muy}\mskip 1.5mu)\,,$
$\displaystyle\left.M_{duud,ij}(p,y)\right|_{y^{0}=0}$
$\displaystyle=C^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{1,duud}(\vec{\mskip
0.5muy}\mskip 1.5mu)+C^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{2,d}(\vec{\mskip
0.5muy}\mskip 1.5mu)+C^{ji,\vec{\mskip 0.5mup}\mskip 1.5mu}_{2,u}(-\vec{\mskip
0.5muy}\mskip 1.5mu)+S^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{2}(\vec{\mskip
0.5muy}\mskip 1.5mu)\,,$
$\displaystyle\left.M_{uddu,ij}(p,y)\right|_{y^{0}=0}$
$\displaystyle=C^{ji,\vec{\mskip 0.5mup}\mskip 1.5mu}_{1,duud}(-\vec{\mskip
0.5muy}\mskip 1.5mu)+C^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{2,u}(\vec{\mskip
0.5muy}\mskip 1.5mu)+C^{ji,\vec{\mskip 0.5mup}\mskip 1.5mu}_{2,d}(-\vec{\mskip
0.5muy}\mskip 1.5mu)+S^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{2}(\vec{\mskip
0.5muy}\mskip 1.5mu)\,,$ (27)

where we use

$\displaystyle C^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}_{1,uudd}(\vec{\mskip
0.5muy}\mskip 1.5mu)=2V\sqrt{m^{2}+\vec{\mskip 0.5mup}\mskip
1.5mu^{2}}\left.\frac{C^{ij,\vec{\mskip 0.5mup}\mskip
1.5mu}_{1,uudd}(\vec{\mskip 0.5muy}\mskip
1.5mu,t,\tau)}{C_{2\mathrm{pt}}^{\vec{\mskip 0.5mup}\mskip
1.5mu}(t)}\right|_{0\ll\tau\ll t}$ (28)

and similarly for all other Wick contractions.

#### Renormalization:

The lattice operators are renormalized multiplicatively and converted to the
$\overline{\mathrm{MS}}$-scheme using the factors $Z_{i}$:

$\displaystyle
J_{i}^{\overline{\mathrm{MS}}}(y)=Z_{i}J_{i}^{\mathrm{latt}}(y)\,.$ (29)

For $\beta=3.4$, the corresponding values of $Z_{i}$ are RQCD:2020kuu :

$\displaystyle Z_{V}=0.7128\,,\qquad Z_{A}=0.7525\,,\qquad Z_{T}=0.8335\,,$
(30)

for the choice

$\displaystyle\mu=2\leavevmode\nobreak\ \mathrm{GeV}\,.$ (31)

The renormalization of the two-current matrix elements is given by:

$\displaystyle
M^{\overline{\mathrm{MS}}}_{q_{1}q_{2}q_{3}q_{4},i_{1}i_{2}}=Z_{i_{1}}Z_{i_{2}}M^{\mathrm{latt}}_{q_{1}q_{2}q_{3}q_{4},i_{1}i_{2}}\,.$
(32)

#### Technical details on Wick contractions:

Figure 4: Sketch of all Wick contractions that are considered in this work.
Colors have no meaning regarding the evaluation technique but only indicate to
which source the corresponding propagator belongs. We use two evaluation
methods for the $D$ graph, which are drawn at the bottom center and the bottom
right.

A sketch summarizing the techniques used to evaluate the Wick contractions is
given in figure 4. The proton source is realized by a point source. Both the
proton source and the proton sink are momentum smeared Bali:2016lva . The
corresponding smeared point-to-all propagator at the source $z$ is denoted by
$M_{z}^{\Phi,\vec{\mskip 0.5mup}\mskip 1.5mu}(y)$. The contractions where at
least one of the two currents is directly connected to the proton source or
sink require usage of the sequential source technique. $C_{1}$ and $C_{2}$
additionally involve stochastic wall sources for a given timeslice $t$, for
which we use $Z_{2}\times Z_{2}$ sources. We denote the corresponding
propagated stochastic source (”stochastic propagator”) by $\psi_{t}^{(\ell)}$,
where $\ell$ indicates the stochastic source. In the case of $C_{2}$, the
stochastic source timeslice and the timeslice where the stochastic propagator
is evaluated are identical. Therefore, the propagator is improved by removing
terms in the corresponding hopping parameter expansion that are trivially zero
in the exact case but contribute to the stochastic noise. A similar
improvement is performed for the loop $L_{1}$ appearing in $S_{1}$ and $D$,
which is also evaluated using stochastic propagators. Notice that for the $D$
contraction we have two versions, one using stochastic sources for both loops,
and one where we use point sources only for one of the two loops. More details
on the techniques are given in Bali:2021gel . This reference gives also
explicit expressions for the quantities that are evaluated on the lattice,
except for the contraction $C_{1,duud}$, which contributes only in the context
of flavor interference. We find:

$\displaystyle C_{1,duud}^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}(\vec{\mskip
0.5muy}\mskip 1.5mu,t,\tau)$
$\displaystyle=\frac{a^{3}}{N_{\mathrm{st}}}\left.\sum_{\vec{\mskip
0.5mux}\mskip
1.5mu}\sum_{\ell}^{N_{\mathrm{st}}}\left\langle\left[q_{2,t,j}^{T,\vec{\mskip
0.5mup}\mskip 1.5mu,(\ell)}(x)\ q_{1,t,i}^{\vec{\mskip 0.5mup}\mskip
1.5mu,(\ell)}(x+y)\right]\right\rangle\right|_{x^{4}=\tau,y^{4}=0}\,$ (33)

with

$\displaystyle\left(q_{1,t,i}^{\vec{\mskip 0.5mup}\mskip
1.5mu,(\ell)}\right)^{a}_{\alpha}(y)$
$\displaystyle:=\left(\overline{S}_{321,t,i}^{\vec{\mskip 0.5mup}\mskip
1.5mu,(\ell)}\right)^{a}_{\alpha}(y)-\left(\overline{S}_{123,t,i}^{\vec{\mskip
0.5mup}\mskip 1.5mu,(\ell)}\right)^{a}_{\alpha}(y)\,,$
$\displaystyle\left(q_{2,t,j}^{\vec{\mskip 0.5mup}\mskip
1.5mu,(\ell)}\right)^{a}_{\alpha}(y)$
$\displaystyle:=\left[\psi_{t}^{\dagger,(\ell)}(y)\ \gamma_{5}\Gamma_{j}\
M^{\Phi,\vec{\mskip 0.5mup}\mskip 1.5mu}_{z}(y)\right]^{a}_{\alpha}\,,$ (34)

where $\overline{S}$ is defined in equation (A.19) of Bali:2021gel .

#### Lattice Setup:

We extend our simulation of reference Bali:2021gel  using the same lattice
setup. For completeness we list again the corresponding lattice parameters in
table 1.

id | $\beta$ | $a[\mathrm{fm}]$ | $L^{3}\times T$ | $\kappa_{l}$ | $\kappa_{s}$ | $m_{\pi}[\mathrm{MeV}]$ | $m_{K}[\mathrm{MeV}]$ | $m_{\pi}La$
---|---|---|---|---|---|---|---|---
H102 | $3.4$ | $0.0856$ | $32^{3}\times 96$ | $0.136865$ | $0.136549339$ | $355$ | $441$ | $4.9$
Table 1: Details on the gauge ensemble H102, which is employed for our
simulation. We use 990 configurations.

In addition to the contractions $C_{1,uudd}$, $C_{1,uuuu}$, $C_{2,u}$,
$C_{2,d}$, $S_{1,u}$, $S_{1,d}$, $S_{2}$ and $D$, which have been already
calculated, we compute the contraction $C_{1,duud}$ according to (33) for
proton momentum $\vec{\mskip 0.5mup}\mskip 1.5mu=\vec{\mskip 0.5mu0}\mskip
1.5mu$. Like for the other $C_{1}$ contractions, we choose as source-sink
separation $t=12a$ and evaluate the four-point functions for insertion times
$\tau\in[3a,t-3a]$. The corresponding data is fitted to a constant behavior in
$\tau$.

##  4 Results for invariant functions

###  4.1 Data quality

Before we discuss our results in a physical context, we first investigate the
quality of our data. In particular, we focus on the potential presence of
excited states and of lattice artifacts, such as anisotropy effects. We only
treat the contraction $C_{1,duud}$, since all other contributions have been
already investigated in our previous work Zimmermann:2019quf ; Bali:2021gel .

#### Excited states:

In the following we consider the $\tau$-dependence of the four-point
correlator $C_{4\mathrm{pt}}(t,\tau,\vec{\mskip 0.5muy}\mskip 1.5mu)$ defined
in (24) for fixed $\vec{\mskip 0.5muy}\mskip 1.5mu$. If there is a
contamination of the data by excited states, this leads to a curvature along
$\tau$. It turns out that there is in fact no visible dependence of the data
on $\tau$. As an example we show the data of the $\langle V_{0}V_{0}\rangle$
correlator for the quark separation $\vec{\mskip 0.5muy}\mskip
1.5mu=(-3,4,3)a$ in figure 5. The absence of curvature shows that excited
states are sufficiently suppressed at the scale of our statistical
uncertainties. The final value for the two-current correlator ground state is
obtained by a constant fit in the region $3\leq\tau\leq 9$. It has a
reasonably small error and is also represented in figure 5 by the magenta
band.

Figure 5: Dependence on the insertion time $\tau$ of the contractions
$C_{1,uuuu}$ (green), $C_{1,uudd}$ (red), and $C_{1,duud}$ (magenta) for the
current combination $V_{0}V_{0}$ and fixed quark separation $\vec{\mskip
0.5muy}\mskip 1.5mu=(-3,4,3)$. The orange band represents the constant fit
that returns the final value for the given $\vec{\mskip 0.5muy}\mskip 1.5mu$.

#### Anisotropy:

Since we compute on a lattice with periodic boundary conditions, the effect of
so-called “mirror charges” becomes important if the distance between the two
current approaches half the lattice size Burkardt:1994pw ; Alexandrou:2008ru .
As a consequence, there are anisotropy effects whose size depends on the angle
$\theta(\vec{y}\mskip 1.5mu)$ between $\vec{y}$ and the closest diagonal of
the lattice. Moreover, the lattice propagator itself exhibits anisotropy
effects, which become large at small distances Bali:2018spj ; Cichy:2012is .
This is the case for contractions where both operators are connected directly
via a quark propagator, i.e. the $C_{2}$ and $S_{2}$ contributions. In order
to reduce the violation of Lorentz invariance, only regions with
$|\vec{y}\mskip 1.5mu|=y\geq 4a$ are considered. As discussed in section 4.1
of Bali:2021gel , we limit ourselves to distances $y\leq 16a$ and take into
account only data points close to the lattice diagonal, in order to reduce
discretization artifacts and finite volumes effects. Specifically, only data
points which fulfill

$\displaystyle\cos(\theta(\vec{y}\mskip 1.5mu))>0.9$ (35)

are considered. For the contraction $C_{1,duud}$ the most important source of
anisotropy is given by mirror charges, which leads to the saw tooth pattern
already observed in other contractions with $C_{1}$ topology. This is
illustrated in figure 6, where the data points fulfilling the constraint (35)
are plotted in red.

Figure 6: The anisotropy of the $C_{1,duud}$ contribution to $\langle
V_{0}V_{0}\rangle$, revealed by a saw tooth pattern for larger angles $\theta$
between the $\vec{y}$ vector and the nearest lattice diagonal.

###  4.2 Invariant functions

(a) $A_{qq^{\prime}}$ (b) $A_{q\delta q^{\prime}}$ (c) $A_{\delta q\delta
q^{\prime}}$

Figure 7: Twist-two functions for three different polarization combinations
and the flavor combinations $uu$, $ud$, $dd$ and $(du)(ud)$. Only the graphs
$C_{1}$ and $C_{2}$ are taken into account. The hadron momentum used is always
zero, $\vec{\mskip 0.5mup}\mskip 1.5mu=\vec{\mskip 0.5mu0}\mskip 1.5mu$.

In the following, we consider physical matrix elements given by (27). We take
only into account the connected contributions $C_{1}$ and $C_{2}$ for quark
separations fulfilling $4a\leq y\leq 16a$ and $\cos(\theta(y))>0.9$. The
disconnected contribution $S_{1}$ was found to be small (but noisy) compared
to $C_{1}$ and $C_{2}$ in Bali:2021gel , whereas the errors on $D$ were too
large to make any useful statement about its size. As pointed out in
Bali:2021gel , $S_{2}$ is seen to violate Lorentz invariance for $y<7a$,
whereas it is orders of magnitudes smaller than the other contractions for
larger quark distances. For that reason, we do not take into account this
contraction in our physical results. The twist-two functions $A_{a_{1}a_{2}}$
and $B_{\delta q\delta q^{\prime}}$ are obtained from the data by solving the
overdetermined equation system (2.2) for momentum $\vec{\mskip 0.5mup}\mskip
1.5mu=\vec{\mskip 0.5mu0}\mskip 1.5mu$. The system of equations is solved by
$\chi^{2}$-minimization (see section 3.2 of Bali:2020mij ). The corresponding
results are shown in figure 7 for the functions $A_{qq^{\prime}}$, $A_{q\delta
q^{\prime}}$, and $A_{\delta q\delta q^{\prime}}$.111Our flavor generic
notation $q$ and $q^{\prime}$ for indices includes the flavor interference
cases $(ud)$ and $(du)$. Here we show the results from Bali:2021gel  for the
flavor diagonal combinations $uu$, $ud$, and $dd$, and compare them to our new
results for the flavor interference channel. The latter are significantly
different from zero and tend to be of similar size as the $dd$ contribution.
Moreover, in the case of $A_{qq^{\prime}}$ and $A_{q\delta q^{\prime}}$ we
observe a change in sign for the $(du)(ud)$ combination as a function of $y$.
Hence, we can conclude that flavor interference can indeed be sizable for
DPDs. We note that our results for $B_{\delta q\delta q^{\prime}}$ have rather
large statistical uncertainties for all flavor combinations other than $ud$
and are hence not shown here.

(a) Ratio $A_{(du)(ud)}/A_{ud}$ (b) Ratio $A_{(du)(ud)}/A_{uu}$ (c) Ratio
$A_{uu}/A_{ud}$

Figure 8: Ratios of unpolarized twist-two functions. The blue points
correspond to only considering the $C_{1}$ graphs, and the green ones to
including both $C_{1}$ and $C_{2}$ contributions. The orange dashed line is
the $SU(6)$ model prediction.

(a) Ratio $A_{\delta(du)\delta(ud)}/A_{ud}$ (b) Ratio
$A_{\Delta(du)\Delta(ud)}/A_{ud}$ (c) Ratio $A_{\Delta u\Delta d}/A_{ud}$ (d)
Ratio $A_{\delta u\delta d}/A_{ud}$

(e) Ratio $A_{\Delta u\Delta u}/A_{uu}$ (f) Ratio $A_{\delta u\delta
u}/A_{uu}$

Figure 9: As figure 8, but for ratios of polarized and unpolarized functions.
Note that there is no $C_{2}$ contribution to the flavor combination $ud$ in
panels (c) and (d)

###  4.3 Comparison with $SU(6)$ predictions

In the following, we compare the ratios (21) derived for an $SU(6)$-symmetric
three-quark wave function with the corresponding data we obtain from the
lattice. Note that the evaluation of the DPDs from a three-quark wave function
corresponds to the contraction $C_{1}$. whereas the interpretation of $C_{2}$
involves wave functions with additional quark-antiquark pairs, as discussed in
section 4.1 of Bali:2020mij . We therefore compare the $SU(6)$ predictions
with either $C_{1}$ alone, or with the sum of $C_{1}$ and $C_{2}$. This is
shown in figure 8 for the unpolarized flavor ratios (21) and in figure 9 for
the polarization ratios (22).

In the case of two unpolarized quarks, the $SU(6)$ model results are fairly
consistent with the corresponding $C_{1}$ results. For large quark distances
$y$, this also holds for the sum of $C_{1}$ and $C_{2}$. However, for smaller
$y$, discrepancies become rather large. By contrast, lattice results for
polarization ratios, which are shown in figure 9, disagree strongly with the
$SU(6)$ model.

##  5 Factorization tests

A common strategy for modeling DPDs is to express them in terms of single-
parton distributions, assuming that correlations between the two partons can
be neglected. A corresponding ansatz can also be formulated at the level of
twist-two functions Bali:2020mij ; Bali:2021gel . We now extend this
formulation to the flavor interference case, which we will find to be special
in this context. We limit ourselves to unpolarized quarks in the following

Technically, the factorization ansatz is obtained by inserting a complete set
of eigenstates between the two operators in the DPD matrix element in (2.1)
and neglecting all intermediate states except for the ground state:

$\displaystyle\sideset{}{{}^{\prime}}{\sum}_{\lambda}\bra{p,\lambda}\mathcal{O}_{a_{1}}(y,z_{1})\
\mathcal{O}_{a_{2}}(0,z_{2})\ket{p,\lambda}$
$\displaystyle\quad\approx\sideset{}{{}^{\prime}}{\sum}_{\lambda,\lambda^{\prime}}\int\frac{\mathrm{d}p^{\prime+}\,\mathrm{d}^{2}\boldsymbol{p}^{\prime}}{2p^{\prime+}(2\pi)^{3}}\,e^{-iy(p^{\prime}-p)}\bra{p,\lambda}\mathcal{O}_{a_{1}}(0,z_{1})\ket{p^{\prime},\lambda^{\prime}}\bra{p^{\prime},\lambda^{\prime}}\mathcal{O}_{a_{2}}(0,z_{2})\ket{p,\lambda}\,.$
(36)

In the flavor diagonal case, the matrix elements on the r.h.s. can be directly
identified with GPD matrix elements
$f_{a}^{\lambda\lambda^{\prime}}(\bar{x},\xi,\boldsymbol{p}^{\prime},\boldsymbol{p})$,
which leads to the factorization formula for DPDs Diehl:2011yj :

$\displaystyle
F_{a_{1}a_{2}}(x_{1},x_{2},\zeta,\boldsymbol{y})\stackrel{{\scriptstyle?}}{{=}}\frac{1}{2(1-\zeta)}\int\frac{\mathrm{d}^{2}\boldsymbol{r}}{(2\pi)^{2}}\,e^{-i\boldsymbol{r}\boldsymbol{y}}\sum_{\lambda\lambda^{\prime}}f_{a_{1}}^{\lambda\lambda^{\prime}}(\bar{x}_{1},-\xi,\boldsymbol{0},-\boldsymbol{r})\
f_{a_{2}}^{\lambda^{\prime}\lambda}(\bar{x}_{2},\xi,-\boldsymbol{r},\boldsymbol{0})\,,$
(37)

where $f$ is defined in equation (5.5) of Bali:2021gel . The notation
$\stackrel{{\scriptstyle?}}{{=}}$ indicates that (37) is an assumption, which
we investigate in the following. $\bar{x}_{i}$ and $\xi$ are functions of
$\zeta$ and/or $x_{i}$:

$\displaystyle\bar{x}_{i}(x_{i},\zeta):=\frac{2x_{i}}{2-\zeta}\,,\qquad\xi(\zeta):=\frac{\zeta}{2-\zeta}\,.$
(38)

(a) $F_{ud}$ (b) $F_{(ud)(du)}$

Figure 10: Depiction of the factorized expression (37) for the flavor-diagonal
case (a), as well as the interference case (b). Figure 11: Support regions for
different values of $\zeta$, blue: physical support region of the DPD, red
dashed: support region of the factorized expression

For flavor-changing operators, the matrix elements in the second line of (36)
are related to GPDs for the transition between a proton and a neutron. In
figure 10 we show a pictorial representation of the factorized expression in
the flavor diagonal and non-diagonal cases. Using isospin symmetry, one can
relate these transition GPDs to ordinary ones Mankiewicz:1997aa :

$\displaystyle\bra{p}\mathcal{O}_{(ud)}\ket{n}=\bra{p}\mathcal{O}_{u}\ket{p}-\bra{p}\mathcal{O}_{d}\ket{p}=\bra{n}\mathcal{O}_{(du)}\ket{p}\,,$
(39)

which leads to the following decomposition of the flavor interference matrix
element in terms of flavor conserving matrix elements:

$\displaystyle\bra{p}\mathcal{O}_{(ud)}\ket{n}\bra{n}\mathcal{O}_{(du)}\ket{p}\
=\Bigl{[}\bra{p}\mathcal{O}_{u}\ket{p}-\bra{p}\mathcal{O}_{d}\ket{p}\Bigr{]}^{2}$
$\displaystyle\quad=\bra{p}\mathcal{O}_{u}\ket{p}\bra{p}\mathcal{O}_{u}\ket{p}-2\bra{p}\mathcal{O}_{u}\ket{p}\bra{p}\mathcal{O}_{d}\ket{p}+\bra{p}\mathcal{O}_{d}\ket{p}\bra{p}\mathcal{O}_{d}\ket{p}\,.$
(40)

Inserting this into the factorization hypothesis (36), we find that we can
relate the flavor interference DPDs to a combination of flavor diagonal ones:

$\displaystyle
F_{(ud)(du)}(x_{1},x_{2},\zeta,\boldsymbol{y})=F_{uu}(x_{1},x_{2},\zeta,\boldsymbol{y})-2F_{ud}(x_{1},x_{2},\zeta,\boldsymbol{y})+F_{dd}(x_{1},x_{2},\zeta,\boldsymbol{y})\,.$
(41)

Moreover, taking Mellin moments and performing a Fourier transform w.r.t.
$\zeta$, we find for the invariant functions:

$\displaystyle
A_{(ud)(du)}(py,y^{2})=A_{uu}(py,y^{2})-2A_{ud}(py,y^{2})+A_{dd}(py,y^{2})\,.$
(42)

Notice that the ordering of the operators is important. If we considered
$\mathcal{O}_{(du)}\,\mathcal{O}_{(ud)}$ instead of
$\mathcal{O}_{(ud)}\,\mathcal{O}_{(du)}$, the intermediate state with lowest
energy would not be a nucleon and the corresponding transition GPDs could not
be related with those in the proton. Hence, we are restricted to the operator
ordering given in (40). In Bali:2021gel , we pointed out that the regions of
support of the two sides of (37) differ to a degree that depends on the value
of $\zeta$. In particular, the mismatch is worst for $\zeta<0$ (see figure
11). Hence, we used (37) only for $\zeta\geq 0$ whilst for $\zeta<0$ we
derived a factorized expression for the order
$\mathcal{O}_{a_{2}}\mathcal{O}_{a_{1}}$. As we just explained, we cannot do
that in the flavor interference case, where we have to keep one ordering of
the operators and integrate over the entire $\zeta$-region. This leads to:

$\displaystyle
A_{(ud)(du)}(py=0,y^{2})\stackrel{{\scriptstyle?}}{{=}}\frac{1}{4\pi^{2}}\int_{-1}^{1}\mathrm{d}\zeta\;\frac{(1-\frac{\zeta}{2})^{2}}{2(1-\zeta)}\int\mathrm{d}r\
r\,J_{0}(yr)\sum_{\lambda\lambda^{\prime}}\int\mathrm{d}x_{1}\
\int\mathrm{d}x_{2}\ $
$\displaystyle\qquad\times\left[f_{u}^{\lambda\lambda^{\prime}}(x_{1},-\xi,\boldsymbol{0},-\boldsymbol{r})f_{u}^{\lambda^{\prime}\lambda}(x_{2},\xi,-\boldsymbol{r},\boldsymbol{0})-2f_{u}^{\lambda\lambda^{\prime}}(x_{1},-\xi,\boldsymbol{0},-\boldsymbol{r})f_{d}^{\lambda^{\prime}\lambda}(x_{2},\xi,-\boldsymbol{r},\boldsymbol{0})\right.$
$\displaystyle\qquad\quad\left.+f_{d}^{\lambda\lambda^{\prime}}(x_{1},-\xi,\boldsymbol{0},-\boldsymbol{r})f_{d}^{\lambda^{\prime}\lambda}(x_{2},\xi,-\boldsymbol{r},\boldsymbol{0})\right]\,,$
(43)

where $y=|\boldsymbol{y}|$ and $r=|\boldsymbol{r}|$. Notice that the
expressions in the square brackets are rotationally invariant w.r.t. to the
momentum $\boldsymbol{r}$, which allowed us to perform the angular part of the
integration over $\boldsymbol{r}$. After integration over $x_{1}$ and $x_{2}$,
they can be expressed in terms of Pauli and Dirac form factors, $F_{1}$ and
$F_{2}$. For the second term in the square brackets of (43), we obtain:

$\displaystyle\frac{1}{2}\sum_{\lambda\lambda^{\prime}}\int\mathrm{d}x_{1}\int\mathrm{d}x_{2}\
f_{u}^{\lambda\lambda^{\prime}}(x_{1},-\xi,\boldsymbol{0},-\boldsymbol{r})\
f_{d}^{\lambda^{\prime}\lambda}(x_{2},\xi,-\boldsymbol{r},\boldsymbol{0})$
$\displaystyle\qquad=K_{1}(\zeta)F_{1}^{u}(t)F_{1}^{d}(t)-K_{2}(\zeta)\left[F_{1}^{u}(t)F_{2}^{d}(t)+F_{1}^{d}(t)F_{2}^{u}(t)\right]$
$\displaystyle\qquad\quad+\left(K_{3}(\zeta)+\frac{\boldsymbol{r}^{2}}{4m^{2}}K_{4}(\zeta)\right)F_{2}^{u}(t)F_{2}^{d}(t)\,,$
(44)

with

$\displaystyle t(\zeta,\boldsymbol{r}^{2})$
$\displaystyle:=-\frac{\zeta^{2}m^{2}+\boldsymbol{r}^{2}}{1-\zeta}\,,\qquad
K_{1}(\zeta):=1-K_{2}(\zeta)\,,\qquad
K_{2}(\zeta):=\frac{\zeta^{2}}{(2-\zeta)^{2}}\,,$ $\displaystyle K_{3}(\zeta)$
$\displaystyle:=\frac{\bigl{(}K_{2}(\zeta)\bigr{)}^{2}}{K_{1}(\zeta)}\,,\qquad
K_{4}(\zeta):=\frac{1}{1-\zeta}\,.$ (45)

Analogous expressions hold for the other terms.

For the Dirac and Pauli form factors, we take the results of lattice
simulations obtained with the same ensemble used in the present work. We
fitted these data to several parametrizations in Bali:2021gel . In the
following plots, bands correspond to the range obtained when inserting these
different parametrizations into the factorization formula.

(a) $C_{1}+C_{2}$ for flavors $uu$ (b) $C_{1}+C_{2}$ for flavors $ud$

Figure 12: Factorization of the unpolarized flavor conserving twist-two
functions $A_{ud}$ and $A_{uu}$. The blue band shows the result for the
factorization for the reduced support region that has been published in
Bali:2021gel . The orange band represents the factorization result where (37)
has been directly applied without exchanging the order of operators.

In order to check the impact of integrating (36) over the entire
$\zeta$-region instead of exchanging the order of operators, we repeat the
factorization analysis of the flavor diagonal DPDs given in Bali:2021gel  by
using the correspondingly modified ansatz. The result of this is shown in
figure 12, along with our original result from Bali:2021gel , where (36) was
only used for $\zeta\geq 0$. We see that the difference due to the increased
unphysical support region is considerable when going to quark distances below
about $0.8\leavevmode\nobreak\ \mathrm{fm}$. This shows that integrating (36)
over the entire $\zeta$-region yields a very poor approximation for small
quark distances. Unfortunately, for the flavor interference contribution this
is the only possibility we have found.

Figure 13: Factorization of the unpolarized flavor interference twist-two
function $A_{(ud)(du)}$. The orange band represents the result obtained from
(43). This is compared to the corresponding results obtained from our four-
point function calculation. Here we show the data points for $A_{(ud)(du)}$
(green) and the r.h.s. of the flavor relation (42) (red).

The result for flavor interference DPDs is plotted in figure 13 and compared
to the corresponding result of the direct calculation (green). Based on the
comparison in figure 13, we consider the factorization formula that uses the
full $\zeta$-region as unreliable for $y<0.8\leavevmode\nobreak\ \mathrm{fm}$.
In the region $0.8\leavevmode\nobreak\ \mathrm{fm}<y<1.1\leavevmode\nobreak\
\mathrm{fm}$, we observe that the factorization ansatz yields the wrong sign.
The only approximate agreement with the four-point result can be observed for
very large quark distances, where the result is close to zero. Hence, we
conclude that within a wide range of $y$, the factorization ansatz is either
not reliable because of the mismatch of the support regions, or it fails.

Let us finally discuss the flavor relation (42), which follows from our
factorization formula but is more general. We see in figure 13 that the left-
and right-hand sides of this relation typically differ by a factor around $2$
in the full $y$ range of the plot. The relation is thus satisfied somewhat
better than the factorization hypothesis for $A_{(ud)(du)}$ in terms of three-
point correlation functions.

##  6 Conclusions

We extended our work on nucleon DPDs Bali:2021gel  by considering flavor
interference distributions. We calculated the corresponding Wick contractions
and extracted the twist-two functions for the flavor combination $(ud)(du)$
for all relevant quark polarizations. This we did for zero proton momentum,
which implies $py=0$. We find that the resulting signal for flavor
interference is non-negligible and of the same order of magnitude as the $dd$
contributions. This suggests that at large parton momentum fractions, flavor
interference DPDs may not be negligible.

Several relations between the different flavor and polarization combinations
can be derived in a simple $SU(6)$ quark model, which can then be compared to
the corresponding lattice quantities. We observe good agreement between the
lattice data and the $SU(6)$ results for unpolarized quarks if we only take
into account the contribution of the $C_{1}$ contraction, which is most
consistent with the picture of a three-quark wave function. If one considers
the full connected contribution ($C_{1}+C_{2}$) there is a significant
mismatch, which increases for smaller quark distances. The $SU(6)$ model
completely fails in the case of two-parton matrix elements for polarized
quarks. Finally, we considered the factorization of flavor interference DPDs
in terms of GPDs. We found that, depending on the considered quark distance,
the factorization ansatz is either not reliable or it does not work in the
region where the signal is clearly non-zero.

## Acknowledgments

This work was in part supported by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) by project SCHA 458/23 and FOR 2926, grant number
40965613. The project leading to this publication received funding from the
Excellence Initiative of Aix-Marseille University - A*MIDEX, a French
“Investissements d’Avenir” programme, AMX-18-ACE-005. We gratefully
acknowledge helpful discussions with Gunnar Bali.

## References

  * (1) P.V. Landshoff and J.C. Polkinghorne, _Calorimeter Triggers for Hard Collisions_ , _Phys. Rev. D_ 18 (1978) 3344.
  * (2) R. Kirschner, _Generalized Lipatov-Altarelli-Parisi Equations and Jet Calculus Rules_ , _Phys. Lett. B_ 84 (1979) 266.
  * (3) H.D. Politzer, _Power Corrections at Short Distances_ , _Nucl. Phys. B_ 172 (1980) 349.
  * (4) N. Paver and D. Treleani, _Multi - Quark Scattering and Large $p_{T}$ Jet Production in Hadronic Collisions_, _Nuovo Cim. A_ 70 (1982) 215.
  * (5) V.P. Shelest, A.M. Snigirev and G.M. Zinovev, _The Multiparton Distribution Equations in QCD_ , _Phys. Lett. B_ 113 (1982) 325.
  * (6) M. Mekhfi, _Multiparton Processes: An Application to Double Drell-Yan_ , _Phys. Rev. D_ 32 (1985) 2371.
  * (7) T. Sjostrand and M. van Zijl, _Multiple Parton-parton Interactions in an Impact Parameter Picture_ , _Phys. Lett. B_ 188 (1987) 149.
  * (8) B. Blok, Y. Dokshitzer, L. Frankfurt and M. Strikman, _The Four jet production at LHC and Tevatron in QCD_ , _Phys. Rev. D_ 83 (2011) 071501 [1009.2714].
  * (9) J.R. Gaunt and W.J. Stirling, _Double Parton Scattering Singularity in One-Loop Integrals_ , _JHEP_ 06 (2011) 048 [1103.1888].
  * (10) M.G. Ryskin and A.M. Snigirev, _A Fresh look at double parton scattering_ , _Phys. Rev. D_ 83 (2011) 114047 [1103.3495].
  * (11) B. Blok, Y. Dokshitser, L. Frankfurt and M. Strikman, _pQCD physics of multiparton interactions_ , _Eur. Phys. J. C_ 72 (2012) 1963 [1106.5533].
  * (12) Diehl, Markus and Ostermeier, Daniel and Schäfer, Andreas, _Elements of a theory for multiparton interactions in QCD_ , _JHEP_ 03 (2012) 089 [1111.0910].
  * (13) A.V. Manohar and W.J. Waalewijn, _A QCD Analysis of Double Parton Scattering: Color Correlations, Interference Effects and Evolution_ , _Phys. Rev. D_ 85 (2012) 114009 [1202.3794].
  * (14) A.V. Manohar and W.J. Waalewijn, _What is Double Parton Scattering?_ , _Phys. Lett. B_ 713 (2012) 196 [1202.5034].
  * (15) M.G. Ryskin and A.M. Snigirev, _Double parton scattering in double logarithm approximation of perturbative QCD_ , _Phys. Rev. D_ 86 (2012) 014018 [1203.2330].
  * (16) J.R. Gaunt, _Single Perturbative Splitting Diagrams in Double Parton Scattering_ , _JHEP_ 01 (2013) 042 [1207.0480].
  * (17) B. Blok, Y. Dokshitzer, L. Frankfurt and M. Strikman, _Perturbative QCD correlations in multi-parton collisions_ , _Eur. Phys. J. C_ 74 (2014) 2926 [1306.3763].
  * (18) M. Diehl, J.R. Gaunt and K. Schönwald, _Double hard scattering without double counting_ , _JHEP_ 06 (2017) 083 [1702.06486].
  * (19) B. Cabouat, J.R. Gaunt and K. Ostrolenk, _A Monte-Carlo Simulation of Double Parton Scattering_ , _JHEP_ 11 (2019) 061 [1906.04669].
  * (20) B. Cabouat and J.R. Gaunt, _Combining single and double parton scatterings in a parton shower_ , _JHEP_ 10 (2020) 012 [2008.01442].
  * (21) Axial Field Spectrometer collaboration, _Double Parton Scattering in $pp$ Collisions at $\sqrt{s}=63$-GeV_, _Z. Phys. C_ 34 (1987) 163.
  * (22) CDF collaboration, _Study of four jet events and evidence for double parton interactions in $p\bar{p}$ collisions at $\sqrt{s}=1.8$ TeV_, _Phys. Rev. D_ 47 (1993) 4857.
  * (23) CDF collaboration, _Measurement of double parton scattering in $\bar{p}p$ collisions at $\sqrt{s}=1.8$ TeV_, _Phys. Rev. Lett._ 79 (1997) 584.
  * (24) D0 collaboration, _Double parton interactions in $\gamma$+3 jet events in $pp^{-}$ bar collisions $\sqrt{s}=1.96$ TeV_, _Phys. Rev. D_ 81 (2010) 052012 [0912.5104].
  * (25) D0 collaboration, _Evidence for simultaneous production of $J/\psi$ and $\Upsilon$ mesons_, _Phys. Rev. Lett._ 116 (2016) 082002 [1511.02428].
  * (26) LHCb collaboration, _Observation of double charm production involving open charm in pp collisions at $\sqrt{s}$ = 7 TeV_, _JHEP_ 06 (2012) 141 [1205.0975].
  * (27) ATLAS collaboration, _Measurement of hard double-parton interactions in $W(\to l\nu)$\+ 2 jet events at $\sqrt{s}$=7 TeV with the ATLAS detector_, _New J. Phys._ 15 (2013) 033038 [1301.6872].
  * (28) CMS collaboration, _Study of Double Parton Scattering Using W \+ 2-Jet Events in Proton-Proton Collisions at $\sqrt{s}$ = 7 TeV_, _JHEP_ 03 (2014) 032 [1312.5729].
  * (29) ATLAS collaboration, _Study of hard double-parton scattering in four-jet events in pp collisions at $\sqrt{s}=7$ TeV with the ATLAS experiment_, _JHEP_ 11 (2016) 110 [1608.01857].
  * (30) LHCb collaboration, _Production of associated Y and open charm hadrons in pp collisions at $\sqrt{s}=7$ and 8 TeV via double parton scattering_, _JHEP_ 07 (2016) 052 [1510.05949].
  * (31) LHCb collaboration, _Observation of Enhanced Double Parton Scattering in Proton-Lead Collisions at $\sqrt{s_{NN}}$ =8.16 TeV_, _Phys. Rev. Lett._ 125 (2020) 212001 [2007.06945].
  * (32) CMS collaboration, _Observation of same-sign WW production from double parton scattering in proton-proton collisions at $\sqrt{s}$ = 13 TeV_, 2206.02681.
  * (33) G.S. Bali, L. Castagnini, M. Diehl, J.R. Gaunt, B. Gläßle, A. Schäfer et al., _Double parton distributions in the pion from lattice QCD_ , _JHEP_ 02 (2021) 067 [2006.14826].
  * (34) G.S. Bali, M. Diehl, B. Gläßle, A. Schäfer and C. Zimmermann, _Double parton distributions in the nucleon from lattice QCD_ , _JHEP_ 09 (2021) 106 [2106.03451].
  * (35) H.-W. Lin et al., _Parton distributions and lattice QCD calculations: a community white paper_ , _Prog. Part. Nucl. Phys._ 100 (2018) 107 [1711.07916].
  * (36) J.-H. Zhang, _Double Parton Distributions from Euclidean Lattice_ , 2304.12481.
  * (37) M. Jaarsma, R. Rahn and W.J. Waalewijn, _Towards double parton distributions from first principles using Large Momentum Effective Theory_ , _JHEP_ 12 (2023) 014 [2305.09716].
  * (38) X. Ji, _Parton Physics on a Euclidean Lattice_ , _Phys. Rev. Lett._ 110 (2013) 262002 [1305.1539].
  * (39) X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang and Y. Zhao, _Large-momentum effective theory_ , _Rev. Mod. Phys._ 93 (2021) 035005 [2004.03543].
  * (40) H.-M. Chang, A.V. Manohar and W.J. Waalewijn, _Double Parton Correlations in the Bag Model_ , _Phys. Rev. D_ 87 (2013) 034009 [1211.3132].
  * (41) M. Rinaldi, S. Scopetta and V. Vento, _Double parton correlations in constituent quark models_ , _Phys. Rev. D_ 87 (2013) 114021 [1302.6462].
  * (42) W. Broniowski and E. Ruiz Arriola, _Valence double parton distributions of the nucleon in a simple model_ , _Few Body Syst._ 55 (2014) 381 [1310.8419].
  * (43) M. Rinaldi, S. Scopetta, M. Traini and V. Vento, _Double parton correlations and constituent quark models: a Light Front approach to the valence sector_ , _JHEP_ 12 (2014) 028 [1409.1500].
  * (44) W. Broniowski, E. Ruiz Arriola and K. Golec-Biernat, _Generalized Valon Model for Double Parton Distributions_ , _Few Body Syst._ 57 (2016) 405 [1602.00254].
  * (45) T. Kasemets and A. Mukherjee, _Quark-gluon double parton distributions in the light-front dressed quark model_ , _Phys. Rev. D_ 94 (2016) 074029 [1606.05686].
  * (46) M. Rinaldi, S. Scopetta, M.C. Traini and V. Vento, _Correlations in Double Parton Distributions: Perturbative and Non-Perturbative effects_ , _JHEP_ 10 (2016) 063 [1608.02521].
  * (47) M. Rinaldi and F.A. Ceccopieri, _Relativistic effects in model calculations of double parton distribution function_ , _Phys. Rev. D_ 95 (2017) 034040 [1611.04793].
  * (48) M. Rinaldi, S. Scopetta, M. Traini and V. Vento, _A model calculation of double parton distribution functions of the pion_ , _Eur. Phys. J. C_ 78 (2018) 781 [1806.10112].
  * (49) A. Courtoy, S. Noguera and S. Scopetta, _Double parton distributions in the pion in the Nambu–Jona-Lasinio model_ , _JHEP_ 12 (2019) 045 [1909.09530].
  * (50) W. Broniowski and E. Ruiz Arriola, _Double parton distribution of valence quarks in the pion in chiral quark models_ , _Phys. Rev. D_ 101 (2020) 014019 [1910.03707].
  * (51) W. Broniowski and E. Ruiz Arriola, _Double parton distributions of the pion in the NJL model_ , _PoS_ LC2019 (2020) 031 [2001.00883].
  * (52) RQCD collaboration, _Nonperturbative Renormalization in Lattice QCD with three Flavors of Clover Fermions: Using Periodic and Open Boundary Conditions_ , _Phys. Rev. D_ 103 (2021) 094511 [2012.06284].
  * (53) G.S. Bali, B. Lang, B.U. Musch and A. Schäfer, _Novel quark smearing for hadrons with high momenta in lattice QCD_ , _Phys. Rev. D_ 93 (2016) 094515 [1602.05525].
  * (54) RQCD collaboration, _Two-current correlations and DPDs for the nucleon on the lattice_ , _PoS_ LATTICE2019 (2019) 040 [1911.05051].
  * (55) M. Burkardt, J.M. Grandy and J.W. Negele, _Calculation and interpretation of hadron correlation functions in lattice QCD_ , _Annals Phys._ 238 (1995) 441 [hep-lat/9406009].
  * (56) C. Alexandrou and G. Koutsou, _A Study of Hadron Deformation in Lattice QCD_ , _Phys. Rev. D_ 78 (2008) 094506 [0809.2056].
  * (57) G.S. Bali, V.M. Braun, B. Gläßle, M. Göckeler, M. Gruber, F. Hutzler et al., _Pion distribution amplitude from Euclidean correlation functions: Exploring universality and higher-twist effects_ , _Phys. Rev. D_ 98 (2018) 094507 [1807.06671].
  * (58) K. Cichy, K. Jansen and P. Korcyl, _Non-perturbative renormalization in coordinate space for $N_{f}=2$ maximally twisted mass fermions with tree-level Symanzik improved gauge action_, _Nucl. Phys. B_ 865 (2012) 268 [1207.0628].
  * (59) L. Mankiewicz, G. Piller and T. Weigl, _Hard leptoproduction of charged vector mesons_ , _Phys. Rev. D_ 59 (1999) 017501 [hep-ph/9712508].