We employ {PhysMaster}, an autonomous agentic AI system integrating theoretical reasoning, numerical computation, and exploitation strategies towards ultra-long horizon automation, to tackle long-standing challenges in non-perturbative lattice analyzes, including low signal-to-noise ratio at large transverse separation, complex systematic uncertainties, and labor-intensive manual workflows. Using the extraction of the CS kernel from quasi-transverse-momentum-dependent wave functions (quasi-TMDWFs) via large-momentum effective theory (LaMET) as a showcase, we demonstrate that \textsc{PhysMaster} automates high-dimensional fitting, renormalization, continuum-chiral extrapolation, and non-perturbative reconstruction in a fully autonomous manner. This framework drastically reduces the duration of the workflow from months to hours without compromising precision, stabilizes signals in the large-$b_\perp$ region to $1~\rm fm$, and produces results consistent with perturbative QCD and state-of-the-art traditional lattice calculations. This work validates the effectiveness of physicist-AI collaboration for first-principles QCD research and establishes a generalizable, reproducible paradigm for automated studies of parton structure and other non-perturbative observables from lattice QCD.

A quantitative investigation into the modification of ground-state field structures in two-color QCD (QC$_2$D) is presented at finite chemical potential. Using lattice simulations with Wilson gauge and fermion actions, we explore the chromo-electromagnetic field strengths under varying matter densities. To ensure accurate measurements, we develop and calibrate two highly improved topological charge operators and evaluate four gradient flow actions. Our results reveal a finite-volume crossover in the regime of the anticipated phase boundary at $\mu = m_\pi/2$, with both chromo-electric and chromo-magnetic field strengths suppressed before recovering and exceeding vacuum values at higher chemical potentials. We find the difference between the squared chromo-electric and chromo-magnetic field strengths, $E^2-B^2$, to increase in magnitude monotonically with increasing chemical potential. At $a\mu=0.7$, we find an $11\%$ suppression of $E^2$, a relatively small effect. A systematic analysis using sigmoid fits of lattice simulations in the crossover regime is performed to confirm the critical chemical potential obtained from the field structure is in agreement with the phase boundary at $m_\pi / 2$. These findings provide new insight into non-Abelian ground-state vacuum field structures and offer a foundation for future studies in real QCD.

The topological susceptibility is one of the quantities that has a large discretization error, and the error can be sensitive to the choice of fermion action. We report on our results from physical point simulations with 2+1 flavor Möbius domain wall fermion at finite temperature. We also present the chiral condensate and disconnected susceptibility. The temporal lattice size is Nt=12 and 16, and the temperature range is around 140 MeV to 250 MeV for the chiral condensate and susceptibility. A coarse lattice with Nt = 10 covers up to 500 MeV to measure the topological susceptibility.

Spectral functions play a central role in the characterization of a wide range of physical systems, including strongly interacting quantum field theories and many-body systems. Their non-perturbative determination from Euclidean correlation functions constitutes a well-known ill-posed inverse problem and has motivated the development of numerous reconstruction techniques. In this work, we propose a systematic and prior-free framework to represent spectral functions using an orthogonal functional basis derived directly from the kernel of Euclidean two-point correlation functions. By differentiating and partially integrating the Euclidean correlator over restricted time intervals, we identify a set of lattice-accessible constraints and the associated basis functions. These functions can be reorganized into an orthogonal basis with respect to which the spectral function may be approximated in a controlled manner. We demonstrate, using several model spectral functions, that the proposed expansion captures global features of the spectral function and reproduces low-energy transport coefficients with good accuracy when the spectral function is sufficiently smooth. While numerical implementation requires high-precision input for the Euclidean correlator, the present framework is intended as a tool to extract robust constraints and overall structures of spectral functions, rather than as a direct reconstruction method. The approach may thus serve as a complementary ingredient or preprocessing step for existing spectral reconstruction techniques.

We present an extension of the Evolving density matrices on Qubits (E$\rho$OQ) framework that enables efficient fault-tolerant preparation of fermionic quantum states. The original method circumvents state preparation by stochastic sampling, but faces a sign problem in fermionic systems leading to a large number of circuits necessary. We resolve this by combining classical stochastic sampling with a linear combination of unitaries method that avoids the exponential circuit scaling that plagued naïve implementations. The resulting algorithm requires $\mathcal{O}(M^2)$ $R_Z$ rotations for circuit preparation, where $M$ is the number of retained basis states. We validate the method for ground and excited states in the Thirring model, including by computing two-point correlation functions relevant to scattering. In this model for fixed accuracy $\varepsilon$, $M$ is found to scale empirically as $M \propto \frac{1}{mg}\log(1/g)\log(1/m)$.

We consider Kogut-Susskind fermions (also known as staggered fermions) in a $(3+1)$-dimensional Hamiltonian formalism and examine a chiral transformation and its associated chiral anomaly. The Hamiltonian of the massless Kogut-Susskind fermion has symmetry under the shift transformations in each space direction $S_k \, (k=1,2,3)$, and the product of the three shift transformations in particular (the odd shifts in general) may be regarded as a unitary discrete chiral transformation, modulo two-site translations. The hermitian part of the transformation kernel $\Gamma = i S_1 S_2 S_3$ can define an axial charge as $Q_A = (1/2)\sum_x \chi^\dagger(x) \left(\Gamma+\Gamma^\dagger \right)\chi(x)$, which is non-onsite, nonquantized, and commutative with the vector charge, analogous to $\tilde{Q}_A = (1/2) \sum_n ( \chi^\dagger_n \chi_{n+1} + \chi^\dagger_{n+1} \chi_{n} )$ for the $(1+1)$ dimensional Kogut-Susskind fermion. However, our $Q_A$ cannot be expressed in terms of any quantized charges in a generalized Onsager algebra. Although $Q_A$ does not commute with the fermion Hamiltonian in general when coupled to background link gauge fields, we show that they become commutative for a class of $U(1)$ link configurations carrying nontrivial magnetic and electric fields. We then verify numerically that the vacuum expectation value of $Q_A$ satisfies the anomalous conservation law of axial charge in the continuum two-flavor theory under an adiabatic evolution of the link gauge field.

The nature of the finite temperature phase transition of QCD depends on the particle density and the mass of the dynamical quarks. We discuss the properties of the phase transition at high-density, considering an effective theory describing the high-density heavy-quark limit of QCD. This effective theory is a simple model in which the Polyakov loop is a dynamical variable, and the quark Boltzmann factor is controlled by only one parameter, $C(\mu,m_q)$, which is a function of the quark mass $m_q$ and the chemical potential $\mu$. The Polyakov loop is an order parameter of $Z_3$ symmetry, and the fundamental properties of the phase transition are thought to be determined by the $Z_3$ symmetry broken by the phase transition. By replacing the Polyakov loop with $Z_3$ spin, we find that the effective model becomes a three-dimensional three-state Potts model ($Z_3$ spin model) with a complex external field term. We investigate the phase structure of the Potts model and discuss QCD in the heavy-quark region. As the density varies from $\mu=0$ to $\mu=\infty$, we find that the phase transition is first order in the low-density region, changes to a crossover at the critical point, and then becomes first order again. This strongly suggests the existence of a first order phase transition in the high-density heavy quark region of QCD.

We propose a novel random matrix model that describes the QCD Kondo phase. The model correctly implements both the chiral symmetry of light quarks and the SU(2) spin symmetry of heavy quarks. We analytically take the large-N limit with N the matrix size and show that the model has three phases: the pure Kondo phase with no chiral condensate, the pure chirally broken phase with no Kondo condensate, and the coexistence phase. The model predicts that the pairing form of the Kondo condensate in the coexistence phase is significantly altered compared to the pure Kondo phase. For each phase, we rigorously derive the low-energy effective theory of Nambu-Goldstone modes and obtain compact closed expressions for the partition function with external sources. We also include a chiral chemical potential into the model and examine the vacuum structure.

We propose the implementation of two ingredients in the phenomenological applications of the unitary approach based on the $z$-expansion of hadronic form factors, commonly referred to as the Boyd-Grinstein-Lebed (BGL) $z$-expansion [1-4]. The first ingredient is the explicit addition of a unitarity filter applied to a given set of input data for the hadronic form factors. This further constraint is not usually taken into account in the phenomenological applications of the BGL $z$-expansion. We show that it follows from the equivalence between the BGL approach and the Dispersion Matrix (DM) method [5]}, which also describes hadronic form factors in a completely model-independent and non-perturbative way. The second ingredient is represented by the introduction of suitable kernel functions in the evaluation of unitarity bounds, leading to the application of multiple dispersive bounds to hadronic form factors, whenever data and/or (non-)perturbative techniques allow to do so. This idea may be useful for the investigation of many physical processes, from the analysis of the electromagnetic form factors of mesons and baryons to the study of weak semileptonic decays of hadrons. An explicit numerical application will be presented in the companion paper [6], where the effects of sub-threshold branch-cuts are analyzed.

We apply the strategy proposed in the companion paper [1] for dealing with multiple dispersive bounds, to the case of sub-threshold branch-cuts, which is a topic addressed extensively in the literature (see, e.g., Refs. [2-8]). We consider the simultaneous application of a double dispersive bound as a proper way to take into account unitarity constraints within phenomenological analyses of hadronic form factors in the presence of sub-threshold branch-cuts. Accordingly, the standard $z$-expansion of hadronic form factors, commonly referred to as the Boyd-Grinstein-Lebed approach [4, 9-11], is modified by including simultaneously the dispersive bounds related to the pair-production and to the sub-threshold regions. For the latter one the effects of above-threshold poles are described through a simple resonance model and the possible choices of the outer function outside the pair-production region are discussed. A detailed numerical analysis of the experimental data or lattice QCD results in the spacelike region for the charged kaon form factor is presented as a direct application of the procedure of double dispersive bound. The comparison with other methodologies present in literature and with the $z$-expansion based on the single, total dispersive bound clearly shows that the $z$-expansion including the double dispersive bound provides the most precise extrapolation at large momentum transfer as well as the most stable results with respect to the choice of the outer function outside the pair-production region.