Double parton distribution of valence quarks in the pion in chiral quark models

The valence double parton distribution of the pion is analyzed in the framework of chiral quark models, where in the chiral limit factorization between the longitudinal and transverse degrees of freedom occurs. This feature leads, at the quark-model scale, to a particularly simple distribution of the form D(x1, x2,y) = δ(1−x1−x2)F (y), where x1,2 are the longitudinal momentum fractions carried the valence quark and antiquark, and y is the separation of their transverse positions. For y = 0 this result complies immediately to the Gaunt-Sterling sum rules. The DGLAP evolution to higher scales is carried out in terms of the Mellin moments. We then explore its role on the longitudinal correlation quantified with the ratio of the double distribution to the product of single distributions, D(x1, x2,y)/D(x1)D(x2). We point out that the ratios of moments 〈x1x2 〉/〈x1 〉〈x2 〉 are independent of the evolution, providing particularly suitable measures to be tested in the upcoming lattice simulations. The transverse form factor F (y) is obtained in variants of the Nambu–Jona-Lasinio model with the spectral and Pauli-Villars regularizations. Interestingly, with the spectral regularization of the model, the effective cross section for the double parton scattering of pions is exactly equal to the geometric cross section, σeff = π〈y〉 and yields about 20 mb.


I. INTRODUCTION
The pion, which according to the constituent quark model is composed of a quark q and antiquarkq pair, enjoys being the would-be Goldstone boson of the spontaneously broken chiral symmetry -a feature which explains its particularly low mass as compared to other hadrons. There exist many observables which probe indirectly the consequences of the pion being both a pseudo-Goldstone boson and a composite qq state, indicating its extended non-elementary structure at intermediate energies. These include the electromagnetic and gravitational form factors (FF), transition form factors, electromagnetic polarizabilities, as well as features of their the mutual interactions, such as the Gasser-Leutwyler coefficients. Likewise, for high-energy processes where hard-soft factorization holds, the soft matrix elements within the pion state yield the single parton distribution functions (sPDF), the corresponding generalized parton distributions (GPD), transverse momentum distributions (TMD), parton distribution amplitudes (PDA), etc.
A strong motivation to investigate the double parton distribution functions (dPDFs) in the pion comes from the more and more promising prospects of the lattice QCD studies [1,2]. In this paper we address the valence dPDF of the pion in the simplest covariant fieldtheoretic approach to its structure, namely the Nambu-Jona-Lasinio (NJL) model [3,4], where the pion arises as a Goldstone boson of the spontaneously broken chiral symmetry and as a qq relativistic bound state of the Bethe-Salpeter equation. This model is not renormalizable and needs a suitable regularization. Here we use the Pauli-Villars (PV) regularization which preserves chiral symmetry and gauge invariance [5] (for a review see, e.g., [6]). We also explore the Spectral Quark Model (SQM) [7,8], where particularly simple analytic formulas can obtained.
The story of dPDFs is rather old [9], with first experimental traces of the possible double parton scattering (DPS) in multi-jet events in pp collisions reported by the Axial Field Spectrometer Collaboration at the CERN ISR [10] and in pp collisions by the CDF Collaboration at Fermilab [11,12]. With a renewed interest in the LHC era, several enlightening theoretical works (see [13][14][15][16][17] and references therein) preceded the measurement of W production in association with two jets by the AT-LAS Collaboration [18], which provided a strong case for the need of DPS. Thereby dPDFs have become more directly accessible to experimental scrutiny and quantitative analysis (for reviews see, e.g., [19] for implications in production processes). For a recent state of affairs see [20].
As any partonic quantity, dPDFs are scale dependent and their evolution is dictated by a generalized form of the DGLAP equations [21,22], where gluon radiation is explicitly implemented. In the case of the pion, the correlation structure is particularly relevant, as it indicates to what extent the quark and anti-quark are entangled in the presence of additional gluons. While it will be hard to measure dPDF for the pion experimentally, some of its features will become soon available from lattice QCD [1,2], thus offering a unique chance of testing the internal structure of the pion and, more specifically, the correlations between its qq constituents. However, the quest for dPDFs is non-perturbative, whereas the evolution only relates these quantities perturbatively at different scales, saying nothing about their absolute determination at a given scale. Hence the importance of non-perturbative modeling.
Similarly to the previous studies of sPDFs, we assume that there exists a reference scale µ 0 where the pion is just a qq state. This is suggested by the observation that the momentum fraction carried by the quarks at the scale µ = 4 GeV is about 42−44% [23,24] and that it increases according to LO and NLO DGLAP evolution with decreasing µ. Then at about µ 0 ∼ 320 MeV the valence quark and antiquark carry all the momentum. While this is admittedly a low scale, it has been checked that switching from LO to NLO evolution in model calculations does not change more than 10% the pion sPDFs [25], justifying the approach.
Chiral quark models offer a view of the pion as a relativistic qq Goldstone boson, departing from the naive quark model. One of the important features of these models is that they are constructed to satisfy the gauge invariance and the chiral Ward-Takahashi identities. This is important to recover the current algebra results at low energies and to guarantee proper normalization of form factors and PDFs. Moreover, the imposition of relativity via the Bethe-Salpeter equation implements proper support for partonic quantities. Despite these strong constraints, there is still freedom in enforcing proper regularization schemes. In this paper we focus for definiteness on the NJL model with the PV regularization [5,6] for which sPDFs [26], TMDs [27], PDAs [28], GPDs [29] have been computed. We also explore the SQM, which is particularly simple in the resulting formulas and offers a qualitatively good description both at low and at higher energies [8,30,31]. These variants of the NJL model have been used in the past to describe the FF, PDFs, PDAs, GPDs, TMDs, comparing favorably both to experimental data as well as available lattice results.
After some of the results presented here were advertised [32], a relevant preprint by Courtoy, Noguera, and Scoppeta has been released [33], addressing in a comprehensive way the dPDF calculation in the NJL model. Our results confirm closely the basic findings [33], despite different regularization schemes. 1 Here we also fully present the results of the QCD evolution for dPDF and discuss in detail the issue of the longitudinal partonic correlation.
It is worth noting some other non-perturbative studies of dPDFs in quark models for the nucleon, such as in the MIT bag [34] or in the constituent quark model model [35], also with the light-front wave functions [36,37]. We stress that the longitudinal correlation features in these models are qualitatively understood in the framework of the valon model [38,39], as shown in [40,41]. 2 The paper is organized as follows. After introducing the definitions and our notation in Sec. II, we present the NJL model and the calculation of the valence dPDF of the pion in Sec. III. A recurrent issue in this kind of calculations regards the interplay between regularization and positivity, a topic which we address in Section IV. The QCD DGLAP evolution and the corresponding matching condition for the model are discussed in Sec. V, with numerical results presented in Sec. VI. In Sec. VII we present ratios of moments of dPDF to product of moments of sPDFs, which are scale independent for the valence distributions. These ratios could most directly be tested in lattice calculations. The transverse form factor depends somewhat on the variant of the model, as shown in Sec. VIII, where also a few related observables of interest are evaluated and discussed. Finally, in Sec. IX we come to our summary.

II. BASIC DEFINITIONS AND PROPERTIES
In this section we list for completeness (and to establish our notation) the basic definitions. The spin-averaged sPDF and dPDF (see [43] and references therein) of a hadron with momentum p involve forward matrix elements of parton bilinear operators, namely with j 1,2 indicating the parton species, the bold face denoting the transverse vectors, and y playing the role of the transverse distance between the two partons. The light-cone coordinates are a ± = (a 0 ±a 3 )/ √ 2. The longitudinal (i.e., + ) momentum fractions carried by the partons are denoted customarily as x 1 and x 2 .
For the quarks and antiquarks, considered in this paper, the bilocal color-singlet operators are where the summation over color is implicit and the flavor indices are omitted for brevity. As the quark and antiquark coordinates in the individual operators O j in Eq. (1) are not split in the transverse direction, the lightfront gauge A − a = 0 along straight line paths eliminates the Wilson gauge link operators W , setting them to identity (the splitting in y occurs between the color singlet bilinears, which is innocuous). Thus the complications in maintaining the gauge invariance and the choice of the Wilson line, encumbering TMDs [44] and the corresponding QCD evolution [45,46] (for a review see, e.g., [47] and references therein), do not occur in the case of dPDFs.
One may pass to the momentum representation via the Fourier transform where the same symbol D is used, with the argument distinguishing the function from its transform. The double distributions for the case where q = 0, denoted for brevity as acquire a special significance, as they satisfy the Gaunt-Stirling (GS) sum rules [48]. These identities follow straightforwardly [49] when the decomposition of the parton operators in a basis of the light-front wave functions is made, and are essentially statements on completeness of the Fock states as well as flavor and longitudinal momentum conservation. The sum rules read where i val indicates the difference of the parton (i) and antiparton (ī) distribtions, and To satisfy the GS sum rules, we have argued in [40,41] that a practical approach is the top-down method, where one starts with n-body parton distributions with the longitudinal momentum fractions constrained with δ(1 − x 1 − ... − x n ), and subsequently generates distributions with a lower number of partons via marginal projections. This avoids complications of bottom-up attempts in constructing dPDFs from known sPDFs, which cannot be unique and where one encounters problems [48,50]. Also, note a recent study [51] devoted to a verification of the GS sum rules in covariant perturbation theory and in light-cone perturbation theory.
The above discussion presumes a probabilistic interpretation of dPDFs, similarly to sPDFs. As has been known, however (see, e.g., [52]), the need for renormalization of the bare definitions (1) may in principle invalidate positivity, as it inherently involves subtraction. The issue is subtle, as violation of positivity precludes a probabilistic interpretation. Moreover, it should also be recalled that the QCD evolution equations (see Sec. V below) do not preserve positivity; whereas the DGLAP evolution produces positive distributions when evolving upward from Q 0 to Q > Q 0 , for downward evolution, Q < Q 0 , the positivity property may actually be violated [53] (mostly at small x). For an explicit nucleon study and a practical distinction between the (scheme dependent) positivity of parton distributions and the necessary positivity of physically measurable cross sections, we refer the reader to [54]. The generation of negative components with downward evolution is also a feature of the dPDF evolution equations applied here.
Secondly, there is a question of retaining positivity at nonzero y, or a corresponding Fourier conjugate momentum q. The mentioned derivation of the GS sum rules for q = 0 via the Fock-state decomposition involves a summation of moduli squared of the n-parton wave functions, |φ n ({k i }, {x i })| 2 , which is positive and may remain so upon certain renormalization procedures, e.g., the Fock space truncation. With q = 0, the correspond- , which mathematically need not be positive definite, as the wave functions may possess nodes. Of course, in processes where hard-soft factorization holds, thus where the dPDFs are useful, the pertinent cross sections must obviously be positive, which constitutes the true positivity constraints. The issue of positivity is further discussed in Sec. IV.
We note that the factorization of the longitudinal and transverse degrees of freedom, has been a standard working assumption in studies of DPS over the recent years, which so far finds support from dynamical model calculations [34,35]. We remark that in this case the positivity property can be considered separately for D j1j2 (x 1 , x 2 ) and for the form factor F j1j2 (q).

III. THE MODEL
The large-N c evaluation in chiral quark models amounts to evaluating one-quark-loop integrals. For definiteness, we consider the positively charged pion π + , with other states related by the isospin symmetry. The relevant diagram for the valence dPDF in the momentum representation is shown in Fig. 1, where the loop momentum integration is constrained with k + = x 1 p + , as we assign x 1 as the longitudinal momentum fraction carried by the valence u quark. Note that with just two constituents this constraint simultaneously fixes the longitudinal momentum carried by the valence antiquarkd to be x 2 = 1 − x 1 . With the Feynman quark propagator The diagram for the chiral quark model evaluation of the valence dPDF of π + . The loop consists of the constituent quark propagators. The blobs indicate the probing operators γ + , and the pseudoscalar quark-pion coupling iγ5. The integration over dk− d 2 k ⊥ dq− enforces the kinematic constraints of the definition (1).
where M is the constituent quark mass due to the spontaneous breaking of the chiral symmetry, the diagram of Fig. 1 reads where the trace is over color and Dirac indices. We use the pseudoscalar quark-pion coupling M/f × iγ 5 , where f = 86 MeV is the pion decay constant in the chiral limit. An explicit and straightforward evaluation yields [32,33] with the form factor which is ultraviolet log-divergent and needs to be properly regularized, as indicated. For q = 0, upon regularization, the integral reduces to which is equal to 1, as follows from a corresponding expression for the square of the pion decay constant f 2 (cf. for instance [29]). Then the result advocated in [32,33]. As shown in [33], other large-N c diagrams, which formally appear in the effective quark-meson model, do not contribute to D ud .
Thus, at the quark model scale µ 0 , the double parton distribution of the pion takes a product form of the individual longitudinal momentum fractions carried by each valence quark, multiplied by δ(1 − x 1 − x 2 ), coming from the momentum conservation and the fact that at the quark model scale the only constituents in the pion are the two valence partons.
We recall that for the valence sPDF the corresponding result is [26] Taking into account the fact that our model result factorize according to Eq. (10), we discuss separately the longitudinal and transverse structures in the following sections.

IV. REGULARIZATION VS POSITIVITY
In this section we address some technical aspects concerning the interplay between the necessary regularization in a chiral quark model and the positivity which proves the basis for the insightful probabilistic interpretation. As already mentioned, even in QCD this is a subtle issue [52]. Here we will discuss sharp cut-offs either in coordinate space or in momentum spaces. More sophisticated (smooth) schemes are worked out in Section VIII.
On a formal level, it is noteworthy that Eq. (11) has a convolution-like structure that can be rewritten in the form where Ψ(k) = 1 Formally, we may pass to the coordinate space by defining with b the two-dimensional transverse coordinate. From here we get such that the form factor appears as the Fourier transform of a manifestly positive function in b-space. In our case Φ(b) ∼ e −M b /b, hence the integral diverges at small distances. If we put a short distance transverse cut-off, it is not guaranteed that the form factor in momentum space remains positive. 3 In fact, if we impose the normalization condition, F (0) = 1, we get a short distance cutoff b 0 ∼ 0.5 fm and the form factor presents a multinodal structure with alternating signs, the first zero ocurring at |q| ∼ 0.75 MeV. Alternatively, since the formula for F (q) is divergent, a regularization must be imposed. A simple method within momentum space would be to use a purely transverse cut-off Λ ⊥ and fix it according to the normalization condition given by Eq. 12. This is the essence of the light front regularization method applied in Ref. [33], which, requires additional assumptions to be compatible with soft physics and, in particular, to provide a non-vanishing vacuum quark condensate (see e.g. [55,56]). Such a prescription, however, does not preserve positivity [33].
An interesting aspect is that, as noted in [27] (see Appendix A for an analysis in the case of the pion em form factor), usual convolution formulas, as the one discussed above, are only formally correct but not necessarily compatible with gauge invariance or chiral symmetry. 4 The above discussion shows once again that implementation of a finite cut-off regularization in chiral quark models is a nontrivial issue, particularly if a finite regularization is imposed separately on the individual Feynmann diagrams. The effective action method suggested by Eguchi [57] is a symmetry preserving scheme which allows for a proper discussion of chiral symmetry breaking even after the implementation of a regularization method [5,6] and provides a common regularization for all diagramms involving pions, photons, or W and Z bosons. Therefore, her we only consider regularization methods which comply with soft pion physics and chiral symmetry, but then they need not preserve positivity of the dPDF form factor away from the soft limit (particularly above the finite cut-off). Thus, chiral quark models are designed to describe soft physics, whereas large |q|, or small b, evade this requirement, even after regularization, which we encounter below.

V. EVOLUTION AND MATCHING
Parton distribution functions describe nonperturbative properties of a hadron, namely its quark and gluon content, as functions of their longitudinal momentum fraction. These quantities depend on the renormalization scale, which usually is taken to be µ = Q, a typical momentum appearing in the experimental process. Their running with the scale µ follows from the invariance of the physical cross sections. The evolution can only be implemented perturbatively, thus a non-perturbative input is required at a reference scale µ 0 to provide initial conditions for the pertinent evolution equations.
Phenomenological analyses of the pion parameterize the valence, sea, and gluon sPDFs at relatively large scales, µ ∼ 2 GeV [23], or intermediate scales µ ∼ 0.5 MeV [24], with the result that at µ ∼ 2 GeV the momentum fraction carried by the quarks in the pion is 0.42 − 0.44. In fact, in a hadronic model such as NJL applied here, one only has valence quarks and antiquarks, hence it makes sense to fix µ 0 by imposing that these constituents saturate the momentum sum rule, x µ0 = 1. This provides a hadronic scale of µ 0 ∼ 320 MeV (see, e.g., [41]), which we refer to as the quark model scale.
Once the reference scale µ 0 is fixed, the matching condition between the model and QCD observables is imposed by The QCD evolution equations for multi-parton distributions have been derived long ago [21,22]. A simple and numerically efficient method is based on the Mellin moments, similarly to the case of sPDFs (for results of a practical implementation see, e.g., [40] for the nucleon and [41] for the pion). While for ease of notation we do not write explicitly the y dependence, the evolution equations quoted below are still valid for that case [58]. Clearly, for a factorized ansatz, the transverse dependence also factorizes out from the evolution. As we have discussed in the previous section, this is actually the case in chiral quark models in the chiral limit.
One introduces the moments of the sPDFs and dPDFs, and the moments of the QCD splitting functions Then the DGLAP evolution equations are [21,22] (for simplicity, the possible two scales µ 1 and µ 2 for the two operators are set to be equal to a single scale µ) where Partons i, j 1 , and j 2 may in general represent the valence quarks/antiquarks, the sea, or the gluons, whose distributions are coupled. Moreover, the last term in Eq. (22) couples dPDFs to sPDFs, serving as an inhomogeneous term. The evolution acquires a simpler form when one probes just the valence quarks, which is the case considered in this work, where we take π + with j 1 = u and j 2 =d. In that case there are no partons i contributing to the splitting P i→j1j2 , and the inhomogeneous term vanishes. In addition, the nonsiglet valence quarks do not mix with the singlet component (sea quarks and gluons), leaving for the considered π + state the equation Figure 2 represents a typical diagram for the applied evolution in terms of the developing parton cascades along the t-channel.

VI. NUMERICAL RESULTS OF EVOLUTION
The DGLAP evolution described above in an obvious manner preserves the longitudinal-transverse factorization, as all moments are multiplied conventionally with F (q). For that reason we may discuss the evolution of D ud (x 1 , x 2 ).
It is known that the DGLAP evolution preserves the GS sum rules [48]. In our case, the sum rules are trivially satisfied at the quark model scale µ 0 , as with Eqs. (10) and (14) one immediately gets Thus the GS sum rules are satisfied at any scale.
Our numerical results for the evolution of the valence dPDF of the pion (multiplied with x 1 x 2 ) are presented in Fig. 3, where we show it for increasing evolution scale µ. We note the gradual drift towards the asymptotic fixed point δ(x 1 )δ(x 2 ).
The longitudinal correlations can be conveniently quantified with the ratio D ud (x 1 , x 2 )/D u (x 1 )Dd(x 2 ), which would be equal to 1 if no correlations were present. Our results for this measure are shown in Fig. 4. We note that for large scales µ and simultaneously low x 1 and x 2 (which is where dPDF is large) the correlation ratio is within a 20% band around unity, hence there the effect is not very substantial, justifying the product ansatz for dPDF. The largest effect occurs for asymmetric kinematics, with x 1 large and x 2 small, or vice versa, where correlations are strong and positive. The explanation of this behavior is related to the overall conservation of the longitudinal momentum. The qualitative argument here is that at large µ, with many partons present, the constraint of the longitudinal momentum is effective only when one of the considered partons takes a large momentum fraction, thus leaving significantly less of available phase space for the other parton. We remark that a qualitatively similar effect has been found for gluodynamics in [59][60][61].

VII. RATIOS OF MOMENTS
The LO DGLAP evolution of the valence moments, due to the absence of the inhomogeneous term, leads to a simple fact that the ratios x n 1 x m 2 / x n 1 x m 2 for the valence distributions do not depend on the evolution scale, as the evolution ratio factors with the anomalous dimensions cancel out. More explicitly, the evolution for the valence moments reads M nm (µ 0 ), (27) with γ i denoting the anomalous dimensions, thus the cancelation is obvious. Despite its simplicity, this feature is rather remarkable, as it allows for an insight into lower scales, where non-pertubative dynamics sets in, from the information at higher scales.
In particular, in the NJL model we immediately find the following scale independent ratios: The first few ratios are shown in Table VII. The largest ratio is for m = n = 1. Naturally, the moments are the quantities to be probed with the upcoming lattice studies and it is interesting to see if the pattern of Eq. (28) holds. The pattern of the moments is specific to a given (nonperturbative) model, allowing for its scrutiny provided the lattice data would be accurate enough.
The LO DGLAP evolution of the lowest moments of the valence dPDF of the pion is displayed in Fig. 5. We note a fall-off from the quark model scale controlled by the corresponding anomalous dimensions in Eq. (27).

VIII. TRANSVERSE STRUCTURE
The NJL model needs to be regularized to get rid of the ultraviolet divergences, leaving only the soft-momentum degrees of freedom in the dynamics. Unlike the longitudinal structure, the transverse structure depends specifically on the regularization. This is a subtle issue, since as already discussed, the chiral and gauge symmetries must be preserved by the regularization procedure, as has been discussed at length in [5] (for a review see, e.g., [6]). We give here the final recipes as applied to the unregularized result of Eq. (11), both for the PV and SQM implementations. We present first the SQM results, as the final formulas are particularly simple in this case.

A. Spectral Quark Model
In the SQM, one replaces the constituent quark mass M in Eq. (11) with a spectral mass, w, and integrates over w with a spectral weight, over a suitably chosen complex contour C [8]. One may adjust ρ(w) and C in such a way that the construction implements an exact vector-meson dominance principle of the pion electromagnetic form factor, namely F em (Q 2 ) = 1/(1 + Q 2 /m 2 ρ ), with m 2 ρ = 24π 2 f 2 /N c , yielding m ρ = 764 MeV at f = 86 MeV in the chiral limit. In this scheme, the normalization condition is C dwρ(w) = 1.
After computing the spectral integral as well as the k integral, the form factor is given by a very simple formula which is a combination of the dipole (with positive sign) and monopole (with negative sign) form factors. The form factor F (q) is correctly normalized at q = 0, namely F (0) = 1, and vanishes as O(|q| −2 ) for large |q|.
In Fig 6(a) we show the form factor F (q) for SQM (solid line). As we can see, in accordance with Eq. (30), this function becomes negative for |q| > m ρ . Thus the result does nor obey positivity (cf. the discussion at the end of Sec. II). Passing to the configuration space via the Fourier-Bessel transform we get with K n (x) denoting the modified Bessel functions of order n.
The form factor combination bf (b) for SQM is is depicted in Fig 6(b) with a solid line. By definition, this function is normalized to unity, d 2 bf (b) = 1, and the corresponding mean squared radius (msr) is given by This radius is a 2D quantity in relative parton-parton transverse coordinates, and naively one expects a geometric relation to the 3D radius, which would be r 2 = 3 b 2 /2 = (1.07 fm) 2 . It is interesting to compare this measure with another msr regarding the hadron size, for instance the electromagnetic (em) msr, which in SQM is given by r 2 em = 6/m 2 ρ = b 2 dPF /2 = (0.62 fm 2 ). Thus, in this model the transverse size (a 2D object) is about √ 2 larger than the em msr (a 3D object). At large distances f (b) behaves as which falls off exponentially with the mass m ρ , which in a constituent picture is typically m ρ 2M corresponding to a double parton property, opposite to the m ρ /2 M scales in single parton properties [8].
At short distances f (b) becomes negative and divergent, since where γ = 0.5772 . . . is the Euler-Mascheroni constant. Using this asymptotics we get a zero of the function f (b) for the value b 0 ∼ 2e −1−γ /m ρ ∼ 0.1 fm. In Fig. 6(b) we show the b−dependence and, as we can see, the function becomes negative at b ∼ 0.15 fm, which by the way is comparable to the current day spacing of fine QCD lattices, a ∼ 0.1 fm. The effective cross section for DPS is a phenomenologically important quantity [62], defined as For SQM we get which coincides exactly with the geometric cross section (see the discussion in [35]). This number is somewhat larger from similar estimates for the nucleon case, where one obtains about 15 mb (see for instance [63][64][65]) Actually, if we restricted the integration limit to the value where the effective form factor f (b) is positive, we would get about twice the value for the effective cross section, 47 mb.
Within this context there exists a nested bound for σ eff , found recently by Rinaldi and Ceccopieri [66], namely π b 2 ≤ σ eff ≤ 3π b 2 . The derivation of these interesting inequalities rests on the positivity of f (b) (such that b n ≥ b n ). The lower bound is saturated by a monopole (which is positive) as well as our solution (which becomes negative for |q| ≥ m ρ or for b 0.5/m ρ ). To compare, a dipole yields σ eff = 3 2 π b 2 , and a Gaussian profile gives σ eff = 2π b 2 .

B. NJL with Pauli-Villars regularization
The most straightforward regularization scheme satisfying the necessary formal requirements is PV regularization with two subtractions in the coincidence limit [5], where the one-quark-loop quantity A(M 2 ) is replaced with A detailed analysis shows that the rule is only applied to the mass dependence under the integral in Eq. (11), but not to the M 2 factor in front, originating from the coupling constant. The normalization yields the wellknown condition linking the constituent quark mass and the cutoff Λ at a given value of f . Following earlier works, we take M = 300 MeV and f = 86 MeV in the chiral limit, which fixes the PV cut-off to be Λ = 731 MeV. The expression for the valence dPDF form factor is Its asymptotic behavior is similar to SQM. Moreover, the msr is given by 2π 2 f 2 (M 2 + Λ 2 ) 2 = (0.77 fm) 2 (40) and the effective cross section is numerically close to the geometric one, σ eff ∼ π b 2 . The corresponding form factors in the momentum and configuration spaces can also be seen in Fig. 6 (dashed lines). As we can appreciate, they are very similar to the SQM case, including negative values at high momenta ∼ 1 GeV or, equivalently, at short transverse distances b 0.1 fm.

IX. CONCLUSIONS
We have presented a detailed study of the valence double parton distributions in the pion in chiral quark models, with the following basic results: • The model leads in the chiral limit to a factorization of the longitudinal and transverse structure of the valence pion dPDF, in accordance to the findings of [32,33]. At the quark model scale, the longitudinal distribution is of the simple form δ(1 − x 1 − x 2 ), reflecting the momentum conservation with just two constituents.
• The LO DGLAP evolution is carried out in the Mellin space, leading to radiative generation of partons at higher scales, and resulting with distributions at higher scales, accessible in experiments or lattice simulations.
• The Gaunt-Stirling sum rules are explicitly satisfied in our approach.
• The longitudinal correlations, quantified as the departure of dPDF from the product of two sPDFs, are studied in detail. An assessment of the validity of the frequently used product ansatz is made, with the result that at high evolution scales and small momentum fractions of both constituents it works to a good approximation.
• Simple expressions are obtained for the transverse form factor of the valence dPDF in the applied regularizations of the chiral quark models. The issue of regularization and positivity is discussed.
• Specific ratios of moments of valence dPDF to product of moments of valence sPDFs follow from the applied model. Such quantities, invariant of the evolution scale, may be used with future lattice data to scrutinize non-perturbative models of the pion structure.