Neural solutions of coupled ghost and gluon Dyson–Schwinger equations in Landau gauge Rodrigo Carmo Terin King Juan Carlos University, Faculty of Experimental Sciences and Technology, Department of Applied Physics, Av. del Alcalde de Móstoles, 28933 Madrid, Spain [email protected]
arXiv:2607.21548v1 [hep-ph] 23 Jul 2026
Abstract The coupled ghost and gluon Dyson–Schwinger equations (DSEs) of four-dimensional Landaugauge Yang–Mills (YM) theory are solved with a neural representation trained only from renormalized equation residuals. The neural and fixed-point solutions agree at the percent level and remain stable under changes of initialization, network size, integration grid, and infrared boundary condition. Variations of the three-gluon vertex model produce substantially larger effects than the neural error. The MiniMOM ultraviolet running and the sign change of the gluon Schwinger function are also reproduced within the limitations of the truncation.
Keywords: Dyson–Schwinger equations, Yang–Mills theory, Landau gauge, physics-informed neural networks, ghost propagator, gluon propagator
1
Introduction
Quantum chromodynamics (QCD) describes the strong interaction between quarks and gluons [1, 2]. At large momentum scales its coupling becomes weak and perturbation theory can be used [3, 4]. At hadronic scales, however, the theory has to be treated nonperturbatively. Confinement and dynamical chiral symmetry breaking are two fundamental properties in this regime. Their description is closely related to the correlation functions of the theory [5–8]. These correlation functions also enter continuum calculations of hadrons and their properties [9–11]. Different methods are available for the nonperturbative regime, e.g., Lattice calculations furnish a first-principles formulation on discretized space-time and have been used successfully for many quantities [12, 13]. Functional equations constitute a continuum alternative. DSEs, functional renormalization-group equations, and nPI equations can be derived from the gauge-fixed action however in practical computations, the infinite systems have to be truncated [5, 14–17]. Consequently, two different problems have to be considered: the retained equations must be solved with sufficient numerical accuracy, and the effect of the neglected or modeled correlation functions must be understood. During the last years, larger systems of functional equations have become accessible. Model input was replaced by dynamically calculated vertices and the effect of additional diagrams was investigated [17, 18]. This progress led, among other results, to a parameter-free calculation of Landau-gauge correlation functions and to a glueball spectrum in quantitative agreement with lattice simulations [18, 19]. It should be noted that some aspects which initially appear to be only technical become important at this level. Examples are spurious divergences in the gluon equation and the realization of the perturbative resummation in a truncated system [17, 20]. The Landau-gauge YM propagators provide a simple system in which these issues can be studied. The ghost and gluon equations are coupled nonlinear integral equations. In four dimensions they have to describe the infrared solution and the logarithmic ultraviolet running at the same time. The infrared solutions form a family of decoupling solutions with an endpoint scaling solution, whereas 1
four-dimensional lattice calculations yield the decoupling type [17, 21]. Quantitative results also depend on the three-gluon vertex. In reduced propagator truncations this vertex is normally modeled, and different choices can change the gluon propagator considerably [22–24]. Neural representations offer another way to solve differential and integral equations. Early approaches used trial functions in which boundary conditions were incorporated explicitly [25]. In physics-informed neural networks (PINNs) the equation residual is part of the loss function [26, 27]. Related techniques were applied to many-body wave functions, quantum-mechanical systems, and operator learning [28–30]. For functional equations, a neural representation is attractive because it is continuous and differentiable and can be inserted directly into nonlocal integrals. On the other hand, a small loss on the training points is not sufficient. The residual must also be checked with an independent quadrature, and observables obtained from the solution may require a distinct convergence analysis. Physics-informed neural solutions of Euclidean QED Dyson– Schwinger equations were studied in our one of previous works. [31]. Additionally, our recent Minkowski-space analysis showed that constraints which are useful in one regime can become inconsistent when the analytic structure changes [32]. Here the method is applied to the YM ghost and gluon propagators. A reduced one-loop truncation is used deliberately. The same equations are solved with a direct fixed-point iteration and with a neural fixed-point procedure. Our purpose is to determine whether the neural representation reaches the same solution without propagator data and to compare its numerical error with discretization and vertex-model effects. The ultraviolet behavior and the Schwinger function are considered as further tests. The coupled equations and their renormalization are introduced in Sec. 2. The direct and neural solution methods are described in Sec. 3. Numerical results, including the convergence tests and the dependence on the three-gluon vertex, are presented in Sec. 4. The results are discussed in Sec. 5, and Sec. 6 contains the conclusions. Additional neural convergence tests are collected in Appendix A, while the physical and numerical parameters used in the calculations are summarized in Appendix B.
2
Coupled ghost and gluon equations
We consider pure SU (Nc ) YM theory in four-dimensional Euclidean space. The Landau gauge is fixed by the condition ∂µ Aaµ = 0. The gluon propagator is then transverse [33]. The ghost and gluon propagators are written as G(p2 ) , p2 Z(p2 ) ab Dµν (p) = δ ab Pµν (p) 2 , p ab DG (p) = −δ ab
(1) Pµν (p) = δµν −
pµ pν . p2
(2)
The scalar part of the gluon propagator is denoted by D(p2 ) =
Z(p2 ) . p2
(3)
The solutions considered below are of the decoupling type. Thus, the ghost dressing is finite in the infrared and the gluon propagator approaches a nonzero value. It should be noted that positive dressing functions on the Euclidean momentum axis do not imply a positive spectral representation. The propagator equations are truncated at one loop and the ghost equation contains the ghost–gluon loop. In the gluon equation, the ghost and gluon loops are retained. The ghost–gluon vertex is kept at tree level and for the three-gluon vertex, the tree-level tensor is multiplied by a scalar dressing. The tadpole is absorbed into the infrared condition and the sunset and squint diagrams are not included. This system is a reduced version of the truncations discussed in Refs. [17, 22, 34]. Its size is still sufficient to test a coupled neural solution, while the
2
same integral equations can be solved independently by a conventional iteration. With x = p2 , y = q 2 , and z = (p + q)2 , the scalar equations read 1 = Ze3 + Nc g 2 Z(y)G(z)KG (x, y, z), (4) G(x) q Z h i C 1 sub 2 = Z3 + Nc g , (5) G(y)G(z)KZgh (x, y, z) + Z(y)Z(z)KZgl (x, y, z)C AAA (x, y, z) − Z(x) x q Z
where q = d4 q/(2π)4 . The angle between the external and loop momenta is denoted by u = cos θ. With the transverse projector, the kernels are R
R
1 − u2 , yz 1 − u2 KZgh (x, y, z) = , 3xz i 2(1 − u2 ) h 2 √ 2 2 . KZgl (x, y, z) = − u xy + 6u xy(x + y) + 3x + 8xy + 3y 3xyz 2 KG (x, y, z) = −
(6) (7) (8)
These expressions are equivalent to the one-loop kernels collected in Ref. [17]. For the baseline calculation, the three-gluon-vertex dressing is chosen as AAA Cbase (x, y, z) =
G(p̄2 ) p̄2 , Z(p̄2 ) p̄2 + Λ2s
p̄2 =
x+y+z , 2
(9)
with Λ2s = 1.54 GeV2 . The factor G/Z implements a simple renormalization-group improvement. Without an additional infrared damping, this form leads to a strong enhancement in the gluon loop. The rational factor is therefore included. Similar effective constructions were used in propagator calculations [17, 22]. They should not be confused with a dynamically calculated three-gluon vertex. In order to estimate the effect of this input, three further choices are used: AAA Cbare = 1, AAA Csupp =
(10) p̄2
p̄2 + Λ2s G(p̄2 ) AAA Cno = . IR Z(p̄2 )
,
(11) (12)
The bare and suppression-only models remove the G/Z factor in different ways. The last choice keeps this factor without infrared suppression and is used to test the stability of the truncated equations. The ghost equation is renormalized by a momentum subtraction. When a hard cutoff is used, the gluon equation also contains a spurious quadratic divergence. Different prescriptions for its subtraction are known [17, 20]. Here a mass counterterm is fixed by an infrared propagator condition. The equations are specified by G(xm ) = Gm ,
Z(xs ) = 1,
D(xm ) = D0 ,
(13)
where xm is an infrared point and xs is the subtraction point. Let ΣG and ΣZ denote the loop contributions in the ghost and gluon equations. The subtracted ghost equation is 1 1 = + ΣG (x) − ΣG (xm ). G(x) Gm
(14)
The two gluon conditions give Csub =
xm xs xm xs xs [ΣZ (xm ) − ΣZ (xs )] + − D−1 . xs − xm xs − x m xs − x m 0 3
(15)
For the numerical solution, it is useful to write the gluon equation for the inverse scalar propagator ΓD (x) = x/Z(x), x −1 . (16) ΓD (x) = x + x [ΣZ (x) − ΣZ (xs )] + Csub xs This form avoids a cancellation of terms proportional to 1/x in the infrared. It should be noted that the value of Csub depends on the regulator and on the quadrature. Only the renormalized propagator is compared below. The MiniMOM coupling is defined by αMM (p2 ) = α(xs )G2 (p2 )Z(p2 ),
(17)
where the finiteness of the Landau-gauge ghost–gluon vertex in Taylor kinematics is used [35–37].
3
Numerical methods
The standard parameters are α(xs ) = 0.05,
Gm = 10,
D0 = 15.54 GeV−2 ,
xm = 10−5 GeV2 ,
xs = 7720 GeV2 ,
(18)
with Nc = 3. External momenta are taken from 10−5 to 106 GeV2 . The radial integration extends from 10−8 to 108 GeV2 . After the transformation v = ln y, the radial integral is evaluated with Gauss–Legendre quadrature. For the angular integral, Z 1
du
p
1 − u2 f (u),
(19)
−1
Gauss–Chebyshev quadrature of the second kind is used. The initial grid contains 140 external points, 140 radial nodes, and 36 angular nodes. The equations are solved by an under-relaxed fixed-point iteration and logarithmic variables are used to keep the dressing functions positive. If (n) (n) the right-hand sides at iteration n give GDSE and ZDSE , the update is (n)
ln G(n+1) = (1 − ω) ln G(n) + ω ln GDSE ,
(20)
(n) ln Z (n+1) = (1 − ω) ln Z (n) + ω ln ZDSE ,
(21)
with ω = 0.08. The iteration is stopped when the largest logarithmic change is smaller than 2 × 10−6 . The input of the network is t = ln x. It has two outputs, NG (t) and NZ (t). These functions multiply analytic baseline functions which satisfy the renormalization conditions. The baseline network has two hidden layers of width 40 and hyperbolic-tangent activation functions. Double precision is used throughout. The implementation is based on PyTorch [38]. For the ghost dressing, the standard function is Gbase (x) = 1 + (Gm − 1)
1 + (xm /sG )νG , 1 + (x/sG )νG
sG = 0.1 GeV2 ,
1 νG = , 2
(22)
and Gbase (xm ) = Gm . For the gluon dressing, we use Zbase (x) =
x/(x + m20 ) , xs /(xs + m20 )
m20 =
xs (1 − D0 xm ) , D0 xs − 1
(23)
which gives Zbase (xs ) = 1 and Zbase (xm )/xm = D0 . The neural dressing functions are defined by Gθ (x) = Gbase (x) exp [hG (t)NG (t)] ,
(24)
Zθ (x) = Zbase (x) exp [hZ (t)NZ (t)] ,
(25)
t − tm , 2
(26)
with hG (t) = tanh
4
hZ (t) = tanh
t − tm t − ts tanh , 2 2
(27)
where tm = ln xm and ts = ln xs . In this way, the three conditions in Eq. (13) are fulfilled exactly. The exponential form keeps the Euclidean dressing functions positive. As noted above, this does not impose reflection positivity or positivity of a Källén–Lehmann density. A direct minimization of the equation residuals from a generic initialization turns out to be poorly conditioned. The reason is the cancellation of large terms in the gluon equation. Therefore, a neural form of the fixed-point iteration is used first. At outer iteration n, the neural functions are inserted into the (n) (n) right-hand sides of Eqs. (14) and (16). This gives the functions GDSE and ZDSE . Mixed targets are then defined by (n)
(n)
(n)
(n)
(n)
(28)
(n)
(29)
ln Gtar = (1 − ρ) ln Gθ + ρ ln GDSE , ln Ztar = (1 − ρ) ln Zθ
+ ρ ln ZDSE ,
with ρ = 0.10. At every outer iteration, 50 Adam steps are used to project the network onto these targets. The reference setup calculation contains 70 outer iterations and the targets are generated by the equations themselves. A converged propagator solution is not used during training. After this preconditioning, the equation residuals are minimized directly. Let RG and RD denote the right-hand sides of the subtracted ghost and inverse-propagator equations. The normalized residuals are rG (x) = rD (x) =
G−1 θ (x) − RG (x)
,
(30)
x/Zθ (x) − RD (x) , 1 + 12 (|x/Zθ (x)| + |RD (x)|)
(31)
1 + 12 |G−1 θ (x)| + |RG (x)|
and the loss is
D
2 LDSE = rG
E x
D
2 + rD
E x
.
(32)
The final optimization consists of 60 Adam steps with learning rate 3 × 10−4 . The training grid has 80 external points, 72 radial nodes, and 18 angular nodes. For the validation, 120 external points, 112 radial nodes, and 30 angular nodes are used. Thus, the residuals quoted below are not evaluated with the training quadrature. The direct solution is only used after the neural calculation with the purpose to give an independent comparison. Relative errors are defined by ∥Fθ − Fdir ∥2 . ∥Fdir ∥2
ε2 [F ] =
(33)
For the ultraviolet analysis, the dressing functions are fitted with "
p2 + Λ 2 G(p ) ∝ ln Λ2
#δ
2
"
p2 + Λ 2 Z(p ) ∝ ln Λ2 2
,
#γ
,
(34)
where Λ2 = 0.04 GeV2 . In one-loop pure YM theory, δ = −9/44, γ = −13/22, and 2δ + γ = −1 for the MiniMOM coupling. The zero-spatial-momentum Schwinger function is ∆(t) =
1 π
Z pmax
dp cos(pt)D(p2 ).
(35)
0
A negative part of this function is sufficient to establish a violation of reflection positivity [17, 39, 40]. This criterion concerns the gluon two-point function, however, it is not a complete criterion for confinement.
5
(a) Ghost and gluon dressing functions.
(b) Scalar gluon propagator.
Figure 1: Direct baseline solution. The dressing functions G(p2 ) and Z(p2 ) are shown in the left panel. The right panel shows D(p2 ) = Z(p2 )/p2 . The infrared value is fixed by D(xm ) = 15.54 GeV−2 .
Figure 2: MiniMOM coupling obtained from the direct baseline solution. The coupling vanishes in the infrared, has a maximum near p2 = 0.19 GeV2 , and decreases logarithmically at large momenta.
4
Numerical results
The direct reference setup calculation converges after 165 iterations. The conditions at the infrared point are obtained as G(xm ) = 9.9999999,
D(xm ) = 15.5400000 GeV−2 .
(36)
The gluon dressing reaches its maximum, Zmax = 2.954, at p2 = 0.806 GeV2 . The gluon propagator is finite in the infrared and has a shallow maximum of 15.714 GeV−2 at p2 = 8.47 × 10−3 GeV2 . The MiniMOM coupling goes to zero at the infrared endpoint of this decoupling solution. Its max = 2.001 at p2 = 0.188 GeV2 . maximum is αMM In order to estimate the discretization error, the simulation was repeated with three grids. The errors in Table 1 are given relative to the finest result. For the baseline grid, the relative error is below 5 × 10−4 for the ghost dressing and about 3.5 × 10−3 for the gluon dressing. The larger sensitivity of the gluon equation will also be seen in the neural calculation. Starting from the analytic standard functions, the neural fixed-point iteration converges without using the direct 6
calculation
Nx
Ny
Nu
(ε2 [G], ε2 [Z])
coarse baseline fine
100 140 190
90 140 200
24 36 52
(1.44 × 10−3 , 1.06 × 10−2 ) (4.24 × 10−4 , 3.52 × 10−3 ) reference
Table 1: Direct solutions for three discretizations. The relative errors are calculated with respect to the fine solution on a common logarithmic momentum grid.
(a) Ghost dressing.
(b) Gluon dressing.
Figure 3: Direct and neural solutions for the baseline truncation. The direct solution is not used during neural training. The agreement extends over the full momentum interval. solution. For the representative baseline run, the relative errors are ε2 [G] = 2.82 × 10−3 ,
ε2 [Z] = 9.76 × 10−3 .
(37)
The largest pointwise differences are 0.71% for the ghost dressing and 2.07% for the gluon dressing. As in the direct refinement test, the gluon is more sensitive to the numerical treatment. On the independent validation grid, the largest residuals are max |rG (x)| = 8.34 × 10−4 ,
max |rD (x)| = 6.80 × 10−3 .
x
x
(38)
The corresponding mean squared values are 8.8 × 10−8 and 5.0 × 10−6 . The residual test is needed in addition to the comparison with the direct curves. Otherwise, agreement could result from interpolation without an equally accurate solution of the integral equations. Five parameter initializations were used for the width-40 network. The last layer was set to zero, so the initial dressing functions were the same, but the hidden representations were different. The resulting errors are ε2 [G] = (3.25 ± 0.77) × 10−3 ,
(39)
−3
(40)
ε2 [Z] = (9.48 ± 0.29) × 10
.
Over the full momentum interval, the largest relative standard deviation is 0.81% for the ghost and 0.75% for the gluon. Three additional runs were started from nonzero perturbations of the last layer and they converge to the same branch as well. Their ghost errors lie between 3.78 × 10−3 and 5.41 × 10−3 , and the gluon errors lie between 7.15 × 10−3 and 1.14 × 10−2 . The hidden-layer width was changed from 24 to 64. The ghost error remains close to 3 × 10−3 . For the gluon it changes from 9.06 × 10−3 to 8.19 × 10−3 and no monotonic dependence on the width is found. At this accuracy, the effects of optimization, quadrature, and initialization are of comparable size. The width scan therefore tests stability but does not provide an asymptotic 7
(a) Pointwise relative differences.
(b) Residuals on the validation grid.
Figure 4: Numerical tests of the baseline neural solution. The left panel shows the difference from the direct result. The right panel shows the ghost and inverse-gluon-propagator residuals evaluated with a finer quadrature than the one used in training.
(a) Ghost dressing.
(b) Gluon dressing.
Figure 5: Neural solutions for five parameter initializations. The equations and renormalization conditions are identical in all runs.
8
(a) Variation of the network width.
(b) Variation of the training quadrature.
Figure 6: Relative neural errors for different network widths and integration grids. The gluon solution changes more under quadrature refinement than under the variation of the width.
(a) Ghost dressing.
(b) Gluon dressing.
Figure 7: Direct solutions for three models of the three-gluon vertex. The same renormalization conditions are used in all cases. The main change occurs in the gluon dressing around its maximum. convergence law. A clearer dependence is obtained for the integration grid. Increasing the radial– angular quadrature from 48×12 to 96×24 changes the gluon error from 1.36×10−2 to 8.36×10−3 . The ghost error stays between 2.8 × 10−3 and 3.3 × 10−3 . This confirms that the gluon equation is the more demanding part of the coupled system. The largest changes are attained when the three-gluon-vertex model is varied. The bare and suppression-only models both lead to positive Euclidean solutions. Compared with the baseline calculation, the gluon dressing changes by about 33–36% in relative L2 norm. The maxima of the gluon dressing and of the MiniMOM coupling are given in Table 2. For each of the three truncations, the neural solution remains below the percent level relative to the corresponding direct solution. vertex model baseline G/Z with IR suppression bare IR suppression only
Zmax 2.954 1.602 1.510
max αMM
ε2 [Gθ ]
ε2 [Zθ ]
2.001 1.643 1.536
2.82 × 10−3
9.76 × 10−3 5.14 × 10−3 2.07 × 10−3
3.16 × 10−3 3.57 × 10−3
Table 2: Results for the three convergent vertex models. The neural errors are calculated with the direct solution of the same truncation as reference.
9
(a) Ghost dressing.
(b) Gluon dressing.
Figure 8: Direct solutions for four infrared ghost conditions. The infrared gluon propagator is kept fixed. The solutions have the same ultraviolet normalization but differ in the transition region. The model without infrared suppression does not converge to the same positive branch. Before the fixed-point iteration breaks down, the maximum of the gluon dressing grows to about 9.5 and the MiniMOM coupling reaches about 4.69. The instability occurs already in the direct calculation. It is therefore caused by the truncated equations and not by the neural representation. For the computations considered here, the vertex dependence is much larger than the direct discretization error or the neural reconstruction error. This observation is important for the interpretation of a numerically accurate neural result. A more accurate solution of a fixed operator does not reduce the uncertainty caused by the modeled vertex. The decoupling solutions of the functional equations are selected by an infrared boundary condition [17, 21]. We use G(xm ) = 5, 10, 15, 20 (41) while D(xm ) = 15.54 GeV−2 is kept fixed. The direct iteration converges for all four values. The main change in the gluon dressing occurs around its maximum. The maximum of the MiniMOM coupling increases from 1.16 to 2.58. The neural calculation was also performed for G(xm ) = 5, 15, and 20. The relative ghost errors are 2.47 × 10−3 , 3.21 × 10−3 , and 4.26 × 10−3 . The corresponding gluon errors are 7.76 × 10−3 , 7.47×10−3 , and 1.09 ×10−2 . Hence, the neural iteration is not restricted to the baseline boundary condition. It follows the different members of the decoupling family. For the ultraviolet (UV) analysis, the external momentum interval is extended to p2 = 1010 GeV2 and the radial cutoff to 1012 GeV2 . In the highest fit interval, the direct solution gives δdir = −0.2744, γdir = −0.4467. (42) The distinct values differ from −9/44 and −13/22. The reduced one-loop truncation and the effective vertex therefore do not reproduce the two anomalous dimensions separately. Their combination is 2δdir + γdir = −0.9955, (43) which is close to the one-loop value −1 of the MiniMOM coupling. For the extended neural solution, the relative errors are ε2 [G] = 7.35 × 10−3 ,
ε2 [Z] = 1.80 × 10−2 .
(44)
The combination extracted from the highest interval is 2δθ + γθ = −1.012.
10
(45)
Figure 9: Ultraviolet exponents of the extended direct solution for different lower boundaries of the fit interval. The horizontal lines show δ = −9/44 and γ = −13/22. solution
first zero t0 [GeV−1 ]
minimum of ∆(t)
7.459 7.408 7.054
−3.53 × 10−2 −3.60 × 10−2 −5.11 × 10−2
direct baseline direct fine neural, representative seed
Table 3: First zero and minimum of the gluon Schwinger function at pmax = 1000 GeV. Thus, our neural simulation reproduces both features of the direct result: the correct MiniMOM combination and the displaced individual exponents. It should be noted that this is the expected behavior of a solver which uses only the truncated equations. The missing UV structure has to be improved in the truncation and cannot be supplied by the numerical method. Finally, the Schwinger function is evaluated by a direct cosine integration. The cutoffs pmax = 100, 300, and 1000 GeV are used. The first zero is stable under this variation. Results for pmax = 1000 GeV are listed in Table 3. The two direct discretizations give first zeros which differ by less than one percent. The neural zero is shifted by about five percent relative to the fine direct result. Nevertheless, the negative part is present in all cases. Thus, the violation of reflection positivity is reproduced qualitatively. The position and depth of the negative region are more sensitive than the Euclidean dressing functions themselves.
5
Discussion
The results enable the different sources of uncertainty to be compared. The first one is the discretization of the direct integral equations. For the reference setup grid it is below the permille level for the ghost and at a few per mille for the gluon. The second one is our neural solution of these discretized equations. Its error is somewhat larger, but remains at the per-mille to percent level. The tests with different initializations, network widths, and integration grids indicate that the same fixed point is obtained. The residuals on the independent grid support this conclusion. A much larger effect is caused by the three-gluon vertex. Replacing the initial model by a bare or suppression-only dressing changes the gluon solution by more than thirty percent. Our neural solution remains accurate for each of these choices. For the present truncation, the sizes of the
11
Figure 10: Gluon Schwinger function from the baseline direct, fine direct, and neural propagators. All three functions become negative. The first zero of the neural result occurs at a smaller Euclidean time. three effects can therefore be summarized as direct discretization error ≲ neural reconstruction error ≪ three-gluon-vertex dependence. (46) This comparison also shows why the numerical and truncation errors should be kept separate. Increasing the network size cannot compensate for missing or modeled correlation functions. The UV calculation provides an example of this point. The direct truncation does not yield the correct anomalous dimensions of the two propagators separately. It does, however, reproduce their MiniMOM combination with good accuracy. Our neural simulation gives the same result. This is not a shortcoming of our neural solver, it reflects the perturbative content of the equations which are solved. The Schwinger function is more sensitive to the propagator than the pointwise comparison on the Euclidean axis. The direct and neural propagators differ by only a few percent, whereas the first zero changes by about five percent and the minimum changes more strongly. Oscillatory transforms can therefore require a higher accuracy than the one needed for the original dressing functions. The negative part itself is stable. In this respect, the result is consistent with the observation that positivity properties should be checked after the solution and not imposed in a way which removes the signal of their violation [32]. Several restrictions of the calculation should be kept in mind. The ghost–gluon vertex is bare. The three-gluon vertex contains only one modeled scalar dressing. The four-gluon vertex and the two-loop diagrams of the gluon equation are not included. The mass counterterm removes the leading cutoff dependence, but its treatment is part of the truncation. Furthermore, only Euclidean pure YM theory is considered. The results are therefore a study of the numerical solution of a fixed propagator truncation and not a new state-of-the-art calculation of YM correlation functions. The vertex test indicates the most useful next step. A neural representation of the three-gluon vertex could be coupled to its DS or 3PI equation. This would reduce the model input in the propagator system. It would also provide a more demanding test of the neural procedure because a three-point function depends on several momentum variables. A parameter-dependent network for the infrared family is another possible extension. Complex momenta and Minkowski space require additional work. In that case, sign-indefinite spectral functions and complex singularities have to be enabled from the beginning. 12
6
Conclusions
The coupled ghost and gluon DSEs were solved with our neural representation. The calculation was performed in four-dimensional YM theory in Landau gauge. A one-loop propagator truncation with a modeled three-gluon vertex was used. The same renormalized equations were also solved by a conventional fixed-point iteration. It is important to remark that our neural computation does not use converged propagator data. Targets from the equations are employed for a fixed-point preconditioning, and the renormalized residuals are minimized afterwards. For the baseline system, the relative errors are 2.8 × 10−3 for the ghost dressing and 9.8 × 10−3 for the gluon dressing. Similar results are attained for different initializations, network widths, integration grids, and infrared boundary conditions. The residuals evaluated with a finer quadrature show that the agreement is not restricted to the training points. It turns out that the model for the three-gluon vertex is more important than the remaining numerical error. The alternative vertex models change the gluon dressing by more than thirty percent, while our neural solution of each fixed truncation stays below the percent level. The UV test gives a similar lesson. The MiniMOM combination of anomalous dimensions is obtained, but the distinct exponents are not. Our neural result follows this behavior of the truncated equations. The Schwinger function becomes negative for the direct and neural propagators, although the position of its first zero is more sensitive to small propagator differences. Neural representations can therefore be used for coupled and renormalized functional equations with an accuracy comparable to conventional numerical methods. Their usefulness does not remove the usual truncation problem. For the system considered here, an improvement of the vertex sector is more important than a further increase of the propagator-network size. A dynamical neural calculation of the three-gluon vertex is a natural next step.
A
Summary of neural convergence tests hidden width
ε2 [G]
ε2 [Z]
24 40 64
3.39 × 10−3 2.82 × 10−3 2.99 × 10−3
9.06 × 10−3 9.76 × 10−3 8.19 × 10−3
Table 4: Network-width scan at fixed depth and quadrature.
radial–angular quadrature
ε2 [G]
ε2 [Z]
48 × 12 72 × 18 96 × 24
3.31 × 10−3 2.82 × 10−3 3.16 × 10−3
1.36 × 10−2 9.76 × 10−3 8.36 × 10−3
Table 5: Dependence of the neural reconstruction on the quadrature used during training. All residuals are subsequently evaluated on a common finer validation quadrature.
13
G(xm )
Zmax
max αMM
(ε2 [G], ε2 [Z])
5 10 15 20
2.553 2.954 3.036 3.065
1.164 2.001 2.376 2.585
(2.47 × 10−3 , 7.76 × 10−3 ) (2.82 × 10−3 , 9.76 × 10−3 ) (3.21 × 10−3 , 7.47 × 10−3 ) (4.26 × 10−3 , 1.09 × 10−2 )
Table 6: Direct observables and neural errors across the infrared family. The baseline row is included for comparison.
B
Numerical parameters quantity
baseline value
color number subtraction-point coupling infrared ghost condition infrared gluon condition infrared point subtraction point three-gluon damping scale direct external/radial/angular points neural training external/radial/angular points neural validation external/radial/angular points hidden layers and width outer neural Picard steps Adam projection steps per outer iteration residual fine-tuning steps
Nc = 3 α(xs ) = 0.05 G(xm ) = 10 D(xm ) = 15.54 GeV−2 xm = 10−5 GeV2 xs = 7720 GeV2 Λ2s = 1.54 GeV2 140/140/36 80/72/18 120/112/30 2 × 40 70 50 60
Table 7: Baseline physical and numerical parameters.
References [1] William J. Marciano and Heinz Pagels. Quantum chromodynamics: A review. Physics Reports, 36:137–276, 1978. doi: 10.1016/0370-1573(78)90208-9. [2] N. Brambilla et al. Qcd and strongly coupled gauge theories: Challenges and perspectives. European Physical Journal C, 74:2981, 2014. doi: 10.1140/epjc/s10052-014-2981-5. [3] David J. Gross and Frank Wilczek. Ultraviolet behavior of non-abelian gauge theories. Physical Review Letters, 30:1343–1346, 1973. doi: 10.1103/PhysRevLett.30.1343. [4] H. David Politzer. Reliable perturbative results for strong interactions? Letters, 30:1346–1349, 1973. doi: 10.1103/PhysRevLett.30.1346.
Physical Review
[5] Craig D. Roberts and Anthony G. Williams. Dyson–schwinger equations and their application to hadronic physics. Progress in Particle and Nuclear Physics, 33:477–575, 1994. doi: 10. 1016/0146-6410(94)90049-3. [6] Reinhard Alkofer and Lorenz von Smekal. The infrared behavior of qcd green’s functions: Confinement, dynamical symmetry breaking, and hadrons as relativistic bound states. Physics Reports, 353:281–465, 2001. doi: 10.1016/S0370-1573(01)00010-2.
14
[7] Craig D. Roberts and Sebastian M. Schmidt. Dyson–schwinger equations: Density, temperature and continuum strong qcd. Progress in Particle and Nuclear Physics, 45:S1–S103, 2000. doi: 10.1016/S0146-6410(00)90011-5. [8] Christian S. Fischer. Qcd at finite temperature and chemical potential from dyson–schwinger equations. Progress in Particle and Nuclear Physics, 105:1–60, 2019. doi: 10.1016/j.ppnp. 2019.01.002. [9] Adnan Bashir, Lei Chang, Ian C. Cloet, Bruno El-Bennich, Yu-Xin Liu, Craig D. Roberts, and Peter C. Tandy. Collective perspective on advances in dyson–schwinger equation qcd. Communications in Theoretical Physics, 58:79–134, 2012. doi: 10.1088/0253-6102/58/1/16. [10] Gernot Eichmann, Helios Sanchis-Alepuz, Richard Williams, Reinhard Alkofer, and Christian S. Fischer. Baryons as relativistic three-quark bound states. Progress in Particle and Nuclear Physics, 91:1–100, 2016. doi: 10.1016/j.ppnp.2016.07.001. [11] Minghui Ding, Craig D. Roberts, and Sebastian M. Schmidt. Emergence of hadron mass and structure. Particles, 6(1):57–120, 2023. doi: 10.3390/particles6010005. [12] Kenneth G. Wilson. Confinement of quarks. Physical Review D, 10:2445–2459, 1974. doi: 10.1103/PhysRevD.10.2445. [13] S. Aoki et al. Review of lattice results concerning low-energy particle physics. European Physical Journal C, 77:112, 2017. doi: 10.1140/epjc/s10052-016-4509-7. [14] Freeman J. Dyson. The s matrix in quantum electrodynamics. Physical Review, 75:1736–1755, 1949. doi: 10.1103/PhysRev.75.1736. [15] Julian Schwinger. On the green’s functions of quantized fields. i. Proceedings of the National Academy of Sciences, 37:452–455, 1951. doi: 10.1073/pnas.37.7.452. [16] Julian Schwinger. On the green’s functions of quantized fields. ii. Proceedings of the National Academy of Sciences, 37:455–459, 1951. doi: 10.1073/pnas.37.7.455. [17] Markus Q. Huber. Nonperturbative properties of yang–mills theories. Physics Reports, 879: 1–92, 2020. doi: 10.1016/j.physrep.2020.04.004. [18] Markus Q. Huber. Correlation functions of landau gauge yang–mills theory. Physical Review D, 101:114009, 2020. doi: 10.1103/PhysRevD.101.114009. [19] Markus Q. Huber, Christian S. Fischer, and Helios Sanchis-Alepuz. Spectrum of scalar and pseudoscalar glueballs from functional methods. European Physical Journal C, 80:1077, 2020. doi: 10.1140/epjc/s10052-020-08649-6. [20] Markus Q. Huber and Lorenz von Smekal. Spurious divergences in dyson–schwinger equations. Journal of High Energy Physics, 2014(6):015, 2014. doi: 10.1007/JHEP06(2014)015. [21] Christian S. Fischer, Axel Maas, and Jan M. Pawlowski. On the infrared behavior of landau gauge yang–mills theory. Annals of Physics, 324:2408–2437, 2009. doi: 10.1016/j.aop.2009. 07.009. [22] Markus Q. Huber and Lorenz von Smekal. On the influence of three-point functions on the propagators of landau gauge yang–mills theory. Journal of High Energy Physics, 2013(4): 149, 2013. doi: 10.1007/JHEP04(2013)149. [23] Adrian Blum, Markus Q. Huber, Mario Mitter, and Lorenz von Smekal. Gluonic threepoint correlations in pure landau gauge qcd. Physical Review D, 89:061703, 2014. doi: 10.1103/PhysRevD.89.061703. 15
[24] Gernot Eichmann, Richard Williams, Reinhard Alkofer, and Milan Vujinovic. The threegluon vertex in landau gauge. Physical Review D, 89:105014, 2014. doi: 10.1103/PhysRevD. 89.105014. [25] Isaac E. Lagaris, Aristidis Likas, and Dimitrios I. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE Transactions on Neural Networks, 9:987–1000, 1998. doi: 10.1109/72.712178. [26] Maziar Raissi, Paris Perdikaris, and George E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686–707, 2019. doi: 10.1016/j.jcp.2018.10.045. [27] George E. Karniadakis, Ioannis G. Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3:422–440, 2021. doi: 10.1038/s42254-021-00314-5. [28] David Pfau, James S. Spencer, Alexander G. D. G. Matthews, and W. M. C. Foulkes. Ab initio solution of the many-electron schrödinger equation with deep neural networks. Physical Review Research, 2:033429, 2020. doi: 10.1103/PhysRevResearch.2.033429. [29] Lorenzo Brevi, Antonio Mandarino, and Enrico Prati. Addressing the non-perturbative regime of the quantum anharmonic oscillator by physics-informed neural networks. New Journal of Physics, 26:103015, 2024. doi: 10.1088/1367-2630/ad7f4b. [30] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George E. Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3:218–229, 2021. doi: 10.1038/s42256-021-00302-5. [31] Rodrigo Carmo Terin. Physics-informed neural networks viewpoint for solving the dyson–schwinger equations of qed. SciPost Physics Core, 8:054, 2025. doi: 10.21468/ SciPostPhysCore.8.3.054. [32] Rodrigo Carmo Terin. Spectral functions in minkowski quantum electrodynamics from neural reconstruction. Journal of High Energy Physics, 2026. in press. [33] Ludvig D. Faddeev and Victor N. Popov. Feynman diagrams for the yang–mills field. Physics Letters B, 25:29–30, 1967. doi: 10.1016/0370-2693(67)90067-6. [34] Lorenz von Smekal, Andreas Hauck, and Reinhard Alkofer. A solution to coupled dyson– schwinger equations for gluons and ghosts in landau gauge. Annals of Physics, 267:1–60, 1998. doi: 10.1006/aphy.1998.5806. [35] Lorenz von Smekal, Reinhard Alkofer, and Andreas Hauck. The infrared behavior of gluon and ghost propagators in landau gauge qcd. Physical Review Letters, 79:3591–3594, 1997. doi: 10.1103/PhysRevLett.79.3591. [36] Lorenz von Smekal, Kim Maltman, and Andre Sternbeck. The strong coupling and its running to four loops in a minimal mom scheme. Physics Letters B, 681:336–342, 2009. doi: 10.1016/ j.physletb.2009.10.030. [37] Philippe Boucaud, Feliciano De Soto, Jean-Pierre Leroy, Alain Le Yaouanc, Jean Micheli, Olivier Pene, and Jose Rodriguez-Quintero. Ghost-gluon running coupling, power corrections and the determination of λms . Physical Review D, 79:014508, 2009. doi: 10.1103/PhysRevD. 79.014508. [38] Adam Paszke et al. Pytorch: An imperative style, high-performance deep learning library. In Advances in Neural Information Processing Systems, volume 32, 2019. 16
[39] Konrad Osterwalder and Robert Schrader. Axioms for euclidean green’s functions. Communications in Mathematical Physics, 31:83–112, 1973. doi: 10.1007/BF01645738. [40] Konrad Osterwalder and Robert Schrader. Axioms for euclidean green’s functions. ii. Communications in Mathematical Physics, 42:281–305, 1975. doi: 10.1007/BF01608978.
17