Combining lattice QCD and phenomenological inputs on generalised parton distributions at moderate skewness

We present a systematic study demonstrating the impact of lattice QCD data on the extraction of generalised parton distributions (GPDs). For this purpose, we use a previously developed modelling of GPDs based on machine learning techniques fulfilling the theoretical requirements of polynomiality, a form of positivity constraint and known reduction limits. A special care is given to estimate the uncertainty stemming from the ill-posed character of the connection between GPDs and the experimental processes usually considered to constrain them, like deeply virtual Compton scattering (DVCS). Mock lattice QCD data inputs are included in a Bayesian framework to the prior model which is fitted to reproduce the most experimentally accessible information of a phenomenological model by Goloskov and Kroll. We highlight the impact of the precision, correlation and kinematic coverage of lattice data on GPD extraction at moderate $\xi$ which has only been brushed in the literature so far, paving the way for a joint extraction of GPDs.


Introduction
Generalised parton distributions (GPDs) were introduced in the 1990s [1][2][3] and have been since then a very active topic for theoretical and experimental studies (see e.g. the review papers [4][5][6] among others).This interest is fuelled both by the interpretation of GPDs as 3D number densities of quarks and gluons within the nucleon [7,8] and by their connection to the energy momentum tensor (EMT) [2,9].
A recent study [10] pointed out the large magnitude of uncertainty underlying the connection between GPDs and the exclusive processes usually considered to constrain these objects, such as deeply virtual Compton scattering (DVCS) [11,12], timelike Compton scattering (TCS) [13] or deeply virtual meson production (DVMP) [14,15].More precisely, the coefficient functions of these processes are not practically invertible, even when taking into account QCD corrections on a large range of energy scale, and despite the constraints brought by Lorentz covariance.The results of [10] were confirmed in the moderate ξ range in [16].In this region, which constitutes the bulk of the most precise experimental data on DVCS, QCD evolution provides little help to alleviate the so-called deconvolution problem.Recently, efforts were conducted to produce GPD models giving a better account of the uncertainty associated to the ill-defined extraction of GPDs from exclusive processes.The study [17] shows that when all the theoretical properties of GPDs (including positivity) are taken into account, the deconvolution uncertainties are mostly present in the so-called Efremov-Radyushkin-Brodsky-Lepage (ERBL) kinematic region.
The ambiguities of the DVCS, TCS and DVMP processes have led to an on-going considerable theoretical and experimental effort to characterise other exclusive processes with richer kinematic structure.The best known case is double DVCS (DDVCS) [18,19] where both the incoming and outgoing photons are deeply virtual.This additional kinematic degree of freedom allows, at leading order and in a specific kinematic range, the extraction of the "deconvoluted" GPDs.A new analysis performed in [20] suggests that DDVCS is measurable in near-future experiments.Recently, other processes have been put forward such as di-photon production [21,22], photon meson pair production [23,24], or the more general single diffractive hard exclusive processes [25,26].
In addition to the experimental inputs, a new source of information regarding GPDs has emerged through lattice QCD simulations.Indeed, new formalisms developed in [27][28][29][30] have offered the possibility to connect non-local spacelike Euclidean matrix elements computed on the lattice to lightlike ones through perturbation theory.We focus here on two different formalisms.First, the large-momentum effective theory [28], known as the "quasi-distribution" formalism, has triggered a lot of interest, both in terms of simulation and perturbative computation regarding the socalled matching kernels.GPD matching kernels have been computed [31,32] and the first lattice QCD studies of GPDs have been performed [33].
The second formalism, and the one we will be focusing on through this paper, is known as "short-distance factorisation" or "pseudodistribution" formalism [29].Based on a Lorentz decomposition of the Euclidean matrix elements, it allows the extraction of parton distributions in Ioffe time ν, that is the Fourier conjugate of the momentum fraction x.This method was first introduced in [29], and steady progresses have been done since then, both on the perturbative matching side [34][35][36], and on the lattice simulation one (see for instance [37,38]).The formalism has two advantages: working with ratios of matrix elements greatly simplifies the renormalisation procedure and allows an easier extrapolation to the continuum limit.It also presents continuous matching kernels.However, this method -like any lattice-QCD attempt to access lightcone distributions -comes with a difficulty in the actual reconstruction of the momentum dependence of the M S distribution.As only a restricted range of Ioffe times can be probed numerically with acceptable noise levels, the inversion of the Fourier transform to reconstruct an x-dependent distribution is an ill-posed inversion problem, also known as an imputation problem.Several attempts have aleady been performed to try to handle this specific ill-posed problem [39][40][41].
The combination of phenomenological and non-local spacelike lattice inputs on parton distributions has already been explored in some recent papers, for instance on the proton PDF [42] and pion PDF [43].For GPD studies, the situation is made significantly more complicated due to the higher dimensionality compared to ordinary PDFs.The recent works in [44,45] have associated lattice data with various phenomenological and experimental inputs, where the authors have mostly considered GPD lattice data at vanishing skewness.An additional difficulty in offering a framework for the joint study of experimental and lattice inputs on GPDs is the number of parameters involved.Indeed, for only two light quark flavours, [45] needs to model 20 different GPDs -two valence quark distributions, two sea quark, plus a gluon distribution, repeated for the four helicity combinations H, E, H and Ẽ.This abundance of functions to extract forced the authors to employ only basic modelling of the skewness relevant in the small ξ regime.
In the present work, we develop a different strategy to combine phenomenology and lattice data focusing on moderate skewness.In this domain, lattice computations offer the perspective of a significant reduction of the uncertainty associated to the deconvolution problem of the usually considered exclusive processes.We use a previously developed GPD model presented in [17] which offers significant flexibility precisely at moderate skewness.It is illusory in the current state of experimental and lattice data to perform a satisfactory flavour and helicity decomposition with this kind of flexible model.Therefore, instead of adjusting it on actual experimental and lattice results, we assume that the flavour and helicity decomposition has been obtained already.On the experimental side, we use the phenomenological information encoded in the Goloskokov and Kroll (GK) model [46][47][48], whereas on the lattice side, we generate mock lattice data.Some tension between lattice and experimental data is hinted at in [42] and [45], whereas [43], using short-distance factorization, states that when taking into account all sources of systematic uncertainty, lattice and experimental data are generally compatible.We therefore use this assumption in the study, which offers the possibility to merge phenomenology and lattice inputs thanks to a Bayesian reweighting procedure.In short, we probe whether experimental data which suffers from a deconvolution problem and lattice data which suffers from an imputation one can be combined advantageously.For this first study, we disregard the evolution properties of GPDs and thus the matching kernels between pseudo-Ioffe distributions and GPDs.The impact of the perturbative QCD corrections is left for future works.
In section 2, we briefly review key aspects of GPDs and their connection with Ioffe-time distributions.We also revisit the main characteristics of the flexible GPD model introduced in [17].In section 3, we show how the model is adjusted to a phenomenological parametrization so as to reproduce the main experimental features of a GPD extraction.Section 4 is dedicated to the generation of the set of mock lattice QCD data.We introduce in section 5 the reweighting procedure to combine our mock lattice QCD data to the GPD model.Finally, we discuss the results of the reweighting in section 6, and conclude.

Kinematics and GPD Modelling 2.1 Definition and properties of GPDs
GPDs can be formally defined as the Fourier transform of non-local matrix elements evaluated on a lightlike distance z − [4]: where we have restricted here our definition to the unpolarized quark GPD in a nucleon whose mass is labelled m.We also omitted the Wilson line in the definition of the operator, that needs to be added when working in gauges other than lightcone one.The average momentum of the incoming and outgoing hadronic states |p 1 and |p 2 is defined as: while x is the usual lightcone momentum fraction defined as: with k + 1,2 the incoming and outgoing momentum of the struck quark.The total four-momentum transfer squared t is given by: and the skewness ξ is defined as: GPDs have to obey several properties which play a crucial role in their modelling.First, the forward limit of GPDs is given in terms of PDFs, such that: where q(x) is the unpolarised quark PDF of flavour q and Θ the Heaviside step function.GPDs are also connected to electromagnetic form factors (EFFs) through: where F q 1 is the contribution of the quark flavour q to the Dirac EFF.GPDs have to obey the socalled polynomiality property, stating that their Mellin moments are polynomials in ξ of a given order [49,50]: where [. . .] is the floor function and mod(n, 2) is 0 for n even, and 1 otherwise.The polynomiality property is understood as a consequence of the Lorentz covariance of GPDs, and can be systematically implemented thanks to the Radon transform [51][52][53] in the double distribution formalism [50].Finally, let us highlight that GPDs are also constrained by the so-called positivity properties [54][55][56][57], the most classical one being given by: Respecting all these constraints at once represents a challenge for the phenomenology of GPDs.GPDs can also be expressed as a function of the Ioffe time parameter ν = p • z [58], which is the Fourier conjugate of the average momentum fraction x.The relation between Ĥ(ν, ξ, t) and H(x, ξ, t) is thus given by [34]: In the following, we drop the "hat" on Ĥ(ν) as the explicit dependence in ν or x is enough to indicate whether we work in Ioffe time or momentum space.
As we have mentioned before, non-local matrix elements computed on the lattice are Euclidean, so the distance z between the operators defining the GPD is spacelike (z 2 < 0) contrary to the lightlike distance z − used in the usual lightcone definition of parton distributions.As a result, Ioffe-time pseudo-distributions computed on the lattice, once the continuum limit has been appropriately taken, are functions M(ν, ξ, t, z 2 ), which can be matched perturbatively to the M S scheme to yield the Ioffe-time distributions such as H(ν, ξ, t, µ 2 ) (see [34,36]).Matching from z 2dependent pseudo-distributions to µ 2 -dependent M S distributions is only a minor correction at small Ioffe time, which we will neglect in this paper.This effect should however be properly accounted for in further studies.
In the following, for simplicity we limit the scope of our study to the singlet sea quark GPD H q(+) (x, ξ, t) = H q (x, ξ, t) − H q (−x, ξ, t).Imposing this parity property means that we will only study the imaginary part of H(ν).We also ignore the t-dependence of GPDs and the entanglement between t and ξ which splits the kinematic domain between physical and unphysical regions, to work with t = 0 all along.This choice is made to primarily focus on the deconvolution of the x and ξ dependence helped by lattice data.The implications of the t dependence are left for a future work.

GPD Modelling with artificial neural networks
Our study is based on methods recently developed in [17], where the use of artificial neural network (ANN) techniques in a direct modelling of GPDs was proposed for the first time.This new way of modelling has been designed to address the problem of model dependence and implementation of the theoretical constraints one encounters in the GPD phenomenology, but also to facilitate a future inclusion of lattice QCD information.To keep our article self-consistent and self-contained, we remind now important details on ANN GPD modelling.
Our GPD model based on ANNs significantly differs from a textbook implementation of machine learning techniques (see e.g.Ref. [59]).The reason for that is the desire to fulfil in the architecture of the neural network a set of theory-driven constraints for a valid GPD model.These constraints are among others linked to the parity properties of GPDs, the polynomiality property (8) and known limits like (6).The positivity constraints (9) are not guarantied at the level of the architecture of the network, but rather enforced numerically during the training procedure.
The model proposed in Ref. [17] is built in the double distribution space involving variables β and α.The relation with the GPD H in momentum space is given by the Radon transform: where dΩ = dβ dα δ(x − β − αξ) and |α| + |β| ≤ 1.
In (11), F (β, α, t) is called a double distribution and is the object we will model.The benefit of using the double distribution space is an automatic fulfilment of the polynomiality property by the resulting GPD H(x, ξ, t).
To achieve a satisfactory flexibility and reproduction of known limits, our double distribution model is composed of three parts: (12) Let us address successively the role of each term in (12).The first one, ensures by design the proper reduction to the forward limit and has the necessary flexibility to model the x = ξ line, which is particularly relevent for the current GPD phenomenology.The prefactor (1 − x 2 ) of F C (β, α) in ( 12) combined with 1/(1 − β 2 ) in ( 13) was introduced to facilitate the fulfilment of the positivity constraint ( 9).The forward limit (unpolarised PDF for the GPD H) is denoted by f (β) ≡ H(β, 0, 0), while h C (β, α) is a profile function generalising that proposed by Radyushkin [54].In this study it is given by: Thanks to a special design of the activation function and the use of the absolute value, the neural network ANN C (|β|, α) is even w.r.t.both β and α variables, and it vanishes at the edge of the support region |β| + |α| = 1.The symmetry in β is introduced to keep the resulting GPD an odd function of x (relevant for phenomenology of DVCS and TCS), while the symmetry in α is mandatory, as it allows fulfillment of the time reversal property, i.e. the invariance under ξ ↔ −ξ exchange.The normalisation by the integral over α, which can be done analytically, allows enforcement of the proper forward limit, while the rest of the model is typically trained to reproduce the diagonal x = ξ probed by amplitudes of processes like DVCS, TCS and DVMP.
As the term F C (β, α) was found to be tightly constrained by the necessity of reproducing both ξ = 0 and x = ξ lines, an additional term F S (β, α) was introduced.This term explicitly vanishes on the ξ = 0 and x = ξ lines, i.e. it does not contribute to the fit of F C (β, α) on the phenomenological inputs.Rather, (x 2 − ξ 2 )F S (β, α) represents a contribution that is entirely unconstrained by LO DVCS data and the knowledge of PDFs, and aims at reproducing the deconvolution uncertainty of exclusive processes.Precisely, it corresponds to a LO shadow distribution as defined and studied in [10].F S is constructed in the following way: where: Since this contribution is not constrained in the fit, the major limit on its size, aside from the maximal flexibility of the neural network, is really imposed by the positivity constraint.During training, we seek to maximise the N S normalisation factor in ( 16) within the limits allowed by positivity so as to leverage the maximal flexibility.Writing the function h S (β, α) as the difference of two different profile functions characterized by ANN S (|β|, α) and ANN S (|β|, α) ensures that F S (β, α) brings no contribution to the forward limit.The f (β) factor helps to enforce the positivity.
Finally, F D (β, α) gives the additional flexibility necessary to model the D-term, a degree of freedom of GPDs associated to the final terms in ξ n+1 in (8) and which plays a crucial role in the characterisation of partonic matter [60,61].One has : and where d i are coefficients of the expansion of Dterm into Gegenbauer polynomials, and where N is an arbitrary truncation parameter.

Representation of experimental data and estimation of uncertainties
Let us now discuss how the model encodes a representation of experimental data for processes like DVCS, TCS and DVMP which we will use as an input for the lattice QCD impact study.We stress again that experimental data for the aforementioned processes mostly probe GPDs at the x = ξ line (with some additional information on the Dterm), and that at t = 0, the forward limit, i.e. the PDF, is very well known from a wealth of inclusive and semi-inclusive processes.
The uncertainty of the model is encoded in a collection of 101 replicas.This way is convenient for propagation of uncertainty to any related quantity and for the use of Bayesian reweighting techniques.A single replica represents the outcome of the independent fit to 200 x = ξ points indicated in the previous paragraph, with a random choice of the initial parameters (weights and biases of ANNs, and normalisation parameter N S in ( 16)).To evaluate the mean and standard deviation of a quantity derived from the GPD at a specific kinematic point -which can be the GPD itself, a 3D number density, an observable, etc.one may use any statistical estimator of the mean and uncertainty of the population X of 101 values returned by replicas.In this analysis we use respectively the median and the median absolute deviation (MAD): where the factor 1.4826 is derived from the assumption of gaussianity.
We note that the population of 101 values returned by replicas may contain outliers, i.e. pathological values with respect to the other entries, distorting the evaluation of mean and uncertainty.This may happen due to instabilities in the numeric procedures, which cannot be entirely avoided due to the complexity of GPD modelling and the constraining procedure.Many methods for detecting and removing outliers were proposed in the literature, like the 3σ-method [63].Alternatively, one can chose robust estimators with respect to outliers, our motivation for using MAD.
Figure 1 shows 101 replicas of H(ν, ξ, t = 0) evaluated at three values of ξ, together with its 1σ band.For small values of ξ, the replica bundle is extremely coherent, or auto-correlated, at small Ioffe time.This is due to the fact that the positivity constraint limits considerably the freedom of the model for x > ξ, on a region that is therefore all the more extended that ξ is small.

Generating mock lattice QCD data
After describing in the previous section the fit of our flexible GPD model with phenomenological inputs, we turn to the question of the generation of plausible mock lattice QCD data.In accordance with our choice of working in the pseudo-distribution formalism, we will generate data at small Ioffe time ν.More precisely, for each value of ξ that we will consider, we generate the data in three independent batches, corresponding to the following regions in ν: We introduce these three independent batches of lattice QCD data to loosely mimic actual simulations where various sets of points with different correlations arise from varying the boosted hadron momenta and the separation between fields in the non-local operator.
We assume the correlation between different batches of data to be zero.On the contrary, we consider the existence of correlations inside the batches, characterized for commodity by a simple coefficient 0 ≤ c < 1, which is the same between any two lattice data points in the same batch, and does not vary from batch to batch.This choice of constant correlation coefficient is obviously an oversimplification, but it already allows us to get a rough estimate of the impact of correlations on the determination of GPDs.
We characterise how uncertainty grows from ν = 0 to ν = 10 by an exponential slope parameter called b. Figure 2 presents the relative uncertainty profiles which we will use in this study, given by g(ν; b, s = 5%, As can be immediately checked, this formula guarantees that when ν = ν max , the relative uncertainty is 100%, and when ν = 0, it saturates to s, here chosen as 5%.We choose the uncertainty to reach 100% at ν max = 10 to replicate the common behaviour of lattice data, which tend to exhibit uncertainty on the order of its central value around this point in Ioffe time (see, for example, Fig. 9.a. of [38]).The parameter b controls the steepness of the uncertainty increase.We will therefore use as uncertainty of our mock lattice data σ Latt (ν) ≡ g(ν; b, s = 5%, ν max = 10)μ(ν) , (21) where μ(ν) is the median of the output of our flexible GPD model at ν and the relevant value of ξ.Each of the three batches of data points contains a number of samples at Ioffe times ν i .We characterize the correlation between these points inside each batch by a coefficient c and express the associated covariance matrix We then draw the central values µ Latt i of our mock lattice data points in a normal distribution centered on the median of the output of the flexible GPD model at each ν i , with the covariance matrix Ω Latt i,i .This means that we choose our mock lattice data to be globally compatible with the central value of our flexible GPD model fitted on phenomenological inputs.Figure 3 gives an example of mock lattice data set (orange points) superposed to the replicas of our GPD model.The four panels demonstrate the effect of various combinations of correlation coefficient c and level of noise parameter b.Increasing c increases the degree to which one of the central values influences the choice of the others within a given batch in ν, while increasing the parameter b results in a data set much more concentrated around the maximum likelihood of the GPD model.

Bayesian Reweighting
Now, we would like to assess the impact of the different sets of mock lattice data generated in the previous section on the GPD model.For this purpose, for each replica R k , we introduce a goodness-of-fit estimator χ 2 k : where µ Latt i is the central value of the lattice data generated at point ν i .With this definition, χ 2 k estimates the relative compatibility of a given replica R k with the mock data within the uncertainty of the latter.Bayesian reweighting consists in affecting to each replica R k a weight ω k expressed from the goodness-of-fit χ 2 k through [64] where Z is a normalization factor such that k ω k = 1.One can further define the effective fraction of replicas R which are compatible with the new data set as where the exponentiated value is the Shannon entropy of the weight set.
In the following, we will generate several sets of mock lattice data at various values of ξ and reweight our GPD model altogether with these new data at several ξ.We call this the "multikinematic" reweighting, as opposed to the "monokinematic" reweighting where only one value of ξ is considered.As we do not model any correlation between mock lattice data at different values of ξ, the total χ 2 will then be the sum of the χ 2 evaluated at each ξ.

Reweighted statistics
Once weights ω k have been attributed to the replicas, we need to define the central value and standard deviation of a weighted distribution.At each kinematic where we want to evaluate the weighted central value and deviation, we first order the replicas monotoneously, and then define the weighted median μω as the element satisfying the condition The weighted median is the element which most accurately splits the weights evenly.The weighted standard deviation estimator σω is obtained by simply replacing the median by its weighted equivalent in (19).

Metrics of the impact of Bayesian reweighting
The effective fraction of replicas compatible with the data set τ defined in (25) tells us how constraining the new data set is compared to the prior model.It is mostly a tool to evaluate the statistical significance of the reweighted distribution: if τ is too small, the weighted mean and standard deviation should not be considered as statistically significant.However, τ brings only indirect information on the reduction of uncertainty.To characterise the latter, which is our main physical objective, we introduce the ratio Σ(y), such as At a given value of y (which can represent any kinematic, in momentum space or Ioffe time), it represents the fraction of uncertainty remaining after reweighting.If Σ(y) = 1, the uncertainty of the replica bundle at y has not decreased via reweighting.If Σ(y) < 1, some reduction of uncertainty has occurred via reweighting at that point.We further define the average retainment of uncertainty in Ioffe time as and a corresponding ratio in momentum space where we adopt a logarithmic scale, We calculate the uncertainty retainment values r ν for the region 0 ≤ ν ≤ 20 and r lnx for the region 5 × 10 −3 ≤ x ≤ 1.These two metrics, intended for ν and x spaces respectively, assign a global numerical value to the reduction of uncertainty following the introduction of the mock lattice data, which will be convenient to compare various scenarios.One should note that although we generate mock lattice data in the range 0 ≤ ν ≤ 6, our uncertainty retainment r ν metric extends up to ν = 20.We do this to assess the ability for lattice data to decrease the uncertainty of our GPD model even in Ioffe time regions where we do not expect to be able to extract statistically significant lattice data.Indeed, [43] concluded that the data at low values of Ioffe time, thanks to their smaller uncertainties, represented in effect the bulk of the constraint even at larger Ioffe times.We also do not include the entire region 0 ≤ x ≤ 1 in our metric of uncertainty retainment r lnx , integrating only from the lower bound x = 5 × 10 −3 .We choose to cap the integration with this lower bound given that Σ(x) becomes relatively constant at lower x.The inclusion of such behaviour in the integration region would completely dominate the discrimination effects at large x.
6 Results of the reweighting

Monokinematic Reweighting
Let us now apply the tools of Bayesian reweighting using our GPD model fitted on phenomenological inputs as a prior, and our mock data as the new information.As a first example, we look at monokinematic reweighting, i.e. we add mock data at a single value of ξ and measure its impact on the GPD extraction at this same value.We recall that, on average, the generated mock lattice data becomes closer to the most likely output of the prior model as b increases.As c increases, the mock lattice data remains more consistently above or below the central value of the prior model.
The result for b = 2 (high precision), c = 0 (uncorrelated data) and ξ = 0.1 is shown on figure 4, while figure 5 shows the result for ξ = 0.5 and the same other parameters.The right hand side plots show the effect of reweighting in Ioffe time: we observe a large reduction of uncertainty, which remains effective far outside the range of the data (the light orange zone) although the case ξ = 0.5 shows more fluctuations linked to the lesser coherence of the replica bundle.The average retainment of uncertainty in Ioffe time is r ν = 0.16 at ξ = 0.1, and r ν = 0.25 at ξ = 0.5.
On the other hand, the left hand side plots show a far less spectacular reduction of uncertainty in momentum space.The situation is actually inverted, with a larger reduction of uncertainty at ξ = 0.5 compared to ξ = 0.1.Indeed, the retainment of uncertainty remains very large at r lnx = 0.78 at ξ = 0.1, whereas it is r lnx = 0.54 at ξ = 0.5.The fact that the uncertainty is less reduced in momentum space compared to Ioffe time is an illustration of the imputation problem evoked in the introduction.The reconstruction of the x-dependence from limited Ioffe time data is an inverse problem which triggers a significant rise in uncertainty because the mock data produces no control over the "high-frequency" behaviour of the momentum distribution.Therefore, reweighting is twice less efficient in momentum space compared to Ioffe time at ξ = 0.5, and five times less efficient at ξ = 0.1.Why does it behave so poorly at ξ = 0.1 in momentum space?The answer comes from the extreme coherence of the replicas of the GPD model in Ioffe time in the region where the reweighting is performed.This means that the GPD model contains in the end very little information in this region, such that the reweighting cannot efficiently select features of the distribution that would make a significant difference in momentum space.
We reiterate that the origin of the large coherence is the fact that positivity constrains tightly the GPD on a large part of the x dependence at ξ = 0.1 (namely for x > 0.1).We add therefore lattice data in a region where another theoretical constraint in momentum space has already considerably reduced the freedom in modelling.Lattice data must be all the better to bear any impact that the initial modelling uncertainties are small.On the other hand, the model is much more flexible at ξ = 0.5, allowing lattice data a relatively better discriminating power in momentum space.
To compare the effect of reweighting at various values of ξ, we show in figure 6 the effective fraction of surviving replicas τ and the retainment of uncertainty in Ioffe time and momentum space as functions of ξ.As the value of ξ increases, the replica bundle decoheres and the effective fraction of replicas τ drops quickly.For ξ > 0.7, a full refit seems necessary due to the low statistics.We observe that the reduction of uncertainty is systematically better in Ioffe time compared to momentum space, as expected due to the imputation problem.Using better and uncorrelated data (b = 2, c = 0) produces generally a significant reduction of uncertainty in Ioffe time compared to the other configurations, but that is not reflected in momentum space.Overall, we observe in the This highlights that with our modelling, it is in the large ξ region, where positivity does not provide significant constraints on the GPD, that we would observe the largest effect due to the inclusion of lattice data.However, we have worked here at t = 0, where we could use the very wellknown unpolarised PDF as an input to our model, and used a simplified version of the positivity constraint neglecting the role of the GPD E for instance.In general, the behaviour of GPDs at ξ = 0 but t = 0 is one of the most crucial aspects of GPD phenomenology and can be accessed readily on the lattice.Although our current study invites to preferentially include data at large ξ, in a more complete setting, the extraction of more precise tdependent PDFs plays a major role.We note that by studying both the limit ξ = 0 and large ξ, one handles two very interesting limits of GPDs: the first is directly involved in the hadron tomography program or the spin sum rule, whereas the latter probes a very different partonic dynamic in the hadron, with a sensitivity to pairs of partons.The small but non-zero values of ξ can on the other hand be constrained rather efficiently by a combination of experimental data, positivity constraints and arguments from perturbation theory if working at a reasonably hard scale [65].It remains however interesting to verify whether lattice data at moderately small value of ξ are in good agreement with positivity constraints considering the words of caution on relying excessively on these inequalities at low renormalisation scale, as underlined for instance in [66].

Multikinematic Reweighting
We would like to explore one last aspect, namely how reweighting is propagated to other values of ξ.
So far, we have only studied the effect of reweighting exactly at the same value of ξ as the data we were simulating.We show in figure 7 the result of a reweighting where the data is added at ξ = 0.1, but we observe the impact on the GPD at ξ = 0.5, with b = 1.1 (low precision) and c = 0.5 (correlated data).With these large uncertainties, the reweighting does not give any significant reduction of uncertainty at ξ = 0.5 in Ioffe time (r ν = 0.93), and even an increase of uncertainty in momentum space (r lnx = 1.15) by smearing the distribution.
Let us now add data for ξ ∈ {0.1, 0.2, 0.3} while keeping b = 1.1 and c = 0.5.The retainment of uncertainty at ξ = 0.5 drops to r ν = 0.82 and r x = 1.0 (figure 8).If we now add data for ξ ∈ {0.1, 0.2, 0.3, 0.4, 0.5}, uncertainty retainment at ξ = 0.5 tightens to r ν = 0.65 and r x = 0.75 (figure 9).It is however not better than a simple reweighting directly at ξ = 0.5 which results in r ν = 0.58 and r x = 0.64 (see Table 1 for all results, including various combinations of uncertainty and correlation that we have not mentioned in the text or figures).This demonstrates that adding some data at a given value of ξ only produces a rather minimal effect on other values of higher ξ within our modelling of GPDs.

Conclusion
We have presented a study of the impact of mock lattice QCD data at moderate value of ξ on a GPD model.The latter is based on machine learning techniques and fitted to the forward limit and diagonal x = ξ of the phenomenological GK model, which represents the typical experimental information available on GPDs.We further enforce a positivity constraint, which considerably limits the freedom of the model in the region x > ξ.We observe as a result that our model uncertainties are largely autocorrelated in the small Ioffe-time region at small ξ, meaning that lattice data only bring minimal additional reduction of uncertainty in momentum space.We observe that the reduction of uncertainty in momentum space is systematically inferior to that in Ioffe time space, as a consequence of the inverse problem of relating the two representations of GPDs.We also observe that the addition of data at some low values of ξ impacts only minimally the GPD at higher values of ξ.However, the latter point happens in a context where the t-dependence is neglected.When restored, one recovers the t-dependent PDF at zero skewness, a quantity which is not directly constrained by experimental data.We thus expect that the impact of lattice QCD data in the low-ξ region will increase with the value of |t|.The Bayesian method employed here appears to be an adequate way to combine both experimental and lattice knowledge on GPDs, when the lattice data is globally in agreement with the prior model.However, our study highlighted the impact of correlations within the lattice data on a potential joint extraction.Real lattice data frequently present a very high degree of correlation, along with systematic effects that are only starting to be under control.This mandates a very careful treatment to avoid significant biases in the assessment of uncertainty reduction, which represents a challenge for the community.

For b → 1 ,
uncertainty increases linearly with Ioffe time, while for larger values of b, it starts with a plateau of good precision and then degrades very quickly.The two cases b = 1.1 and b = 2 presented on figure 2 represent two archetypes of lattice uncertainties, corresponding respectively to data of bad or good signal to noise ratio.

Fig. 2
Fig. 2 The relative uncertainty of the mock lattice data g(ν; b, 0.05, 10) is displayed as a function of the Ioffe time ν for b = 1.1 (green) and for b = 2 (red).
Low Correlation: c = 0, High Correlation: c = 0.5, Low Precision: b = 1.1, High Precision: b = 2, r lnx : Average uncertainty retainment in x, rν : Average uncertainty retainment in ν, τ : Effective fraction of replicas retained after reweighting.Benoît Blossier and Kostas Orginos for stimulating discussions.This work is supported in part in the framework of the GLUODYNAMICS project funded by the "P2IO LabEx (ANR-10-LABX-0038)" in the framework "Investissements d'Avenir" (ANR-11-IDEX-0003-01) managed by the Agence Nationale de la Recherche (ANR), France.HD is supported in part by the U.S. DOE Grant #DE-FG02-04ER41302 and the Gordon and Betty Moore Foundation.The work of P.S. was supported by the Grant No. 2019/35/D/ST2/00272 of the National Science Centre, Poland

Table 1
Results as a function of the reweighting parameters.