Automated extraction of Collins–Soper kernel from lattice QCD using an autonomous AI physicist system

Understanding the three-dimensional structure of nucleons is a central goal of modern hadronic physics. Transverse-momentum-dependent (TMD) parton distribution functions (TMDPDFs) encode both longitudinal motion and intrinsic transverse dynamics, enabling a comprehensive description of nucleon structure [1]. The rapidity evolution of TMD observables—connecting measurements across different energy scales—is governed by the Collins-Soper (CS) kernel $ K(b_{\perp}, \mu) $, a universal quantity that bridges soft and collinear dynamics in QCD [2−4]. A precise determination of the nonperturbative behavior of the CS kernel is critical for accurate TMD factorization, resummation in QCD phenomenology, and the reconciliation of data from experiments such as semi-inclusive deep-inelastic scattering (SIDIS) and Drell-Yan production.

Lattice QCD provides a first-principles framework for nonperturbative QCD calculations, and the large-momentum effective theory (LaMET) [5−7] enables the determination of the CS kernel by relating Euclidean quasi-TMDWFs, which are computable on the lattice, to light-cone TMDs. Despite substantial recent advances [8−21], several challenges remain. First, nonlocal correlation functions decay exponentially with increasing transverse separation, which can lead to a poor signal-to-noise ratio (SNR) and unstable extractions. Second, multiple extrapolations, including the infinite-momentum, continuum, and chiral limits, together with renormalization procedures, introduce intricate systematic uncertainties that require careful treatment. Moreover, extracting the CS kernel from quasi-TMDWFs involves fitting momentum-dependent correlators across multiple lattice ensembles, demanding efficient numerical techniques and robust model validation. Therefore, lattice QCD calculations, such as Ref. [21], remain labor-intensive, with the engineering aspects of data processing, fitting, and extrapolation often requiring substantial effort. These limitations motivate the development of AI-driven tools to automate and optimize lattice QCD analyses.

In this work, we employ the AI-assisted agentic system PHYSMASTER [22] to address the challenges of CS kernel extraction. PHYSMASTER is an LLM-based multi-agent framework for theoretical and computational physics research, integrating Monte Carlo Tree Search (MCTS) for long-horizon task navigation, a layered academic data universe (LANDAU) for reliable knowledge reuse, and seamless integration of theoretical reasoning with code execution.

By applying this framework to state-of-the-art lattice data [21], we demonstrate that PHYSMASTER can automate a labor-intensive engineering workflow once the relevant physics ingredients are supplied. For the present benchmark, the agent performs ground-state fits, renormalization, large-λ extrapolation, Fourier transformation, and CS-kernel extraction within a few hours after the input correlators are prepared; the final results are consistent with traditional lattice calculations and show improved stability in the large-$ b_{\perp} $ region ($ b_{\perp} \gt 0.5 $ fm) through physics-inspired constraints. Our emphasis in this work is therefore on the automation of a physics-informed analysis pipeline rather than the autonomous discovery of the underlying formalism. This study not only advances the state of the art in CS-kernel extraction but also establishes a paradigm for AI-automated lattice QCD analyses of nonperturbative observables in QCD.

In the LaMET framework, the quasi-TMDWF $ \tilde{f}(x, b_{\perp}, \mu, \zeta_{z}) $ for a meson boosted to a large longitudinal momentum $ P^{z} $ is defined as [7, 23]:

where $ \tilde{\Phi}^{0} $ is the bare nonlocal matrix element constructed from gauge-invariant operators with staple-shaped Wilson links, $ Z_{O} $ is the renormalization factor for logarithmic divergences, and $ Z_{E} $ (derived from Wilson loops) cancels the linear divergences from the Wilson links. The rapidity scale $ \zeta_{z} = (2xP^{z})^{2} $ connects quasi-TMDWFs to physical light-cone TMDs.

The CS kernel governs the rapidity evolution of quasi-TMDWFs, with the relation between quasi-TMDWFs at different boost momenta, $ P_{1}^{z} $ and $ P_{2}^{z} $, given by [7, 14]:

Here, H is the perturbative matching kernel, and we adopt the $ b_{\perp} $-unexpanded next-to-leading-order (uNLO) matching scheme to avoid small-$ b_{\perp} $ expansion artifacts [17]. The CS kernel is extracted by extrapolating to the infinite-momentum limit ($ P^{z} \to \infty $), thereby suppressing power corrections of $ {\cal{O}}(1/(P^{z})^{2}) $ [19, 24−28]. This approach has been widely applied to extracting the CS kernel from lattice QCD [8−21].

PHYSMASTER enables end-to-end automated extraction of physical observables in lattice QCD through the three-stage workflow shown in Fig. 1, addressing several limitations of traditional manual analyses:

Figure 1.  (color online) The PHYSMASTER workflow for CS kernel extraction. The framework integrates literature retrieval, MCTS-driven exploration, hierarchical agent collaboration, and result validation.

● Pre-task stage: This stage is divided into query clarification, literature retrieval, and local library construction. More explicitly, PHYSMASTER first decomposes the extraction task—here, the determination of the CS kernel—into subtasks, including correlator fitting, renormalization, extrapolation, and kernel extraction, and identifies key constraints, such as symmetry requirements and perturbative behavior at small $ b_{\perp} $. It then builds a task-specific knowledge base that integrates prior results from traditional lattice QCD [21], LaMET formalism [5−7, 23], and perturbative QCD predictions [29−33]. In the present study, the raw correlators, the LaMET formalism, the renormalization prescription, and the functional ansätze used for tail stabilization and $ b_\perp $ reconstruction are fixed external inputs.

● Task execution stage: This stage includes MCTS-driven exploration, hierarchical agent collaboration, and parametrization/reconstruction. PHYSMASTER uses Monte Carlo Tree Search to navigate high-dimensional fitting and validation decisions given the prescribed physics inputs. A ''supervisor'' agent manages workflow progress and evaluates intermediate results, while a ''theoretician'' agent performs analytical derivations, code execution, and numerical fitting. In particular, PHYSMASTER optimizes fit windows, priors, and consistency checks, and then implements the externally supplied physics-motivated parametrizations of quasi-TMDWFs to stabilize the large-$ b_{\perp} $ region.

● Post-task stage: At this stage, PHYSMASTER validates the results by comparing the extracted CS kernels with conventional lattice results and perturbative predictions, as well as by assessing statistical uncertainties and selected modeling systematics. The current implementation does not yet fully disentangle all sources of systematic uncertainty, such as excited-state contamination, renormalization, tail modeling, Fourier reconstruction, and the $ P^z\to\infty $ extrapolation, in an independent manner.

Further details on PHYSMASTER can be found in the appendix and in Ref. [22].

In this section, we summarize the complete lattice QCD workflow executed by PHYSMASTER to extract the CS kernel from quasi-TMD wave functions of the pion. The system performs an end-to-end chain of operations, starting from Euclidean two-point correlators $ C_2 $ and Wilson-loop data provided by Ref. [21], and ending with a statistically controlled lattice determination of $ K(b_\perp,\mu) $ at fixed transverse separation. An overview of the workflow is given in Fig. 2, where $ b_{\perp}=4a $ and $ \mu=2 $ GeV are used. For illustrative purposes, we choose the C32P23 lattice ensemble. In the present demonstration, the gauge-field generation itself is not automated; PHYSMASTER starts from the precomputed correlators and Wilson-loop inputs.

Figure 2.  (color online) Automated long-chain analysis workflow executed by PHYSMASTER for CS-kernel extraction. Starting from raw lattice two-point correlators and Wilson loop data, the agent performs correlator-ratio construction, ground-state fitting, renormalization, large-λ extrapolation, Fourier transformation, and final CS-kernel extraction in an automated, physics-informed manner.

We begin by computing the ratio of nonlocal to local correlators to reduce the hadron-state dependence and remove the time dependence:

As shown in Fig. 2(a), the agent automatically determines the fitting window and performs either one-state or two-state ground-state fits to extract the bare quasi-TMDWF matrix element $ \Phi_{\rm{bare}}(z, b_{\perp}) = \lim_{t \to \infty} R(z, b_{\perp}; t) $.

To remove the linear and logarithmic divergences present in the bare matrix element, PHYSMASTER applies the following renormalization procedure:

where $ Z_E $ is the Wilson-loop renormalization factor and $ Z_O $ accounts for operator renormalization. The agent correctly identifies the $ z + 2L $ entry for each matrix element and performs pointwise renormalization, resulting in significantly improved large-z behavior compared with the bare results, as shown in Fig. 2(b) and Fig. 2(c).

For large $ \lambda = zP^z $, lattice signals suffer from exponential noise degradation. PHYSMASTER stabilizes the correlator tail using a physics-motivated parametrization [34]:

where the fitting parameters are constrained by Bayesian priors. In the present implementation, this large-λ functional form is supplied as an externally motivated ansatz rather than being autonomously discovered by the agent. The agent performs simultaneous fits to the real and imaginary parts over sliding windows at large λ, replacing noisy data with smoothly continued values, as shown in Fig. 2(d). The associated model dependence is treated as part of the systematic uncertainty in the reconstruction.

To obtain the momentum-space distributions, PHYSMASTER performs a Fourier transformation

The results in Fig. 2(e) are consistent with manual analyses, with minor differences in the real part attributed to different asymptotic-continuation strategies—typically accounted for as systematic uncertainties.

Finally, PHYSMASTER extracts the CS kernel using the LaMET relation in Eq. (2). The residual x and $ P^z $ dependence is removed by fitting the leading $ {\cal{O}}(1/(P^{z})^{2}) $ correction in the central region ($ x \in [0.3, 0.7] $) and extrapolating to $ P^{z} \to \infty $.

An important feature of the final extraction is that the momentum-space quasi-TMDWFs are fitted jointly across all available $ b_{\perp} $ values. As a result, the better-constrained small-$ b_{\perp} $ data help stabilize the large-$ b_{\perp} $ behavior through the shared fit parameters. This differs from the intermediate results shown earlier in the workflow, where each $ b_{\perp} $ is processed separately and therefore carries weaker cross-$ b_{\perp} $ constraints. The smoother behavior and smaller uncertainties in the final CS kernel should thus be understood as a consequence of the global, physics-constrained fit rather than as a simple point-by-point error reduction. At the same time, this stabilization may introduce additional model dependence, which we regard as part of the systematic uncertainty associated with the parametrization strategy.

Despite the complexity of the workflow, PHYSMASTER reproduces the benchmark CS kernel at $ b_{\perp}=4a $ and $ \mu=2 $ GeV within uncertainties, confirming the reliability of the automated analysis chain. For the present benchmark dataset, the sequence of data processing, high-dimensional fitting, Fourier reconstruction, and final kernel extraction can be completed within a few hours once the input correlators are prepared.

To address the signal loss at large $ b_{\perp} $, PHYSMASTER parametrizes the momentum-space quasi-TMDWF $ \tilde{f}(x, b_{\perp}, P^{z}) $ using a functional form constrained by both perturbative QCD and lattice data. In the present implementation, this functional form is imposed externally based on physical considerations; the agent then performs the numerical optimization and validation within this ansatz family. Specifically, we use:

The factors $ x^\alpha $ and $ (1-x)^\beta $ encode the expected small-x and large-x endpoint behavior of the parton distribution, while the exponential term suppresses unphysical oscillations at large $ b_{\perp} $, and the remaining polynomial captures subleading structure in the intermediate region.

PHYSMASTER determines the parameters through a joint fit to lattice quasi-TMDWF data over all available $ b_{\perp} $ values in the intermediate x range ($ x \in [0.3, 0.7] $), where power corrections are minimized. The small-$ b_{\perp} $ region ($ b_{\perp} \lt 0.3 $ fm) is further constrained by perturbative QCD predictions to reduce model dependence. Because the fit is global in $ b_{\perp} $, the better-determined small-$ b_{\perp} $ data help constrain the large-$ b_{\perp} $ behavior, yielding a smoother reconstruction than independent point-by-point treatments. The constrained results are shown in Fig. 3.

Figure 3.  (color online) The renormalized quasi-TMD wave function (quasi-TMDWF) at longitudinal momentum $ P^z=1.47 $ GeV. The real and imaginary parts are shown as functions of the momentum fraction x, demonstrating the stable behavior enforced by PHYSMASTER’s physics-motivated parametrization.

The final Collins–Soper kernel $ K(b_\perp, \mu=2{\rm{GeV}}) $ extracted by PHYSMASTER is displayed in Fig. 4. The result is consistent with state-of-the-art lattice QCD determinations within uncertainties. In the large-$ b_\perp $ region, up to $ 1\,{\rm{fm}} $, the physically constrained parametrization improves stability compared with direct extractions, providing robust nonperturbative constraints that are difficult to achieve with unconstrained fits. This also explains why the final kernel can exhibit smaller and smoother uncertainties than the intermediate quantities: the latter are treated separately at each $ b_\perp $, whereas the final kernel is obtained from a global fit across $ b_\perp $ values. We emphasize, however, that this reduction in statistical fluctuations is accompanied by additional model dependence from the imposed ansatz, which should be regarded as a systematic effect. Overall, once the physics inputs are specified, the fitting and extraction procedure is automated, demonstrating that PHYSMASTER can reliably implement a physically constrained lattice QCD analysis and produce results consistent with benchmark calculations, with minimal manual intervention.

Figure 4.  (color online) The Collins–Soper kernel $ K(b_\perp, \mu= $$ 2\; {\rm{GeV}}) $ is extracted by PHYSMASTER. The results provide robust nonperturbative constraints up to $ b_\perp \sim 1 $ fm, with improved stability relative to traditional lattice extractions [21].

The lattice determination of the CS kernel demonstrates that the automated PHYSMASTER framework is capable of performing reliable lattice QCD analyses in a structured and reproducible manner, offering clear advantages in efficiency, consistency, stability, and extensibility. For the present benchmark, once the input correlators and Wilson-loop data are prepared, the workflow— including data processing, high-dimensional fitting, Fourier reconstruction, and final kernel extraction—can be completed within a few hours. Meanwhile, MCTS-driven strategy selection and objective fit validation reduce subjective choices in fitting ranges and validation steps, while the physics-constrained reconstruction of quasi-TMD wave functions effectively mitigates signal loss at large transverse separation, enabling reliable nonperturbative extraction up to $ b_\perp\sim 1\;{\rm{fm}} $. The resulting CS kernel, extracted via this automated approach, can serve as critical input for global TMD fits, and the framework can be extended to additional lattice ensembles and other TMD observables with modest modifications.

We note, however, that the generation of raw lattice QCD data is not addressed in the present work, as PHYSMASTER takes existing lattice correlator and Wilson-loop data as input. In addition, while the parametrizations in Eqs. (7) and (8) are physically motivated, they are externally imposed and therefore introduce model dependence. Our current study also does not yet disentangle all sources of systematic uncertainty in a fully independent manner; in particular, the effects of excited-state contamination, renormalization, tail modeling, Fourier reconstruction, and the $ P^z\to\infty $ extrapolation are not separately quantified here. These limitations should be kept in mind when interpreting the apparent reduction in the final uncertainties.

To address these limitations, we will advance this AI-automated lattice QCD paradigm along several key directions in future work. We plan to incorporate additional lattice configurations and finer lattice spacings to further reduce statistical and systematic uncertainties, perform dedicated variations of the reconstruction ansatz and fit strategy to quantify model dependence more systematically, extend the automated framework to nucleon TMD observables for realistic phenomenological applications, and integrate neural-network-based parametrizations to capture complex nonperturbative behavior without manual model design. Most importantly, we will expand PHYSMASTER to cover a larger fraction of the lattice QCD pipeline, moving toward a more complete end-to-end automated workflow.

As an AI co-scientist, PHYSMASTER automates labor-intensive numerical routines while preserving physical interpretability, freeing researchers to focus on physical interpretation and theoretical innovation. This work establishes a reproducible and scalable blueprint for AI-assisted studies of parton structure and other nonperturbative observables in lattice QCD.

PHYSMASTER is a multi-agent system designed to assist research in theoretical and computational physics [22]. The system integrates theoretical reasoning with executable numerical computation and is structured to support ultra-long-horizon scientific workflows.

At the beginning of the workflow, the clarifier transforms the original natural-language query into a structured research task. More specifically, it extracts essential information from the query and formulates it as a structured contract, including the research topic, domain, task description, expected input and output formats, and relevant physical constraints such as symmetries, conservation laws, dimensional analysis, and characteristic scales. In addition, the clarifier decomposes the problem into a sequence of executable subtasks with tractable intermediate results that can be dynamically scheduled during later stages of the workflow.

To provide reliable background knowledge for the automated research process, PHYSMASTER is equipped with a persistent knowledge infrastructure, LANDAU, the Layered Academic DAta Universe, which is named in honor of the eminent physicist Lev Landau. As shown in Fig. A1, the LANDAU system organizes scientific knowledge into three major components:

Figure A1.  (color online) Structure of the LANDAU layered academic data universe, consisting of a literature library, validated methodologies, and high-confidence physics priors that support reliable AI-automated theoretical and computational physics research.

where $ {\cal{L}} $ denotes the library of knowledge extracted from retrieved scientific papers, $ {\cal{M}} $ represents the validated methodology consisting of effective reasoning paths and technical workflows, and $ {\cal{P}} $ contains high-confidence priors manually selected from textbooks or authoritative sources. Together, these layers provide both conceptual guidance and quantitative references for subsequent reasoning and computation.

The core problem-solving process is performed during the task-execution stage using a hierarchical agent collaboration framework combined with Monte Carlo Tree Search (MCTS), inspired by Ref. [35].

Real research problems in theoretical and computational physics typically require extensive iteration involving modeling, derivation, coding, and verification. Such problems often require explicit progress tracking and the exploration of multiple candidate solution trajectories, as shown in Fig. A2. In PHYSMASTER, each node in the MCTS search tree represents an attempt to construct or refine a partial solution. Subtasks generated by the clarifier are assigned to these nodes, together with concise summaries of prior exploration results and relevant background knowledge retrieved from LANDAU. The search process expands multiple trajectories in parallel, enabling the exploration of diverse modeling strategies or computational implementations.

Figure A2.  (color online) Monte Carlo Tree Search (MCTS) Exploration in PHYSMASTER.

During this stage, two agents collaborate interactively:

● The ''Supervisor'' serves as both a scheduler and a critic. It manages global scheduling, assigns subtasks for execution at each node, and provides the necessary background knowledge from LANDAU. Meanwhile, it evaluates intermediate results using evidence retrieved from LANDAU, assigns scalar rewards to nodes, summarizes the current exploration state, and generates actionable critiques that guide subsequent refinement steps.

● The ''Theoretician'' is responsible for executing its assigned subtask. Depending on the nature of the task, it may construct theoretical models, perform analytical reasoning, derive equations, or translate the model into executable code for numerical computation.

Reliable feedback is essential for guiding long-horizon exploration. The signals provided by the supervisor are fed into the tree policy to determine future node selection and expansion. Nodes are selected according to the UCT (Upper Confidence Bounds applied to Trees) criterion

where $ Q_v $ and $ N_v $ denote the accumulated reward and visit count of node v, respectively, and C controls the trade-off between exploration and exploitation. Through this structured exploration process, the system gradually identifies high-quality solution trajectories while continuing to explore alternative approaches.

In conclusion, PHYSMASTER integrates structured task formulation, a layered scientific knowledge infrastructure, and MCTS-based hierarchical reasoning to support ultra-long-horizon research workflows that combine theoretical reasoning with executable computation. This architecture enables AI-assisted exploration across research domains that rely on rigorous formal theoretical analysis and advanced numerical computation, including high-energy and particle theory, condensed matter theory, cosmology and astrophysics, and quantum information. In turn, PHYSMASTER establishes a collaborative paradigm between physicists and AI agents: physicists provide physical insight and key theoretical inputs, while the AI undertakes much of the labor-intensive computational and analytical execution, substantially enhancing overall research efficiency.