Determination of proton PDF uncertainties with Markov chain Monte Carlo

This paper presents a Markov Chain Monte Carlo analysis of proton parton distribution functions using diverse experimental data and next-to-next-to-leading order QCD calculations to demonstrate how the Bayesian framework provides a statistically rigorous method for determining uncertainties, handling non-Gaussianities, and establishing well-founded tolerance criteria for the Hessian method.

The Gist

Imagine the proton, the tiny particle at the heart of every atom in your body, not as a solid marble, but as a bustling, chaotic city. Inside this city, there are three main types of citizens: quarks (the up and down ones) and gluons (the glue holding them together).

Physicists call these citizens "Partons." To understand how the universe works, especially in giant particle smashers like the Large Hadron Collider (LHC), we need a map of this city. This map is called a Parton Distribution Function (PDF). It tells us: If I zoom in on a proton, what is the probability of finding a specific type of quark or gluon at a specific speed?

However, we can't see these citizens directly. We have to build our map by throwing things at the proton and watching how they bounce off. Because the data is messy and the math is incredibly complex, there is always a margin of error on our map.

This paper is about how to draw that map more accurately and, more importantly, how to measure the uncertainty of that map without getting it wrong.

The Old Way: The "Perfect Circle" Assumption

For decades, scientists used a method called the Hessian method to calculate these errors.

The Analogy: Imagine you are trying to find the lowest point in a valley (the best map). The Hessian method assumes the valley is a perfect, smooth bowl. If you are near the bottom, the shape is predictable. You can draw a perfect circle around the bottom to say, "The true answer is definitely inside this circle."

The Problem: Real life isn't a perfect bowl. Sometimes the valley is lopsided, has a weird flat spot, or has a long tail stretching out in one direction. If you force a perfect circle onto a lopsided valley, you either miss the true answer or claim you know more than you actually do. The old method also had to guess a "tolerance" (how big the circle should be), which was often just a guess.

The New Way: The "MCMC" Adventure

This paper introduces a new approach using Markov Chain Monte Carlo (MCMC).

The Analogy: Instead of assuming the valley is a perfect bowl, imagine sending out thousands of hikers (computer simulations) to explore the valley.

1. They start at different points.
2. They take steps, sometimes going up, sometimes down, based on how "good" the map looks at that spot.
3. Over time, they naturally cluster around the best areas.
4. By looking at where all the hikers ended up, you get a real, 3D picture of the valley. You can see if it's lopsided, if it has a long tail, or if it's actually two separate valleys.

This method doesn't guess the shape; it samples the shape directly. It uses a technique called Bayesian inference, which is basically a fancy way of saying: "Start with what you think you know, look at the new data, and update your belief."

What Did They Find?

The authors ran this "hiker" simulation using data from huge experiments (like HERA, LHC, and Tevatron). Here are their main discoveries, translated into everyday terms:

1. The Map is Lopsided: They found that for some parts of the proton (specifically the "valence" quarks), the uncertainty isn't a neat circle. It's skewed. The old "perfect circle" method would have given a misleading answer here. The MCMC method showed the true, messy shape.
2. No More Guessing the "Tolerance": In the old method, scientists had to arbitrarily decide how big the error circle should be. With MCMC, the hikers naturally tell you the size of the error based on the data. It's like the hikers saying, "Hey, we only went this far, so the answer is definitely within this range."
3. Handling the "Weird" Parameters: Sometimes, a parameter in the math is so weakly constrained by data that it can go wild. The old method would just freeze it in place to avoid errors. The MCMC method handled this gracefully, showing exactly how much that parameter could wiggle without breaking the physics.

Why Does This Matter?

If you are trying to find New Physics (like a new particle that could explain dark matter), you need to be sure that a weird signal you see isn't just a mistake in your map of the proton.

• If your error bars are too small (because you assumed a perfect bowl), you might think you found a new particle when it was just a calculation error.
• If your error bars are too big, you might miss a real discovery.

By using the MCMC "hiker" method, this paper provides a more honest, statistically robust map of the proton. It admits when the data is messy and doesn't force a perfect shape onto a messy reality.

The Bottom Line

Think of this paper as upgrading from a hand-drawn sketch of a city (which assumes everything is a perfect grid) to a satellite drone flyover (which captures every alleyway, hill, and weird shape).

The authors didn't just draw a new map; they built a better compass to measure how confident we can be in that map. This is crucial for the next generation of experiments where precision is everything.

Technical Summary

1. Problem Statement

Parton Distribution Functions (PDFs) are non-perturbative quantities essential for making precise theoretical predictions in high-energy physics, particularly for the Large Hadron Collider (LHC). While experimental data has improved, theoretical uncertainties often dominate the error budget. A major source of these uncertainties is the methodology used to propagate experimental errors into the PDF parameters.

The standard Hessian method, widely used for over two decades, relies on several idealized assumptions that often fail in practice:

• Gaussianity: It assumes the posterior distribution of parameters is Gaussian (quadratic $\chi^2$ expansion).
• Tolerance Criterion ($\Delta\chi^2$): It requires an arbitrary choice of a tolerance threshold ($T$) to define the uncertainty region, as the strict statistical definition (Wilks' theorem) often breaks down due to data tensions and non-Gaussianities.
• Linearity: It struggles with non-linear dependencies between parameters and observables.

The Monte Carlo (MC) replica method (e.g., NNPDF) offers an alternative but has been criticized for potential distortions in the posterior distribution, particularly for non-linear models, and often relies on neural networks which introduce training/validation splits that can inflate uncertainties.

The Goal: To perform a robust, statistically grounded determination of proton PDF uncertainties using Markov Chain Monte Carlo (MCMC) methods to directly sample the posterior probability distribution, thereby avoiding ad-hoc tolerance criteria and Gaussian approximations.

2. Methodology

A. Data and Theory Setup

• Data: The analysis utilizes a global dataset of 1,984 data points (after kinematic cuts) comprising:
  • Deep Inelastic Scattering (DIS): Combined H1/ZEUS (HERA), BCDMS, and NMC data.
  • Drell-Yan (DY) & Boson Production: $Z$ and $W^\pm$ data from Tevatron (CDF, DØ) and LHC (ATLAS, CMS, LHCb).
• Theory:
  • Calculations are performed at Next-to-Next-to-Leading Order (NNLO) in QCD.
  • DIS data uses the aSACOT-$\chi$ scheme (approximate NNLO) for heavy-quark mass effects.
  • DY data uses NLO predictions with NNLO $K$-factors.
  • Fast interpolation grids (APPLgrid/FK-tables) are used to speed up $\chi^2$ evaluations.

B. PDF Parametrization

• PDFs are parametrized at an initial scale $Q_0 = 1.3$ GeV (charm mass) using a functional form inspired by the CJ15 analysis.
• Free Parameters: 15 free parameters are varied (covering valence, sea, and gluon distributions).
• Constraints: Momentum and number sum rules are enforced. Heavy quarks (charm, bottom) are generated perturbatively.
• Special Handling: The parameter $p_4^{dv}$ (related to the down-valence distribution) is weakly constrained and exhibits non-Gaussian behavior. Instead of fixing it (as often done in Hessian fits), the authors apply a flat uniform prior to allow the MCMC to explore its full range without biasing the $\chi^2$ landscape.

C. MCMC Implementation

• Algorithm: The study employs the Adaptive Metropolis-Hastings (aMH) algorithm.
  • Phase 1: A standard Random Walk MH with a fixed covariance matrix to reach the vicinity of the minimum.
  • Phase 2: Switches to a self-learning proposal distribution where the covariance matrix is updated based on collected samples to improve mixing and reduce autocorrelation.
• Execution:
  • 36 independent chains were generated in parallel.
  • Thermalization: The first ~140,000 steps of each chain were discarded (burn-in).
  • Thinning: To remove autocorrelation, chains were thinned by a factor of $\eta = 3000$ (determined via the $\Gamma$-method).
  • Final Sample: A combined set of 4,068 uncorrelated samples (PDF sets) representing the posterior distribution.

D. Uncertainty Estimation Methods

The authors compare three methods for extracting confidence intervals from the MCMC samples:

1. $\alpha$%-symmetric: Assumes Gaussianity; calculates mean and standard deviation (90% CI).
2. $\alpha$%-asymmetric: Uses the best-fit sample as the central value and defines bounds via 10% and 90% quantiles of the ordered samples.
3. Cumulative $\chi^2$: Defines bounds based on the 90% quantile of the $\chi^2$ distribution of the samples. This method finds the absolute range of an observable consistent with the data at a 90% confidence level, regardless of the shape of the parameter distribution.

3. Key Contributions

1. Statistical Rigor: The paper demonstrates that MCMC provides a direct sampling of the true posterior distribution, naturally handling non-Gaussianities and parameter correlations that the Hessian method misses.
2. Data-Driven Tolerance: By analyzing the distribution of $\chi^2$ values from the MCMC samples, the authors determine that the 90% quantile corresponds to $\Delta\chi^2 \approx 21.4$. This provides a statistically well-founded tolerance criterion ($T^2 \approx 21$) that can be used to calibrate the Hessian method, addressing one of its long-standing arbitrary aspects.
3. Handling Weak Constraints: The study successfully incorporates a weakly constrained parameter ($p_4^{dv}$) using a prior, a scenario where traditional Hessian fits often fail or require arbitrary fixing.
4. Benchmarking: A direct comparison between the MCMC results and a standard Hessian fit using the same data and parametrization.

4. Results

• Convergence: The 36 chains converged to the same global minimum ($\chi^2_{min} \approx 2366$), confirming the stability of the fit.
• Non-Gaussianity: The analysis revealed significant non-Gaussian behavior in the marginal distributions of certain parameters (e.g., $p_4^{uv}$), particularly for valence quark distributions.
• Uncertainty Comparison:
  • For gluon and sea quarks ($\bar{d}+\bar{u}$), where distributions are nearly Gaussian, the Hessian and MCMC (Cumulative $\chi^2$) uncertainties agree well.
  • For valence quarks ($u_v, d_v$) and strange sea ($s+\bar{s}$), the Hessian method produces symmetric uncertainty bands that underestimate or overestimate the true uncertainty compared to the MCMC results.
  • The Cumulative $\chi^2$ method yields the most conservative and robust bounds, capturing the "tails" of the distribution that symmetric methods miss.
• Tolerance Validation: The MCMC-derived tolerance ($\Delta\chi^2 \approx 21.4$) aligns closely with the theoretical expectation for 15 degrees of freedom ($\approx 22.3$), validating the use of the Hessian method if this specific tolerance is applied.

5. Significance and Outlook

• Reliability: The study proves that MCMC is a superior tool for uncertainty estimation in cases where the Gaussian approximation fails, which is common in complex global fits.
• Methodological Improvement: It offers a concrete path to improve the Hessian method by replacing arbitrary tolerance choices with data-driven quantiles derived from MCMC.
• Future Applications: The authors note that while MCMC is computationally expensive, it is essential for next-generation precision physics. They highlight ongoing work applying these methods to Nuclear PDFs (nPDFs), where constraints are sparser and non-Gaussian effects are expected to be even more pronounced.
• Limitations: The current implementation struggles with multimodal landscapes (getting stuck in local minima) and scales poorly with very high-dimensional parameter spaces. Future work may explore Hamiltonian Monte Carlo (HMC) to improve efficiency.

In conclusion, this paper establishes MCMC as a robust, statistically principled framework for PDF determination, providing a critical benchmark for uncertainty quantification in the era of high-precision collider physics.