Light-quark mass dependence of the nucleon axial charge and pion-nucleon scattering phenomenology

The light-quark mass dependence of the nucleon axial isovector charge ($g_A$) has been studied up to next-to-next-to-leading order, $O(p^4)$, in relativistic chiral perturbation theory using extended-on-mass-shell renormalization, without and with explicit $\Delta (1232)$ degrees of freedom. We show that in the $\Delta$-less case, at this order, the flat trend of $g_A (M_{\pi})$ exhibited by state-of-the-art lattice QCD (LQCD) results cannot be reproduced using low energy constants (LECs) extracted from pion-nucleon elastic and inelastic scattering. A satisfactory description of these LQCD data is only achieved in the theory with $\Delta$. From this fit we report $g_A (M_{\pi\rm{(phys)}}) = 1.260 \pm 0.012$, close to the experimental result, and $d_{16}= -0.88\pm 0.88$ GeV$^{-2}$, in agreement with its empirical value. The large uncertainties are of theoretical origin, reflecting the difference between $O(p^3)$ and $O(p^4)$ that still persists at large $M_{\pi}$ in presence of the $\Delta$.


INTRODUCTION
The axial isovector charge g u−d A (g A from now on) is a fundamental property of the nucleon related to the difference in the spin fractions carried by u and d quarks. With a magnitude precisely determined from neutron β decay 1 , the nucleon axial charge represents a benchmark for non-perturbative studies of Quantum Chromodynamics in the lattice (LQCD) alongside with other nucleon properties, such as scalar and tensor charges, electromagnetic form factors and parton distribution functions (see Ref. [3] and Sec. 10 of Ref. [4] for recent reviews). Over the last decade, progress in this direction has been significant [4], owing to better computational resources, improved algorithms and techniques to reduce systematic errors, in particular those induced by excited-state contamination, which can be considerable in the baryon sector for currently available lattice ensembles [5,6].
As the effective theory of QCD in the non-perturbative regime, chiral perturbation theory (ChPT) plays an important role. It allows to determine the light (u, d) quark mass dependence of low-energy hadronic observables in terms of low-energy constants (LECs) and perform model-independent extrapolations of LQCD results to the physical point. ChPT also allows to account for finite lattice-volume and lattice spacing corrections in a systematic way [7,8]. It has also proved helpful to deal with the contamination from excited states [5,6,9]. The interplay between ChPT and LQCD can also be used to determine poorly known LECs, which are difficult to access with experimental data. This is the case of d 16 present in the O(p 3 ) part of the πN Lagrangian. Via a 4d 16 M 2 π term, this LEC controls the slope of the light-quark mass dependence of g A . 2 Therefore, its extraction from LQCD results at low pion masses is only natural. LEC d 16 has been identified as one of the most significant sources of uncertainty in quark mass dependence of nuclear properties such as ground-state and binding energies through long-range nuclear forces [10][11][12].
Armed with relativistic baryon ChPT, we revisit this problem up to next-to-next-to-leading order (NNLO) in the perturbative expansion. The EOMS scheme ensures that not only power counting but also analytic properties of loop functions are properly preserved. After describing in Sec. II the details of the ChPT calculation and introducing the relevant terms of the effective Lagrangian we show in Sec. III the M π dependence of g A predicted with the phenomenological LECs obtained in Ref. [26] from πN elastic and inelastic scattering. In order to obtain a better description of recent LQCD data we have performed fits from which some of the LECs and, in particular, d 16 , are determined. These are presented in Sec. IV. A meaningful representation of the flat trend exhibited by the LQCD results is only achieved with explicit ∆(1232). Differences between orders, which are considerable between NLO and NNLO, are used to provide a measure of the systematic error arising from the truncation of the perturbative series. These differences are the main source of uncertainty in the extracted d 16 value.

II. THE NUCLEON AXIAL CHARGE IN RELATIVISTIC BCHPT
In the isospin limit, the matrix element of the axial isovector current with q = (u, d) T the quark-field doublet, taken between on-shell nucleon states of equal four-momenta p (p 2 = m 2 N ), can be written as The isovector character of the current is manifest given the presence of the Pauli isospin matrices τ a . We calculate the axial charge g A employing relativistic BChPT with pions, nucleons and ∆ as degrees of freedom. We use the standard power counting [36], extended to diagrams with ∆ baryons following small scale expansion (SSE) of Refs. [37,38]. In this power counting, the mass difference δ = m N − m ∆ is O(p).

A. Relevant terms of the effective Lagrangian
In this section the terms in the effective Lagrangian, L eff , required for our calculation are presented. We need where superscripts indicate the chiral order and subscripts show the present degrees of freedom. The terms in L (1,3) πN , L (1) π∆ and L (1) πN ∆ that are relevant for our study have been displayed in Ref. [34] using the same notation adopted here. In addition, O(p 2 ) contributions are entailed to give g A at NNLO. Following the notation of Ref. [39], the required terms of the πN Lagrangian are where Ψ is the isospin doublet of nucleon fields. Including only isovector axial external fields a a µ τ a /2 ≡ a µ , we end up with , and ... denotes a trace over isospin. In L (2) πN ∆ , after redundant terms are eliminated from Eq. (67) of Ref. [40] (see also Sec. 3.1 of Ref. [41] and the Appendix of Ref. [42]) only the following monomials contribute to the nucleon axial charge. Here Ψ µ denotes the Rarita-Schwinger field of the ∆ resonance. States ij = δ ij − τ i τ j /3 are isospin-3/2 projectors, are also isospin doublets, whose explicit expressions in terms of the physical ∆ states are derived, for example, in Appendix A of Ref. [43]. The isovector projection of u µ is represented by ω µ,i = 1 2 τ i u µ (see for instance Ref. [40]). Finally, also at O(p 2 ) [40,44], introduces an M π dependent correction to the ∆ mass, in the same way as the term proportional to c 1 does for the nucleon mass.

B. Perturbative calculation
The set of Feynman diagrams that contribute to the axial form factor and, in particular, to g A are displayed in Fig. 1  and 2. The LECs that each of these diagrams brings about are listed in Table I. At O(p 4 ), there are contributions from O(p 2 ) vertices but also baryon (N, ∆) mass insertions. The later are calculated perturbatively, i.e. we evaluate directly diagrams (i), (j), (l)-(n) of Fig. 2, avoiding Dyson re-summations to all orders in the propagators. Alternative choices have been considered in Refs. [33,35].  Although not diagrammatically represented, nucleon wave-function renormalization is taken into account in the standard way. To the order of our calculation, only the O(p) diagram (a), giving just the axial charge in the chiral limitg A , should be multiplied by the wave-function renormalization constant Z N , calculated from the O(p 4 ) nucleon self-energy.
Additional O(p 3 ) and O(p 4 ) contributions to g A are generated in this way. The LECs that enter in them are also listed in Table I. In order to absorb the ultraviolet divergencies generated by loops, we perform dimensional renormalization in the MS scheme [45], at the scale µ =m, wherem denotes the nucleon mass in the chiral limit. An additional renormalization is then performed to account for the power-counting breaking caused by the presence of baryon masses which do not vanish in the chiral limit. Among the available schemes we employ EOMS renormalization [45] which consists in the absorption of the power counting breaking (PCB) terms in a redefinition of the LECs. In this way, power counting is restored without altering the analytic properties of the loops and preserving covariance. Note that, as in earlier studies [24,34], PCB terms are identified and subtracted in an expansion in powers or M π but not in δ. The EOMS renormalization shifts for the LECs are lengthy, so we have included them in the supplementary material.
In the way outlined above we obtain the axial charge within the EOMS renormalization scheme up to O(p 4 ) with explicit ∆. Our result has the following structure, with superindices indicating the chiral order: The structures arising from the loops are preserved at the price of keeping higher order terms. The O(p 3 ) part of g A is given in Eq. (A4) of Ref. [34].  [34]. At O(p 4 ) we agree with the / ∆ EOMS expressions given in Ref. [22,33] except for the different treatment of nucleon mass insertions, which in our case, do not include resummations. As a result, c 1 M 2 π factors appear only linearly in our calculation but not at all orders. Furthermore, a PCB term proportional to c 4 present in Ref. [22] after the EOMS renormalization is absent in our final result.

III. LIGHT-QUARK MASS DEPENDENCE OF THE AXIAL COUPLING FROM PHENOMENOLOGICAL INPUT
Once the framework and the chiral order are established, the M π dependence of g A is ultimately determined by the LEC values. In Ref. [26], elastic πN and inelastic πN → ππN scattering has been studied up to O(p 4 ) in covariant ChPT using a modified version of the EOMS approach [25]. The LECs that enter g A (M π ) at O(p 4 ) in the ∆-less model were extracted, including d 16 , owing to the inclusion of low energy total and double differential πN → ππN cross section data in the combined analysis. To make a prediction of g A (M π ) based on this phenomenological input, we should transform the LECs from the modified EOMS of Refs. [25,26] to the conventional one adopted here. To the order we are working at, this transformation, whose details are disclosed in Appendix A, affects the numerical value of d 16 but not of c 1−4 . The axial coupling in the chiral limit,g A , which is not extracted in Refs. [25,26], is determined as from the experimental value, precisely known from β decay [1]. Up to higher orders for g A , the remaining parameters, i.e. the pion decay constant and the nucleon mass in the chiral limit, have been fixed to F π F π(phys) = 92.2 MeV andm m N (phys) + 4c 1 M 2 π(phys) , with M π(phys) = 135 MeV and m N (phys) = 939 MeV, respectively. The numerical values for d 16 , c 1−4 as well asg A and their uncertainties are summarized in Table II.   TABLE II: LECs used to obtain the g A (M π ) dependence in the / ∆ case.
The resulting g A (M π ) at O(p 3 ) and O(p 4 ) are displayed in Fig. 3 together with a subset of recent LQCD determinations 4 . The curves are accompanied by 1σ statistical error bands arising from the uncertainties in the LECs of Table II [26] and 1σ error bands. The LQCD data points are from CalLat 18 [46] (black circles), Mainz 19 [47] (red crosses), RQCD 19 [9] (green triangles) and NME 21 [48] (blue squares).
increasing tension as M π grows and a curvature in the central value that is absent in the data. On the other hand, it is immediately apparent that the O(p 4 ) prediction not only does not improve the O(p 3 ) but plainly fails to describe the M π dependence of LQCD data. We have checked that alternative c 1−4 determinations from earlier πN analyses using EOMS [22,24] do not mitigate the steep rise of g A (M π ) at O(p 4 ). Therefore, the inability to reconcile the light-quark mass dependence ofg A at O(p 4 ) with phenomenology, earlier observed in non-relativistic heavy-baryon ChPT [29,30], is also a feature of the relativistic version of the theory.
The problem with the O(p 4 ) result might be caused by an accidental slow convergence of g A in BChPT. Instead, one could argue that additional degrees of freedom such as the ∆ [30] or even the Roper N (1440) resonance might solve the issue. In Ref. [26], a version of the theory incorporating the ∆ pole is also considered. However, in order to predict g A (M π ) up to O(p 4 ) with explicit ∆, additional LECs absent in that study are needed. For this reason we do not tackle the role of the ∆ in the g A (M π ) dependence from purely empirical input but postpone it to the next Sec. IV where our fits to recent LQCD data are reported.

A. Data set and fit strategy
The axial charge has been historically underpredicted in the lattice (see for instance Fig. 2 of Ref. [49]), but, as a result of conceptual and technical improvements, the majority of LQCD results nowadays agree with the experimental value at the level of a few percent [4]. Most significantly, the analysis of excited-state contamination has notably evolved in last years. Therefore, we only include in our data set results from studies with an improved treatment of these effects. We take into account renormalized where a stands for the lattice spacing, from 5 CalLat 18 [46], Mainz 19 [47], RQCD 19 [9] 6 and NME 21 [48] 7 . Our analysis treats 2+1+1 (CalLat 18) and 2+1 (Mainz 19, RQCD 19, NME 21) ensembles on the same footing, assuming that the c-quark sea content plays a negligible role. We disregard LQCD determinations of g A coming from q 2 , finite volume, a or M π extrapolations. Since we do not correct the finite volume effects, we consider only large volumes, satisfying M π L ≥ 3.5, so that the neglected extrapolation is small and can be absorbed in the errors.
In order to gauge the ability of BChPT at O(p 3 ) and O(p 4 ) to describe the M π dependence of g A , and to extract g A and d 16 LECs, we perform fits to the LQCD data set defined above. With this aim we define the following χ 2 : In addition to g A (M π ) from BChPT, our theoretical parametrization incorporates lattice discretization corrections as Free parameters x j , with j = 1, 2, 3, 4 when lattice data point i ∈ {CalLat 18, Mainz 19, RQCD 19, NME 21} respectively, control the leading a-dependence of the LQCD data, which is action specific: n 1,4 = 1, while n 2,3 = 2. These discretization corrections are small and do not substantially change the extracted LECs, but result in an appreciable reduction of the fits' χ 2 /dof. Some of the LECs upon which g A (M π ) depends are well determined in other studies and are kept fixed, while others are treated as free parameters together withg A and d 16 . To improve our description of the LQCD data and reduce correlations [50], for free LECs we assume naturalness Λ n−1 c n ∼ 1, encoded in Gaussian priors Here, c n generically denotes a LEC of chiral order O(p n ); the breakdown scale is set to Λ = 1 GeV∼ 4πF π [51,52]. We anticipate that a prior forg A is superfluous, since its value is always driven to a natural one by low M π g i A LQCD data.
The large relative contribution of O(p 4 ) terms discussed in the previous section and illustrated by Fig. 3 is a clear indication that the theoretical error associated with truncation of the perturbative expansion should be taken into account. We follow the method proposed in Ref. [53]. Let X be an observable with a chiral perturbative expansion ∆X (m) = X (m) − X (m−1) encompasses all the monomials that start at order m. If X is calculated up to order n, X ≈ X (n) , assuming that the truncation error is dominated by order n + 1, its contribution ∆X (n+1) can be estimated in a conservative way as [25,53] where Q denotes the expansion variable.
Connecting to Eq. (9), ∆g A(loop) and ∆g (4) A(loop) . Based on the results in Sec. III it is easy to foresee that at larger M π , ∆g (5) Aχ will be determined by the last term in Eq. (16). In our O(p 3 ) fits we do not assume any prior knowledge about O(p 4 ) and, therefore, the truncation error is given by ∆g (4) In our χ 2 , Eq. (11), this theoretical error is added in quadrature to the one of LQCD points: This recipe assigns larger uncertainties to points at high M π , where the convergence of the chiral expansion is poorer, therefore reducing their impact on the fit. We should also mention that, as discussed in Refs. [25,53,54], although this theoretical error estimate does not have a clear statistical interpretation, results are similar to those obtained in a preferable but beyond the scope of the present study Bayesian approach [50]. Fits are carried out iteratively until convergence is achieved. The first minimization is performed without theoretical errors, which are subsequently evaluated using the LECs determined in the previous step. The lattice discretization parameters x j are fixed in the first iteration because evolving them results in overfitting. With the described procedure, LQCD and truncation errors are not independent and it is not obvious how to combine them in the error for a given observable. Consequently, following Ref. [25], we plot them separately in the error bands for g A (M π ). Furthermore, as a quantitative measure of the agreement of our best-fit curve with the LQCD data, we also provide the χ 2 0 value, defined as Finally, we have tested the convergence range of our calculation by varying the maximum M π of the lattice data included in the fits in the range of M πcut ∈ [200, 402] MeV. A plateau in the χ 2 and the fitted LECs is found towards the end of the interval. In consequence we report results including all points in the chosen data sets up to M π = 402 MeV 8 . As discussed in Sec. IV C 2, the theoretical error becomes larger at higher M π where the convergence is poorer and, therefore, the corresponding LQCD points weight less in the fits.  Table III. At the first sight, this model describes LQCD data well, with a good χ 2 /dof, relatively small errors and naturalg A and d 16 . As apparent from the comparison of the first columns of Tables II and III, the d 16 value is slightly above the phenomenological one extracted in Ref. [26] and with much smaller error. Instead,g A is in tension with the value obtained from experimental input, Eq. (10). As a consequence g A (M π(phys) ) = 1.205 ± 0.010, well below the experimental result. However, an inspection of the results for the O(p 4 ) model described below reveals a large contribution from the O(p 4 ) terms at M π 200 MeV, in line with the findings of Sec. III. The fact that these terms are considerably larger than the error band of the O(p 3 ) result implies that, in this case, the theoretical error estimated from O(p 1 ) and O(p 3 ) terms falls short in accounting for the O(p 4 ) contribution. The agreement of g A (M π ) at O(p 3 ) with LQCD points appears then as misleadingly good, while the uncertainties in the LEC values can be regarded as underestimated. In other words, to be realistic, an O(p 3 ) analysis should be limited to small M π < 200 MeV, which would be unfeasible with the small amount of presently available LQCD data in this region.

O(p 4 )
In the O(p 4 ) / ∆ model one has the additional contribution of NLO LECs c 1−4 . They were initially allowed to evolve in the fits under the constrains set by their empirical determination [26] (second column of Table II) implemented as Gaussian priors. However, as these LECs are quite well determined, this procedure yields substantially the figures favored by the priors. Therefore, we report here the results of a simpler fit with c 1−4 held fixed to their central phenomenological values (second column of Table III).
In any case, as apparent from the lower left panel of Fig. 4 the O(p 4 ) / ∆ model fails to describe the light-quark mass dependence of g A . 9 The small χ 2 is merely a consequence of the large theoretical error, which reduces the impact of high M π points on the fit. Nevertheless, the poor agreement is reflected in the magnitude of χ 2 0 and also in the quite unnatural d 16 in spite of its prior. This large (in absolute value) d 16 is inconsistent with its phenomenological value and nonetheless unable to correct the M π dependence at M π 300 MeV, which is largely dominated by the O(p 4 ) terms and, therefore, very similar to the one displayed in Fig. 3.
The fact that the very wide theoretical error band encompasses the lattice points reflects that the disagreement would be removed by a (large) contribution of O(p 5 ) counterbalancing the O(p 4 ) ones. Actually, within heavy-baryon ChPT it has been shown [29] that the curve can be bent down by additional contributions of orders O(p 5,6 ) with LECs of natural size (see also Fig. 1 of Ref. [11]). Here we take a different avenue and introduce the ∆(1232) explicitly as advocated in Ref. [30] based on the Adler-Weisberger sum rule [55,56] and a heavy-baryon ChPT calculation for g A .  (4,5) Aχ . The LQCD points from CalLat 18 [46] (black circles), Mainz 19 [47] (red crosses), RQCD 19 [9] (green triangles) and NME 21 [48] (blue squares) are shown at their finite a values i.e. without (small) discretization corrections.   coupling h A to its large-N c value h A = 1.35 [26], which is close to its empirical value [24]. For the L (1) π∆ coupling g 1 , whose impact on g A is small, the large-N c limit yields |g 1 | = 2.29 [24,26]. We are unable to discern its sign in our fits to g A LQCD data and opt for g 1 = −2.29, preferred both by πN elastic scattering [24] and our own studies of the nucleon axial form factor at low q 2 to be reported elsewhere. Finally, in analogy to the nucleon case, m ∆ m ∆(phys) − 4 a 1 M 2 (phys) with m ∆(phys) = 1232 MeV and a 1 = 0.90 GeV −1 from the LQCD M π dependence of the ∆ mass [57].

O(p 3 )
With this model, the result of the fit for g A (M π ) closely resembles the O(p 3 ) / ∆ one as can be seen in the comparison of the upper panels of Fig. 3. However, the value of d 16 changes considerably, including its sign, with respect to the O(p 3 ) / ∆ one. The same feature was obtained when this LEC was extracted from LQCD results for the axial form factor at low q 2 with O(p 3 ) relativistic BChPT [34]. The O(p 4 ) fit described in the next section produces O(p 4 ) terms larger than the theoretical error estimated from O(p 1 ) and O(p 3 ) ones (although to a lesser extent than in the / ∆ case) making the extension to O(p 4 ) desirable for a more realistic determination of LECs and their uncertainties.

O(p 4 )
In the fit of the O(p 4 ) model with explicit ∆ we fix the c 1−4 LECs to the values extracted from the πN scattering [25] (see the last column of Table III). They take the ∆ pole into account and are in good agreement with the joint πN +ππN fit from Ref. [26]. In addition tog A and d 16  The result of the fit, depicted in the lower right panel of Fig. 4, satisfactorily describes the trend of g A (M π ) as predicted by LQCD up to relatively large M π . The theoretical error is large and rapidly increasing with M π due to the large contribution of O(p 4 ) terms. Nevertheless, the description of LQCD data and convergence are notably improved with respect to the O(p 4 ) / ∆ model. The extracted LECs are given in Table III. The b 4 value might seem unnatural but one should keep in mind that it is a combination of LECs. As shown in Table IV the correlations among LECs are sizable; they are an indication of degeneracy among the parameters that could be partially lifted by adding a new dimension to the fit (i.e. studying the axial form factor at finite q 2 ). The d 16 value obtained from this model, d 16 = −0.88 ± 0.88 GeV −2 , is in good agreement with the determination from inelastic πN → ππN with explicit ∆ pole [26] which, translated to standard EOMS, is d 16(pheno) = −1.0 ± 1.0 GeV −2 . Although, as argued in the Introduction, the M π dependence of g A is in principle better suited to extract d 16 than the available experimental πN → ππN data, the convergence issues of the former lead to large errors, comparable to those of the phenomenological result. Theg A value is higher than in the O(p 3 ) fits and with a larger error. Furthermore, at the physical point g A (M π(phys) ) = 1.260 ± 0.012 is now close to the experimental value although with a much larger uncertainty. The stability of the results forg A and d 16 as a function of the maximum M π considered, M πcut , is shown in Fig. 5. One can see that for both quantities the numerical values and their errors stabilize for M π 300 MeV. Owing to the consideration of theoretical errors, the analysis could be extended to a broader range of M π because points at higher M π , where the convergence is worse, have a larger combined uncertainty, Eq. (18), and contribute less to the fit.

V. CONCLUSIONS AND OUTLOOK
We have investigated the pion mass dependence of the nucleon axial coupling up to O(p 4 ) ≡ O(M 3 π ) (NNLO) in relativistic BChPT with EOMS renormalization. We have shown that, at this order and without including ∆(1232) explicitly, the pion mass dependence of g A obtained using LECs extracted from phenomenological analyses of pionnucleon elastic and inelastic (πN → ππN ) scattering cannot describe the rather flat behavior predicted by stateof-the-art LQCD simulations. The disagreement is manifest from pion masses right above the physical point. This feature, earlier observed in non-relativistic heavy-baryon ChPT [29,30], is therefore also present in the relativistic theory. The fact that O(p 4 ) terms become large from relatively low M π 200 MeV implies that O(p 3 ) analyses of g A (M π ) likely underestimate theoretical uncertainties.
In line with the findings of Ref. [30] within heavy-baryon ChPT, we can satisfactorily describe the LQCD data for g A (M π ) at O(p 4 ) only after the ∆ is incorporated as an explicit degree of freedom. However, although in a much smaller degree than in the / ∆ case, a fast increase in the relative size of O(p 4 ) terms with M π is observed and reflected by the estimate of theoretical uncertainties. This fact jeopardizes the precision of the ChPT description of g A (M π ) at M π 300 MeV and negatively influences the extraction of LECs based on LQCD data at such M π . Together with the sizable correlations, which can be reduced by considering both the M π and q 2 dependence of the axial form factor in the fits, this feature implies that heavier resonances and/or O(p 5 ) terms are still required to reach a good convergence and reduce theoretical errors. The impact of setting the baryon masses in the loops to the values obtained by the LQCD simulations may be also worth exploring in view of the results of Ref. [35] although this would correspond to the resummation of baryon selfenergy insertions of higher order. For this purpose it would be relevant to have more information about the ∆ pole position for the different lattice ensembles.
From our O(p 4 ) analysis of recent LQCD data we have obtained g A (M π(phys) ) = 1.260 ± 0.012, close to the experimental value, and d 16 = −0.88 ± 0.88 GeV −2 , in agreement with πN phenomenology. As a consequence of the previously discussed issues, errors are still large, particularly for d 16 that is naturally extracted from the light-quark mass dependence of the nucleon axial coupling. New LQCD results at M π 250 MeV will improve the precision of the extracted LECs. Besides, our ongoing effort to extend the analysis to the whole axial form factor (at low q 2 ) may shed more light on d 16 , as well as other LECs such as d 22 .

ACKNOWLEDGMENTS
We thank G. Guerrero, S. Leupold and M. J. Vicente Vacas for useful discussions, and D. Yao for sharing details of his related calculation with us. E. Epelbaum, A. Gasparyan and H. Krebs kindly answered our questions about their renormalization prescription. We are also indebted to the RQCD team for making their axial form-factor results available to us and to R. Gupta for providing relevant information about the NME ones. This research has been partially supported by the Spanish Ministerio de Ciencia e Innovación under contracts FIS2017-84038-C2-1-P and PID2020-112777GB-I00, the EU STRONG-2020 project under the program H2020-INFRAIA-2018-1, grant agreement no. 824093 and by Generalitat Valenciana under contract PROMETEO/2020/023.

Appendix A: LEC conversion
The covariant renormalization prescription adopted in Refs. [25,26] differs from EOMS. In EOMS only powercounting breaking terms are absorbed in the LECs. To obtain an equivalence between covariant and heavy-baryon results, the prescription of Ref. [25] also subtracts infrared-regular terms at the order of the calculation, as well as pieces proportional to log M 2 π /m 2 N . After settingλ = 0 and µ = m N in Eq. (21) of Ref. [25] one is left with a dimensionally regularised LEC, x r : which allows to express the renormalized EOMS LECs x EOMS in terms of the corresponding ones,x cov , from Refs. [25,26]. The β-functions β x required to cancel the mesonic tadpole terms are listed in Appendix B of Ref. [25] as well as finite shifts δx once β c1−4 = 0 [25] and the finite shifts coincide with those in EOMS. Equations A2-A4 are derived for the / ∆ case but hold also for the model with the ∆ pole of Ref. [25] as it does not involve additional renormalizations. To the order of the present study Z N , defined in Eq. 8, can be written as Z N = 1 + 1 16π 2 F 2 π δZ (3)loop N + δZ (4)loop N . The contribution to g A that starts at O(p 4 ) is theng A δZ (4)loop N , with δZ (4)loop N , including ∆, given by δZ (4)