Semiparametric Estimation of Physics Parameters and Parton Distribution Functions

Introduction

The determination of Parton Distribution Functions (PDFs) relies on the rigorous comparison of proton collider data and their associated uncertainties against Standard Model (SM) predictions [1, 2, 3, 4, 5, 6]. These theoretical predictions, calculated via perturbation theory, depend on fundamental SM parameters and possibly on beyond the SM (BSM) parameters.

Conventionally, these parameters are held fixed during a PDF fit, often using values derived from independent, non-hadronic processes such as leptonic production or particle decays. However, this decoupled approach is increasingly insufficient for the precision era of the Large Hadron Collider (LHC), especially in light of its high luminosity upgrade (HL-LHC). A primary objective of the LHC physics programme is the determination of both SM and BSM parameters with high precision. Such a goal demands that PDFs are treated not as static inputs, but rather as dynamic components during the fit.

To mitigate the risk of double-counting information, a common strategy is to partition experimental data into two distinct sets: one to constrain PDFs at a fixed parameter value, and another to determine parameters using a fixed PDF set. Schematically, if $\theta$ are the parameters characterizing the PDFs and $\mathbf{c}$ is a vector of SM and/or BSM parameters that affect the observable predicted by $T$ in a given model, a PDF fit sets the parameters $\mathbf{c}$ to some value determined in a separate fit $\mathbf{c} = \mathbf{c}^{*}$ while varying only the PDF parameters $\theta$,

where $\mathcal{L}$ represents the PDF luminosity and $G$ is the forward map including the DGLAP evolution of PDFs and the partonic cross sections, which connects the PDFs to the theoretical predictions. By contrast, in traditional EFT fits of one or more parameters in $\mathbf{c}$, the PDFs are set to some input set $\theta= \theta^{*}$, which in turn are determined with some $\mathbf{c} = \mathbf{c}^{*}$, while the parameters $\mathbf{c}$ are kept free in the forward map

This methodology is effective when a process exhibits a strong sensitivity to the parameter of interest but only a weak dependence on the underlying PDFs. Classic examples include the Z boson mass $m_Z$, inferred from the position of the dilepton invariant-mass resonance peak in Drell–Yan production, and the Higgs mass, determined from the diphoton invariant-mass peak in $H \to\gamma\gamma$; in both cases the parameter is encoded in a kinematic feature that remains largely invariant to the PDF-dependent overall rate. However, many critical parameters, such as the strong coupling constant $\alpha_s$, defy such a clean separation [7, 8]. Jet cross-sections, for instance, are highly sensitive to both $\alpha_s$ and the gluon distribution, making them difficult to disentangle. In this case reliable unbiased results can only be obtained if $\alpha_s$ and the PDFs are simultaneously determined [9], as opposed to using a fixed PDF set.

A similar challenge arises in the search for new physics. If high-energy datasets, where BSM effects are most likely to manifest, are excluded from PDF fits to serve solely as BSM probes, the constraints on large-$x$ gluons and quarks are severely weakened. This leaves the community reliant on ageing datasets and the potentially unreliable extrapolation of PDF models into high-energy regimes [10, 11, 12, 13, 14, 15, 16, 17, 18]. Ideally, a global fit would incorporate all available data to constrain PDFs and physical parameters simultaneously, fully accounting for experimental, theoretical, and cross-correlation effects. Tools such as SIMUnet [19, 20], Colibri [21], the Theory Covariance Method [9, 22, 23] and HERAfitter [24] have started to go in this direction; SIMUnet in particular provides a simultaneous fit of PDFs and an arbitrary number of Wilson coefficients. Yet, such an integrated approach presents significant challenges.

Beyond the considerable logistical and computational burdens, a fundamental methodological tension exists: while (B)SM parameters take values in a relatively low-dimensional real space, PDFs are continuous functions whose determination involves complex procedures such as cross-validation and the application of tolerances. Combining these distinct statistical frameworks requires a delicate approach to ensure that the precision of parameter extraction is not compromised by the methodological “noise” inherent in PDF modelling.

Semiparametric inference [25, 26, 27] addresses the estimation of a finite-dimensional parameter of interest in the presence of an infinite-dimensional nuisance parameter. A central object in semiparametric inference is the one-step corrected estimator, which corrects an initial estimate obtained with an approximate nuisance parameter by a single update along the so-called efficient score. This is defined as the component of the score of the parameters of interest that is insensitive to variations of the nuisance parameter, and, as such, removes the first-order bias induced by errors in the estimate of the nuisance. In our setting, the parameter of interest could be SM or BSM parameters, while the PDF $f$ plays the role of the infinite-dimensional nuisance parameter.

In the remainder of this section we focus on the Standard Model Effective Field Theory (SMEFT) framework at dimension-six to parametrise BSM effects [28, 29]. This choice is motivated in part by the fact that the dependence on the Wilson coefficients factorises from the SM cross section to a good approximation, as the EFT corrections enter multiplicatively through the partonic cross section while leaving the PDF convolution structure unchanged. This factorisation, made explicit in Sect. 2, is a key property that enables the analytical development of the semiparametric method presented in this work.

Recent work has explored various strategies to disentangle PDF effects in the large-$x$ region from BSM effects parametrised via Wilson coefficients that modify the tails of the high-energy experimental distributions. These strategies include conservative fits, or simultaneous fits of PDFs and parameters. In some cases, none of these strategies is ideal: conservative fits yield less precise PDFs and simultaneous fits may be too computationally burdensome to be included in a simple SMEFT fit, especially as the current tools for simultaneous fits are linked with specific PDF parametrizations; NNPDF in the case of SIMUnet and HERAPDF in the case of XFitter.

We will show that within the semiparametric framework one can estimate the bias to fitting a PDF nonparametrically. In some cases, we can also account for the bias due to fitting the PDF with the parameter of interest set at an incorrect initial value. This can be done without the need for the complex machinery of a simultaneous fit. A correction term accounting for the use of fixed PDFs is computed a posteriori and applied directly to the Wilson coefficient fit. The framework also provides the asymptotic distribution of the debiased parameter estimate, which is used to obtain well-calibrated frequentist uncertainties for the Wilson coefficient.

The paper is organized as follows. Sect. 2 introduces a linear toy model for DIS data in the presence of a simple EFT correction. In Sect. 3, we introduce and develop the methodology of semiparametric inference in this setting, before applying it as a proof of concept to extract an EFT parameter from BCDMS pseudo-data in Sect. 4 and HERA pseudo-data in Sect. 5. We summarise and conclude in Sect. 6.

PDF Model

We start by introducing the forward map that we adopt throughout this work to map from the space of PDFs to the space of observables. We focus on Deep Inelastic Scattering (DIS) data in which the PDFs enter linearly, thus allowing a fully analytic treatment of the semiparametric method. Possible generalisations to the proton-proton case are deferred to our conclusions. Furthermore, we identify the parameter of interest as a Wilson coefficient associated with a specific energy-growing operator in the dimension-six SMEFT expansion. This choice has the advantage that the parameter dependence is entirely contained in the partonic cross section and that, to a good approximation, the dependence on the parameter is encapsulated in a multiplicative factor.

Forward Map

The theoretical prediction $T^J$ for a DIS observable $J$ at kinematic point $(x,Q^2)$ is obtained by convolving¹ a known response kernel $\sigma^J$ with the PDF $f$,

where the sum runs over the $n_{\mathrm{parton}}$ active flavours at the scale $Q$ and $\Gamma_{ij}(Q^2,Q_0^2)$ denotes the evolution kernel operator that evolves the PDF $f$ from the parameterisation scale $Q_0$ to the hard scattering scale $Q$ as determined by the DGLAP evolution equations. The response kernel $\sigma_i^J$ can be computed in perturbation theory and may furthermore depend on an EFT parameter of interest $c \in\mathbb{R}$ that we assume to be real by restricting ourselves to CP even dimension-6 operators.

To speed up the computation of the integral in Eq. (4), one can write $f_i(x,Q^2)$ on a given grid in $x$ and $Q^2$, namely $\{x_\alpha,Q_\beta^2\}$ and write it as a sum of interpolating functions times the function itself computed at the grid points,

which can be expressed with respect to the parameterisation scale by interpolating the DGLAP evolution kernel,

Doing so, the integral in Eq. (4) can be written as

The integral in the parentheses can be pre-computed before doing a fit and is called Fast-Kernel (FK) table in NNPDF jargon [30]. The forward model then takes the form

where the computation of the integral at each iteration of the fitted functions $f_i$ is substituted by a precomputed finite sum over the collection of grid points. Additionally, for the purpose of this work, we are only working with a single PDF and single observable, so we drop the indices $i$ and $J$ in the following sections for simplicity.

EFT Expansion of the Response Kernel

The response kernel $\sigma$ admits the following expansion up to quadratic order in the presence of a single dimension-six EFT parameter $c$,

where we have assumed for simplicity that the functions $k_{\mathrm{lin}}$ and $k_{\mathrm{quad}}$ are independent of $x$ and thus depend only on $Q^2$. A key feature of the decomposition in Eq. (9) is that the correction functions $k_{\mathrm{lin}}$ and $k_{\mathrm{quad}}$ are treated as independent of the PDFs. This is motivated by their interpretation as ratios of EFT to SM partonic cross sections: since they are defined at the level of hard-scattering coefficients, they carry no PDF dependence at leading order. At higher orders, evaluating such ratios numerically requires a reference PDF set, but the PDF dependence largely cancels between numerator and denominator. A detailed study of this cancellation in the context of high-energy Drell-Yan tails, presented in Appendix C of Ref. [11], found the residual PDF sensitivity to be at the permil to percent level, justifying the approximation adopted here.

As a simple toy model that permits energy-growing effects characteristic for a four-fermion operator in the SMEFT (see Ref. [31] for a more realistic study of four-fermion operators on DIS observables), we take

where $\Lambda$ denotes the order of the scale of new physics. Since we will be analysing DIS data in which the highest $Q^2$ value is attained at $Q^2=30000\ \mathrm{GeV}^2$, in our toy model we set $\Lambda=500\ \mathrm{GeV}$ such that the scale of potential new physics is well above the highest $Q^2$ bin. The dependence of the EFT coefficient functions $k_{\mathrm{lin}}$ and $k_{\mathrm{quad}}$ on $Q^2$ for this choice of $\Lambda$ is illustrated in Figure 3. Note, however, that one has freedom in defining the values of both $c$ and $\Lambda$ as they always appear as the fixed combination $c/\Lambda^2$. In the applications considered in Sect. 4 and Sect. 5, the experiments probe different ranges of $Q^2$. When generating BSM pseudo data, we therefore choose $c$ separately for each application so that the induced BSM effect is sufficiently pronounced over the relevant kinematic range to clearly illustrate its impact on the inference results.

Vector (PDF)

Substituting (9) into (4), the forward model factorises as

where

Note that the functions $k_{\mathrm{lin}}$ and $k_{\mathrm{quad}}$ entering the response kernel $\sigma$ may in general depend explicitly on both $x$ and $Q^2$. The bias-correction procedure in this work is based on the simplified model in Eq. (11). Nevertheless, the semiparametric estimation procedure described in the next section can also be applied to this more general setting. The main difference is that the corresponding nuisance tangent space and efficient influence function must be rederived under the more general model.

In the following, it will be useful to consider pseudo-data generated according to a known underlying law. Given $n$ measurements with kinematic points $\{(x_i,Q_i^2)\}_{i=1}^n$ and an experimental covariance matrix $\boldsymbol{\Sigma}\in\mathbb{R}^{n\times n}$, we denote the corresponding theory predictions in the SM by

which, following Eq. (11), receives an overall multiplicative EFT correction as

where $\odot$ denotes the elementwise product and we have defined the vector of EFT corrections as

We then generate pseudo-data $\mathbf{y}$ normally distributed around $\mathbf{T}(c,f)$ with covariance $\boldsymbol{\Sigma}$,

which we will refer to as “data” from now on.

Semiparametric Inference in the SMEFT

To estimate the Wilson coefficient $c$, we make use of a statistical approach known as semiparametric inference [32, 33], which is a general framework for estimating a finite-dimensional parameter in the presence of an infinite-dimensional nuisance parameter. In our setting, the finite-dimensional parameter of interest is the Wilson coefficient $c$, while the PDF $f$ plays the role of an infinite-dimensional nuisance parameter. The procedure requires a preliminary, possibly biased, estimate of $f$ and $c$. The preliminary estimate of $f$ could come from a standard PDF fit such as the CT global fit [34], MSHT20 [35] or NNPDF4.0 [36]. We then use the efficient influence function defined subsequently to correct the bias in the preliminary estimate of $c$. Specifically, the procedure is as follows:

1. Obtain preliminary estimate $\widehat{f}_0$ by fitting the data $\mathbf{y}$ with a standard PDF fit.

2. Obtain preliminary estimate $\widehat{c}_0$ based on the fitted $\widehat{f}_0$ and the data.

3. Compute the efficient influence function $\varphi(y;c_0,f_0)$ for $c$.

4. Construct the debiased estimator $\widehat{c}_{\mathrm{one-step}}=\widehat{c}_0+\varphi(y;\widehat{c}_0,\widehat{f}_0)$.

A key benefit of this approach is that it avoids the need to refit the PDFs in the presence of the BSM contribution: the bias correction is applied once, a posteriori, to the preliminary estimate. The resulting estimator $\widehat{c}_{\mathrm{one-step}}$ is often called the “one-step estimator” since it takes one correction step in the direction suggested by the efficient influence function. Moreover, as we will show, this procedure also yields confidence intervals for $c$ with valid frequentist coverage under suitable regularity conditions.

We begin by giving a brief introduction to semiparametric inference and then describe how to apply the method to our problem setting.

Semiparametric Inference

The goal of semiparametric inference [32, 33] is to estimate a finite-dimensional parameter $c$ when there is an infinite-dimensional nuisance parameter $f$. For example, suppose that $Y_1,\ldots,Y_n$ are an i.i.d. sample from a density $f(y)$ and that we want to estimate

The parameter $c$ is scalar but the unknown density $f$ is infinite-dimensional. A natural first choice is to estimate $f$ by some non-parametric estimator $\widehat{f}$, such as a kernel density estimator, and then form the plugin estimator

This estimator is very biased and is far from optimal. Semiparametric theory improves on this by adding an influence function correction. The resulting one-step estimator takes the form

where $\varphi$ is the efficient influence function, which we define below. This estimator has the property that $\sqrt{n}(\widehat{c}_{\mathrm{one-step}}-c)$ converges to a Normal distribution with zero mean and variance $\sigma^2=\mathbb{E}[\varphi^2(Y,f)]$. This estimator is optimal in the sense that any other well-behaved estimator $\widetilde{c}$ will have an asymptotic variance that is larger than $\sigma^2$.

Figure 1 illustrates how to construct the efficient influence function by means of the efficient score. Typically, one forms one-dimensional subfamilies $f_\epsilon$ passing through the true distribution at $\epsilon=0$. For example, one may write

Vector (PDF)

where $h\in\mathcal{H}$ is a perturbation in the space of allowed perturbative corrections $\mathcal{H}$ so that $f_\epsilon$ remains a valid density.² Each choice of function $h$ defines a different subfamily. In some problems, including our SMEFT example from Sect. 2.2, the finite-dimensional parameter of interest and the infinite-dimensional nuisance parameter enter the statistical model independently. This means that the log-likelihood $\ell$ can be written explicitly as

Taking the derivative of $\ell$ with respect to $c$ gives the score for the parameter of interest

For each parametric submodel $f_\epsilon$ through $f$, we also obtain a nuisance score

The closure of the span of all such nuisance scores over $h$ is called the nuisance tangent space and is denoted by $T_f$. The efficient score for $c$ is the component of $s_c$ orthogonal to the nuisance tangent space:

where $\Pi(s_c\mid T_f)$ denotes the $L^2(P)$ projection of $s_c$ onto $T_f$. The efficient influence function $\varphi$ is then given by

where the expectation is defined with respect to the density $f$. For a more comprehensive review of semiparametric inference and efficient influence functions, see, for example, [25, 26, 27].

Estimation Strategy

For our application, we consider a general PDF function space,

and optionally impose the normalization constraint

where $\kappa$ is a positive real constant. The constrained formulation allows known PDF sum rules to be incorporated directly, such as the valence sum rule for $V,V_{3},V_{8}$, while the unconstrained formulation allows to apply the methodology to the triplet PDFs like $T_{3}$. Based on our model given by (16) with the corresponding space

our goal is to estimate $c$ and $f$ from the data vector $\boldsymbol{y}$. One of the most common approaches is through coupled maximum likelihood estimation (MLE). That is, we seek $(\widehat{c}_{0},\widehat{f}_{0})$ maximizing the log-likelihood (up to a constant)

over all $c$ and all $f\in\mathcal{F}$. Because $\mathcal{M}$ is a large, nonparametric class, $\widehat{c}_{0}$ may have first-order bias induced by estimation of $f$. The debiased one-step estimate of $c$ is

where $\varphi$ is the efficient influence function for $c$, which we will derive in Sect. 3.3. Unlike the usual i.i.d. formulation reviewed in Sect. 3.1, our observation is the full data vector $\boldsymbol{y}$, whose entries may not be i.i.d. due to experimental and theoretical correlations. This means that the log-likelihood cannot be written as a single sum over individual observations (events). Therefore, we compute the one-step correction from the efficient influence function of the full sample $\boldsymbol{y}$ instead of as an empirical average over individual observations in Eq. (19).

Under some regularity conditions [26], the one-step estimator is asymptotically efficient. In particular, we have

as $n\to\infty$, where $n$ is the total number of elements in $\boldsymbol{y}$. The variance $\sigma^{2}$ is the semiparametric efficiency bound for the scaled estimator, which is a non-degenerate limiting constant and does not depend on $n$. Equivalently, under the same normalization, the full-sample influence-function satisfies

The factor $1/n$ appears because $\varphi(\boldsymbol{y};c,f)$ is the influence-function correction for the full data vector $\boldsymbol{y}$. This is analogous to the fact that, in the i.i.d. setting from Eq. (19), the variance of the empirical average of per-observation influence functions is the variance of a single influence function divided by $n$. Therefore, a $1-\alpha$ confidence interval for $c$ can be constructed as

where $\widehat{\sigma}=\sqrt{n\mathbb{E}[\varphi^2(y;\widehat{c}_0,\widehat{f}_0)]}$ and $z_{\alpha/2}$ is the $(1-\alpha/2)$-quantile of the standard normal distribution.

Regarding the preliminary estimators $\widehat{c}_0$ and $\widehat{f}_0$, ideally these are obtained to maximize $\ell(c,f)$ simultaneously over $f$ and $c$. Since $f$ is being treated nonparametrically, the resulting estimator $\widehat{c}_0$ will be biased and our procedure will reduce this bias. However, as discussed in the introduction, this joint optimization can be challenging in practice. In this work, we focus on the more realistic decoupled fitting procedure, where we minimize over $c$ with $f$ fixed at a preliminary estimated $\widehat{f}_0$ based on the Standard Model. Specifically, we find

In practice, this introduces some extra bias to $\widehat{c}_0$. As long as this extra bias is not too large, our procedure will still remove this bias. In Appendix A, we establish the asymptotic validity of the one-step estimator under regularity conditions that do not explicitly require the preliminary estimates to be obtained from the joint MLE. Furthermore, our applications from Sect. 4 and 5 demonstrate that the proposed estimator performs well in practice with decoupled MLE.

Efficient Influence Function

As shown in Sect. 3.2, the efficient influence function $\varphi$ is the key to constructing the debiased estimator of $c$. In this section, we derive the efficient influence function $\varphi$. For notational convenience, define the linear operators

where $T_{\mathrm{SM}}(h)$ denotes the vector of SM predictions evaluated at the kinematic points $\{(x_i,Q_i^2)\}_{i=1}^n$ with input function $h$. With this notation, the mean vector is $T(c,f)=K_cf$.

The first step to derive the efficient influence function is to compute the score functions for $c$ and $f$. The score function for $c$ is given by the partial derivative of the log-likelihood with respect to $c$:

Similarly, the score for $f$ in a perturbation direction $h$ is obtained via the parametric submodel $f_\epsilon=f+\epsilon h$:

The efficient score for $c$ is then given by the projection of $s_c$ onto the orthogonal complement of the nuisance tangent space $\mathcal{T}_f$ (i.e., the space spanned by the score for $f$), which is defined with respect to the $L_2(P)$ inner product $\langle g_1,g_2\rangle=\mathbb{E}[g_1(y)g_2(y)]$. Thus,

where $h^*$ is the perturbation of $f$ whose score is closest to the score $s_c$, i.e., the variation of $f$ that best resembles the effect of a shift in $c$ on the likelihood,

Since

the projection problem is equivalent to minimizing

over the admissible perturbation directions $h$. The set of admissible directions depends on whether the normalization constraint is imposed. In the unconstrained model, $h$ may be any element of the nuisance function space. The minimizer is therefore

where $K_c^*$ denotes the adjoint of $K_c$. In the constrained model, the normalization condition

must also hold along the parametric submodel $f_\epsilon=f+\epsilon h$. Consequently, admissible perturbations must satisfy

The constrained projection problem is thus

Introducing a Lagrange multiplier $\lambda$, we define

Having the Gâteaux derivative of $\mathcal{L}$ with respect to $h$ to vanish in every direction $g$ yields the first-order condition

where $1$ denotes the constant function equal to one. Assuming that $K_c^*\Sigma^{-1}K_c$ is invertible, the solution is

The Lagrange multiplier $\lambda$ can be found as

where the inner product $\langle\cdot,\cdot\rangle$ is defined as $\langle g,h\rangle= \int g(y)h(y)\,dy$. In both the constrained and unconstrained case, substituting $h^{*}$ into the expression for the efficient score gives

and the efficient influence function (EIF) is then given by

The one-step estimate of $c$ then follows by substituting Eq. (51) into Eq. (30).

Basis Expansion

In practice, existing PDFs have been fitted to SM predictions. In our application, to introduce BSM effects, we treat the NNPDF4.0 PDF set [36] as the ground truth, and generate pseudo-data based on this PDF set and a specified Wilson coefficient $c \ne0$. Then, as a first step in the estimation procedure, we will need to obtain some preliminary estimates of $f$ and $c$ from the pseudo-data. Fitting the PDF $f$ directly at the grid points $x_{\alpha}$ can be an ill-posed problem whenever the forward map (FK) is not full rank, in which case some combination of grid values $f(x_{\alpha})$ lies in the null space of the forward map and is therefore invisible to the data. In such cases, the inverse problem needs to be regularised by introducing, for example, $m$ B-spline basis functions [37, 38] to represent the PDF $f$. Specifically, in this work, we fit a single PDF flavor $f$ by expanding $yf(y)$ in terms of $m$ basis functions $\{B_k\}_{k=1}^{m}$ on $[0,1]$ at the scale $Q_0$,

and define the design matrix $A \in\mathbb{R}^{n \times m}$ by

where the number of basis functions $m$ is such that the design matrix $A$ is full rank. The SM prediction can then be written as

Similarly, the operators $K_c$ and $\dot{K}_c$ of Eq. (35) become matrices acting on $\beta$,

while the log-likelihood reduces to

If $f$ is required to satisfy the normalization condition $\int_{0}^{1} f(y)\,dy=\kappa$, then the coefficient vector $\beta$ must satisfy

In this finite-dimensional representation, a nuisance perturbation $h(y)$ is represented by its spline coefficient vector $\gamma\in\mathbb{R}^{m}$, such that

The scores for $c$ and $\beta$ are

The efficient score is $s_{c,\mathrm{eff}}(y)=s_c(y)-s_\beta(y,\gamma^*)$, where, in the unconstrained case,

Provided that $K_c^{\mathsf T}\Sigma^{-1}K_c$ is invertible, the solution is

In the constrained case, admissible nuisance directions must preserve the normalization of $f$, and therefore satisfy

The constrained projection problem is thus

Setting $U:=2K_c^{\mathsf T}\Sigma^{-1}K_c$ and $b:=2K_c^{\mathsf T}\Sigma^{-1}\dot{K}_c\beta$, the optimality conditions form the system

which can be solved for $\gamma^*$. The efficient score and EIF then take the same form as (50) and (51) with $f$ replaced by $\beta$:

On the computational side, the implementation is straightforward. Notice that the design matrix $A$ only needs to be computed once, and it can be directly computed using the FK table with the basis functions evaluated at the specified grid points. Specifically, we replace the PDF $f$ in Eq. (8) with $B_k(y)/y$ for each basis function $B_k$, which gives the corresponding columns of $A$. Once $A$ is constructed, the matrices $K_c$ and $\dot{K}_c$ can be obtained efficiently by element-wise multiplying the rows of $A$ with $\nu(c)$ and $\nu'(c)$, respectively. Thus, most quantities required for the computation can be formed with simple matrix operations. The main remaining computational task is to solve the least squares problem to obtain $h^*$, which can be done efficiently using standard linear algebra routines.

Application: $T_3$ from BCDMS Data

As a first example, we apply the semiparametric estimation methodology described in Sect. 3 to the case of extracting Wilson coefficients in a controlled setup based on simulated DIS data from the BCDMS experiment [39, 40] sensitive to the non-singlet triplet PDF $T_3$. We will show that our method can accurately estimate EFT parameters even when starting from a preliminary estimate based on SM PDFs. For the estimates in this section and Sect. 5, we use $m = 8$ B-spline functions, as we found that this provides a reasonable compromise between flexibility and stability.

Both applications presented in this section and in Sect. 5 are deliberately restricted to a single PDF combination and a single experimental dataset. While a realistic global PDF analysis involves simultaneously constraining 7 or 8 flavour combinations across hundreds of correlated measurements from a variety of processes, this controlled setting serves to validate the semiparametric method before extending to the global case. Crucially, the problem remains genuinely non-trivial even in this simplified setting: the PDF is still an infinite-dimensional nuisance parameter, so standard parametric inference does not apply; the EFT correction enters multiplicatively through the partonic cross section (see Sect. 2); and the one-step estimator must invert a realistic experimental covariance matrix with hundreds of correlated data points. The choice of $T_3$ as the single PDF combination is itself physically motivated — as we show below, the structure function difference $F_2^p-F_2^d$ isolates $T_3$ exactly from the remaining flavours, making it a clean, natural single-PDF observable rather than an artificial construction. Single-flavour determinations of this kind have historically served as a testbed for new methodologies before their extension to global analyses: the NNPDF approach was first applied to the non-singlet distribution from DIS data [41], then to a complete DIS fit [42], before evolving into the current global analyses that include data from proton-(anti)proton colliders [43]. Following the same strategy, $T_3$ has been used as a case study for Gaussian-process-based PDF determinations [44], for studying the neural tangent kernel (NTK) of neural-network PDF fits [45], and for estimating the faithfulness of machine-learning approaches to inverse problems [46]. These applications are therefore preparatory to a realistic global analysis, as discussed in Sect. 6.

Pseudo-Data Generation

We consider the fixed-target neutral current (NC) structure functions $F_2^p$ and $F_2^n$ as measured by the BCDMS experiment. The proton structure function $F_2^p$ can be described at NNLO in QCD by a sum of convolutions of the $g$, $\Sigma$, $T_3$, $T_8$ and $T_{15}$ PDF flavours in the

evolution basis [44],

Similarly, assuming isospin symmetry between the up and down quarks, the neutron structure function $F_2^n$ can be written as

which differs from the expression for $F_2^p$ in Eq. (68) by an opposite sign in front of the coefficient function $C_{T_3}$ due to the asymmetry between the $u^+$ and $d^+$ distributions entering $T_3$. For a deuterium target, being a bound state of a proton and neutron, we therefore find that the contribution coming from $T_3$ to its structure function cancels, leaving

The difference between Eq. (70) and (68) thus isolates $T_3$,

making it an ideal proving ground to test our approach in a simple setting where only a single PDF combination enters. Indeed, the corresponding SM predictions $\mathbf{T}_{\mathrm{SM}}(f)$ can be obtained by multiplying $T_3$ with the FK tables in the evolution basis,

and the EFT predictions follow straightforwardly in our model as

We provide in Table 1 an overview of the BCDMS measurements of $F_2^p$ and $F_2^n$ we include in our fit, indicating for each its publication reference, number of data points, Björken $x$ range, $Q^2$ range, and finally the source of the corresponding theoretical prediction. Both structure functions $F_2^p$ and $F_2^d$ are measured at Björken $x$ values between $x = 0.07$ and $x = 0.75$, with corresponding $Q$ values ranging between $Q = 2.7$ GeV and $Q = 15.1$ GeV, and $Q = 3.0$ GeV and $Q = 15.1$ GeV for the proton and deuteron, respectively. Focusing on their overlap in the $(x,Q^2)$ plane to construct the observable $F_2^p - F_2^d$, and applying the NNPDF4.0 data set selection criteria [36], we find that 254 data points remain. We show the corresponding kinematic coverage in Figure 2. Theoretical predictions are computed with YADISM [48] at NNLO in QCD.

Vector (PDF)

Obs.	Ref.	$N_{\mathrm{dat}}$	$x$ range	$Q$ range (GeV)	Theory
$F_{2}^{p}$	[39]	351	$0.07 \le x \le0.75$	$2.7 \le Q \le15.1$	YADISM [48]
$F_{2}^{d}$	[40]	254	$0.07 \le x \le0.75$	$3.0 \le Q \le15.1$	YADISM [48]

Next, we generate pseudo-data $\mathbf{y}$ from a known underlying law by adding experimental Gaussian noise to the theoretical prediction $\mathbf{T}(c,f)$:

where the covariance matrix $\boldsymbol{\Sigma}$ is the experimental covariance of the difference $F_2^p - F_2^d$,

In the remainder of this section, we adopt pseudo-data generated from Eq. (74) at $c = 500$ and $\Lambda= 500$ GeV, including both linear and quadratic EFT corrections. This seemingly high value of $c$ is chosen in order to induce a sufficiently pronounced BSM effect given our default value of $\Lambda$. However, as argued in Sect. 2.2, $c$ and $\Lambda$ enter only through the combination $c/\Lambda^{2}$, so this choice is equivalent, for example, to setting $\Lambda= 100$ GeV and $c = 20$.

Results

We now present the main results of applying semiparametric inference in the context of BCDMS data. The fitting procedure is as described in Sect. 3, while the pseudo-data is generated as described in Sect. 4.1.

Figure 4 presents the results of the fit to $F_{2}^{p} - F_{2}^{d}$ at the data level, with each panel corresponding to a different $Q^{2}$ bin. The preliminary fitted curve (orange) closely follows the true BSM curve (blue), particularly in the high-$Q^{2}$ bins, where the BSM contribution is more pronounced. This behavior is not surprising since the preliminary estimates are obtained by fitting data that already contain the BSM contribution. Consequently, the preliminary fit can absorb part of the BSM-induced structure through the fitted PDF and reproduce the overall shape of the observed data. However, good agreement at the data level does not necessarily imply an accurate estimate of the parameter $c$. In particular, when the PDF is estimated under the SM-based preliminary fitting procedure, part of the BSM contribution can be absorbed into the estimated PDF, leading to bias in the resulting estimate of $c$.

Indeed, the preliminary estimate $\widehat{c}_{0}$ is biased and systematically underestimates the true value of $c = 500$. The lower orange line in Figure 5 shows $\widehat{c}_{0}$ and its corresponding 95% confidence interval for $c$. The confidence interval is computed using the parametric bootstrap method with 500 bootstrap samples. We estimate the standard error of $\widehat{c}_{0}$ by the standard deviation of the bootstrap estimates of $c$. The bootstrap confidence interval is then given by $\widehat{c}_{0} \pm z_{\alpha/2}\widehat{\sigma}$, where $\widehat{\sigma}$ is the estimated standard error. As shown in Figure 5, the resulting interval is shifted below the true value and fails to cover $c = 500$. In contrast, the one-step estimator $\widehat{c}_{\mathrm{one-step}}$ is much closer to the true value. Its confidence interval is constructed using the asymptotic normality of the one-step estimator; see Eq. (33). The resulting 95% confidence interval successfully covers the true value $c = 500$. For comparison, we also compute the joint MLE by simultaneously optimizing over $c$ and the PDF rather than first estimating the PDF through the preliminary fitting procedure. In this case, the joint MLE point estimate is close to the one-step estimate, and its corresponding 95% confidence interval, constructed using the asymptotic normality of the MLE with the

Vector (PDF)

standard error estimated from the inverse Hessian of the negative log-likelihood, also covers the true value.

Furthermore, we perform a simulation study in which the full data-generation and estimation procedure is repeated 1000 times. For each simulated dataset, we compute the preliminary estimate $\widehat{c}_0$, the one-step estimate $\widehat{c}_{\mathrm{one-step}}$, and the joint MLE. Figure 6 shows the resulting sampling distributions. The sampling distribution of $\widehat{c}_0$ is centered below the true value of $c$, confirming that the preliminary estimator is systematically biased downward. In contrast, the sampling distribution of $\widehat{c}_{\mathrm{one-step}}$ is centered much closer to the true value, which shows that the one-step correction effectively removes most of this bias. The sampling distribution of the joint MLE is likewise centered near the true value and closely resembles that of the one-step estimator.

Vector (PDF)

Figure 7 summarizes the empirical coverage and widths of the confidence intervals for $c$ over the 1000 simulations. The coverage rate is estimated as the proportion of simulated confidence intervals containing the true value $c = 500$, while the interval width is defined as the difference between its upper and lower endpoints. The confidence intervals based on the preliminary estimator $\widehat{c}_0$ have coverage substantially below the nominal 95% level, consistent with the systematic bias observed in its sampling distribution. In comparison, the intervals based on the one-step estimator achieve coverage much closer to the nominal level. The joint-MLE intervals similarly provide substantially improved coverage relative to those based on the preliminary estimator. In terms of width, we can see that the confidence intervals based on $\widehat{c}_{\mathrm{one-step}}$ are wider than those based on $\widehat{c}_0$. This increase should not be interpreted as a deficiency of the one-step procedure. The preliminary intervals are narrower in part because they do not adequately account for the uncertainty associated with estimating the nuisance parameter, which contributes to their substantial undercoverage. The one-step correction incorporates the effect of nuisance estimation through the influence-function correction, resulting in uncertainty quantification that more accurately reflects the uncertainty of the estimator. In this example, the one-step intervals are narrower than the corresponding joint-MLE intervals while achieving coverage close to the nominal level. Moreover, under regularity conditions, the one-step estimator is asymptotically efficient, indicating that it has the smallest asymptotic variance among regular estimators.

Vector (PDF)

Finally, we assess the coverage calibration of the confidence intervals for $c$ by plotting the empirical coverage rate against the nominal confidence level. Figure 8 shows that the confidence intervals based on $\hat{c}_0$ are poorly calibrated, with empirical coverage rates substantially below the nominal levels. By contrast, the confidence intervals based on the one-step estimator track the diagonal much more closely, indicating substantially improved calibration that is comparable to the calibration of the joint MLE.

Vector (PDF)

Application: $V$ from HERA Data

As a second application, we consider the extraction of the Wilson coefficient $c$ from the NC structure function $F_3^{\gamma Z}(x,Q^2)$ as measured by HERA [47]. Unlike the BCDMS experiment in Sect. 4, the proton structure function $F_3^{\gamma Z}$ is predominantly sensitive to the valence quark distributions. As these are subject to the valence sum rules, the goal of this section is to illustrate how our methodology can be extended to hold also in the presence of normalization constraints.

Pseudo-Data Generation

In NC DIS, the $F_3$ structure function receives contributions from photon-$Z$ interference and $Z$ exchange and can be written in terms of electroweak parameters as [47]

where $v_e$ and $a_e$ denote respectively the vector and axial-vector couplings of the electron to the $Z$-boson with mass $M_Z$, while the kinematic prefactor $\kappa_Z$ is defined in terms of the weak mixing angle $\theta_W$ as

In the quark-parton model, which assumes gluons are not present and the longitudinal contribution $F_L=0$, the contribution $F_3^{\gamma Z}$ coming from photon-$Z$ interference can be approximated in terms of the valence up and down distributions $u_V$ and $d_V$,

where $e_u=2/3$ and $e_d=-1/3$ denote the electric charge of the up and down-type quarks, while $a_u=1/2$ and $a_d=-1/2$ denote the axial vector coupling of the up and down-type quark to the $Z$-boson, respectively. Eq. (78) thus takes the form

It then follows straightforwardly that Eq. (79) can be written in terms of the evolution basis as

which indicates how the dominant sensitivity enters through the total valence distribution $V$. Beyond leading order, the naive quark-parton model no longer holds, and Eq. (80) receives additional contributions coming from quark-gluon radiation. However, these contributions are small in comparison and its dominant sensitivity to $V$ remains. The $F_3^{\gamma Z}$ structure function therefore provides a reasonable probe to constrain $V$ also beyond leading order.

Regarding pseudo-data generation, we consider the HERA measurement of the $F_3^{\gamma Z}$ structure function from [47] consisting of 67 data points in $(x,Q^2)$ separated in 10 distinct

$Q^2$ bins ranging from $Q^2 = 1000\ \mathrm{GeV}^2$ to $Q^2 = 30000\ \mathrm{GeV}^2$, with $x$ ranging from $x = 0.013$ to $x = 0.650$. We display this kinematic coverage in Figure 2. Following the procedure of Sect. 2.2, we generate pseudo-data for this dataset with the SM prediction $T_{\mathrm{SM}}$ computed at NNLO in QCD with YADISM [48] and the EFT correction $\nu(c)$ evaluated at $c = 2$. The corresponding uncertainties follow the experimental covariance matrix of the HERA $F_3^{\gamma Z}$ measurement. As $F_3^{\gamma Z}$ is dominantly sensitive to $V$, we fit only $V$ and fix all other contributions to a pre-determined PDF set, which we take as NNPDF4.0 [36]. In addition, we fix $V$ to NNPDF4.0 outside the data region so that we only fit in a region of $x$ where the data provides sufficient sensitivity. We have verified explicitly that without this additional constraint the forward map becomes ill-conditioned and unstable.

Results

The results are similar to those obtained in Sect. 4. That is, the preliminary estimate $\hat{c}_0$ is biased and systematically underestimates the true value of $c = 2$, while the one-step estimator $\hat{c}_{\mathrm{one-step}}$ successfully corrects for this bias and yields confidence intervals with correct coverage. The data-level fit for $F_3^{\gamma Z}$ of one realization of pseudo-data is shown in Figure 9 and the corresponding 95% confidence intervals for $c$ are shown in Figure 10. As we can see, the bias in $\hat{c}_0$ is much more pronounced compared with the BCDMS case. Nevertheless, even in this more challenging setting, the one-step correction effectively removes the bias. The sampling distributions of the estimates of $c$ obtained from the separate MLE, its one-step correction, and the joint MLE are shown in Figure 11. The coverage and width of the

Vector (PDF)

Vector (PDF)

Vector (PDF)

confidence intervals for $c$ based on the three estimation methods are summarized in Figure 12, while the coverage calibration of the confidence intervals is shown in Figure 13.

Vector (PDF)

Vector (PDF)

Overall, these results demonstrate that the preliminary estimator can suffer from substantial bias because changes induced by the BSM contribution can be partially absorbed into the estimated PDF. The one-step estimator effectively corrects this bias without having to refit the PDF and produces confidence intervals with substantially improved coverage and calibration. This is crucial for making valid inference about possible BSM effects, since biased estimates or miscalibrated intervals could lead to misleading conclusions about the size or significance of new physics contributions. Importantly, the performance of the one-step estimator is also close to that of the joint MLE, both in terms of the location of the sampling distribution and the resulting uncertainty quantification. Thus, the one-step approach can recover much of the inferential benefit of joint estimation while avoiding the need to perform a full simultaneous fit over the parameter of interest and the nuisance PDF.

Discussion and Conclusion

In this work, we developed a semiparametric framework for the simultaneous estimation of a physics parameter $c$ and uncertain PDFs, and demonstrated its effectiveness as a proof of concept in a controlled idealised setting. The central ingredient is a one-step estimator that debiases the preliminary estimate of $c$ by adding a correction term constructed from the efficient influence function. This has the large practical benefit that estimates of physics parameters extracted using biased PDFs can be corrected without having to refit the PDFs, so that a single PDF can be reused across many possible SM parameters and BSM scenarios.

We first applied the proposed method to the case of BSM pseudo-data from the BCDMS experiment sensitive to the $T_3$ triplet PDF, followed by the extraction of the valence distribution $V$ from BSM pseudo-data of $F_3^{\gamma Z}$ measured at HERA. The results demonstrate that the one-step estimator effectively corrects the bias in the preliminary estimate and produces confidence intervals with correct coverage. Working with pseudo-data and known ground truth allows us to validate the methodology rigorously and establish a firm basis for applications to real experimental data.

It is worth noting that, in the present study, the joint MLE already provides an accurate estimate of $c$, as the PDF is represented by only eight basis functions, resulting in a relatively low-dimensional nuisance model. As the PDF model becomes more flexible and higher-dimensional, however, nuisance estimation becomes increasingly challenging, and the resulting estimation errors may induce non-negligible bias in the estimate of $c$. In such settings, the one-step correction is particularly valuable, as it can reduce the first-order effect of nuisance estimation errors on the target-parameter estimate. Although the correction is applied to preliminary estimates obtained from decoupled fits of $c$ and $f$ in the present study, it could also be extended to simultaneous fits in future work, if more sophisticated machinery for simultaneous fits becomes available.

More generally, the one-step correction is not tied to a particular method for constructing the preliminary estimate. The key requirement is that the data-level fit of the preliminary estimate should be consistent and sufficiently close to the data. When the preliminary estimate is already accurate, the empirical efficient influence function will be close to zero, and the one-step correction therefore has only a negligible effect. By contrast, when the preliminary estimate is biased due to nuisance PDF estimation error, the one-step correction can reduce this bias and thereby improve both estimation accuracy and the validity of subsequent inference.

Finally, this work represents a proof of concept demonstrating the promise of the semiparametric approach to the joint determination of PDFs and physics parameters. A key strength of the framework is its complete generality: it is not restricted to Wilson coefficients of an EFT, but applies equally to any physics parameter that enters the forward map as a smooth functional. From this perspective, SM precision parameters, such as the strong coupling $\alpha_s(M_Z)$, the top mass $m_t$, the $W$-boson mass $M_W$, or the weak mixing angle $\sin^2\theta_W$, are perhaps the most compelling next application, given the central role that PDF uncertainties play in limiting the precision of such determinations. Future directions also include extending the methodology to more PDF flavours and to more complex BSM scenarios involving multiple Wilson coefficients or more general EFT expansions. A further direction is the application to hadronic processes, where the cross section is bilinear in two initial-state PDFs rather than linear in a single one as in the DIS case studied here.

---

The results presented in this work can be reproduced using the code available at the following repository:

https://github.com/richardzhs/Semiparametric_Parameter_Estimation_PDF

Acknowledgments

We are grateful to the PHYSTAT community for inspiring discussions during the early stages of this work. We especially thank Nina Elmer for providing assistance with setting up the toy model in Sect. 4. J.t.H and M.U would like to thank the members of the NNPDF collaboration for insightful discussions. The work of J.t.H is supported by the Science and Technology Facilities Council (STFC) via grant awards ST/T000600/1 and ST/X000494/1. MU’s work has been supported by the European Research Council under the European Union’s Horizon 2020 research and innovation Programme (PBSP, Grant agreement n.950246) and is partially supported by the STFC consolidated grants ST/X000664/1 and UKRI4224.

Validity of the One-Step Estimator

In this section, we establish the asymptotic validity of the proposed one-step estimator. We also discuss the interpretation and practical implications of the assumptions.

Remark 1. For Assumption (A1), the quantity

represents the efficient information for $c$, which is normalized by the number of measurements. Thus, Assumption (A1) requires that this normalized information converges to a positive finite limit. This condition can be justified through more primitive assumptions on the residual efficient direction and the covariance structure. Specifically, let $\delta_i$ be the $i$-th component of $\delta$. Then

Therefore, a sufficient set of conditions for Assumption (A1) is as follows. First, suppose that the components of $\delta$ are uniformly bounded. Second, assume that the off-diagonal mass of the precision matrix is asymptotically negligible:

Under these two conditions, the off-diagonal term satisfies

Thus, the normalized efficient information is asymptotically determined by its diagonal contribution. If, in addition,

then

In other words, Assumption (A1) holds when the dependence between different measurement points is sufficiently weak, and the average diagonal contribution to the efficient information remains dominant and converges to a positive finite limit.

Remark 2. Assumption (A2) is an empirical process condition that requires the difference between the estimated influence function and the true influence function, after centering by its expectation, to be asymptotically negligible at the $n^{-1/2}$ rate. This condition is a common requirement in semiparametric estimation, which requires the class of functions $\{\varphi(\cdot;c',f'):c'\in\mathbb{R}, f'\in\mathcal{F}\}$ belongs to Donsker class, which intuitively means that the class of functions is not too complex.

Remark 3. For Assumption (A3), although the statement requires $\lvert\widehat{c}_0-c\rvert^2$ to be $o_p(n^{-1/2})$, in practice this might be relaxed. Specifically, if we have $\nu''(\widehat{c}_0)=0$, then $\lvert\widehat{c}_0-c\rvert^2$ can be dropped from the requirement. This can be satisfied for example if $\nu(c)$ is linear in $c$, i.e. BSM effect depends linearly on the Wilson coefficient. Even when $\nu(c)$ is not exactly linear—such as in our setting, where $\nu(c)$ is quadratic in $c$—the approximation can still perform well in practice, provided that the linear term dominates the contribution. Therefore, the main requirement for the preliminary estimator is $K_{\widehat{c}_0}(\widehat{f}_0-f)=o_p(n^{-1/2})$. In other words, the preliminary estimator can be relatively inaccurate in estimating the nuisance function $f$, as long as the resulting data-level fit is sufficiently close to the data.

Proof. Let $P$ denote the distribution of $\mathbf{y}$ under the true parameters $(c,f)$, and $\widehat{P}$ denote the distribution under the estimated parameters $(\widehat{c}_0,\widehat{f}_0)$. Our target parameter $c$ can be viewed as a functional of the distribution $P$, i.e. $c=\psi(P)$. By the von Mises expansion, the efficient influence function $\varphi$ satisfies

where the second line follows from the fact that $\int\varphi(x;\widehat{c}_0,\widehat{f}_0)d\widehat{P}(x)=0$. $R_2(\widehat{P},P)$ is a second-order remainder term which depends only on products or squares of differences between $\widehat{P}$ and $P$. Having $c=\psi(P)$ and $\widehat{c}_0=\psi(\widehat{P})$ as well as the definition of the one-step estimator, we have

where $\varphi=\varphi(\cdot;c,f)$ and $\widehat{\varphi}=\varphi(\cdot;\widehat{c}_0,\widehat{f}_0)$, and

Notice that $\mathbb{E}_P[\varphi(\mathbf{y};c,f)]=0$, and $G[\widehat{\varphi}-\varphi]=o_P(n^{-1/2})$ by Assumption (A2). Therefore, we have

For the remainder term, it is given by

Define $\mu(c,f)=K_cf=\nu(c)\odot T_{\mathrm{SM}}(f)$ and write

Also, denote

Then the remainder term can be written as

Taylor expanding $\mu$ around $(\widehat{c}_0,\widehat{f}_0)$ gives

Notice that

where the second equality follows from the definition of $\widehat{\delta}$ and the fact that $K_{\widehat{c}_0}h^*$ lies in the nuisance tangent space. Similarly,

since $K_{\widehat{c}_0}(f-\widehat{f}_0)$ lies in the nuisance tangent space. Therefore, the first two terms in the expansion (A.15) cancel out with $\widehat{c}_0-c$ in the expression for $R_2(\widehat{P},P)$ in Eq. (A.14), and we are left with

Also, since

if $\nu^{\prime\prime}(\widehat{c}_0)=0$, then the reminder term reduces to

Therefore, we conclude that $R_2(\widehat{P},P)=o_P(n^{-1/2})$ under Assumption (A3). This implies that $\mathbb{E}_P[\widehat{c}_{\mathrm{one-step}}-c]=o(n^{-1/2})$. Finally, since

Assumption (A1) implies

Therefore,