Bjorken polarized sum rule and infrared-safe QCD couplings

TL;DR

This paper compares different QCD models, including perturbative and analytic versions, in fitting polarized Bjorken sum rule data, highlighting the advantages of infrared-safe couplings without Landau singularities.

Contribution

It introduces a comprehensive comparison of four QCD variants, demonstrating improved fits with analytic QCD models that avoid Landau singularities at low Q^2.

Findings

2δ and 3δ analytic QCD models provide the best fit results

Infrared-safe QCD couplings improve the fitting procedure

Fitting quality decreases when elastic contributions are included in total BSR

Abstract

Experimental data obtained for the polarized Bjorken sum rule (BSR) $\Gamma_1^{p-n}(Q^2)$ are fitted by using predictions derived within a truncated operator product expansion (OPE) approach to QCD. Four QCD versions are considered: perturbative QCD (pQCD) in the MSbar scheme, Analytic Perturbation Theory (APT), and 2$\delta$ and 3$\delta$ analytic QCD versions. In contrast to pQCD, these QCD variants do not have Landau singularities at low positive $Q^2$, which facilitates the fitting procedure significantly. The fitting procedure is applied first to the experimental data of the inelastic part of BSR, and the known elastic contributions are added after the fitting. In general, when 2$\delta$ and 3$\delta$ QCD coupling is used the fitted curves give the best results, within the $Q^2$-range of the fit as well as in extended $Q^2$-intervals. When the fitting procedure is applied to the total BSR, i.e., to the sum of the experimental data and the elastic contribution, the quality of the results deteriorates significantly.

Tables

Table 1. Table 1: The values of the extracted fit parameters k = μ 2 / Q 2 𝑘 superscript 𝜇 2 superscript 𝑄 2 k=\mu^{2}/Q^{2} and f 2 p − n ​ ( Q in 2 ) superscript subscript 𝑓 2 𝑝 𝑛 subscript superscript 𝑄 2 in f_{2}^{p-n}(Q^{2}_{\rm in}) ( Q in 2 = 1 ​ GeV 2 subscript superscript 𝑄 2 in 1 superscript GeV 2 Q^{2}_{\rm in}=1\ {\rm GeV}^{2} ), obtained by fitting the OPE expression ( 10 ), truncated at D = 2 𝐷 2 D=2 ( μ 6 = 0 subscript 𝜇 6 0 \mu_{6}=0 ), for the considered four QCD variants, and for the minimal Q 2 superscript 𝑄 2 Q^{2} values of the fit being 0.66 0.66 0.66 , 0.47 0.47 0.47 and 0.268 ​ GeV 2 0.268 superscript GeV 2 0.268\ {\rm GeV}^{2} . The uncertainties of the extracted f 2 p − n superscript subscript 𝑓 2 𝑝 𝑛 f_{2}^{p-n} values are separated to those coming from the statistical (small) and systematic (large) uncertainties of the experimental data (cf. Appendix D ). The resulting fit quality parameters χ 2 / d . o . f . formulae-sequence superscript 𝜒 2 d o f \chi^{2}/{\rm d.o.f.} are given as well (see the text for explanation).

QCD variant	$Q_{\min }^2\left(\mathrm{fit}\right)$	$k$	$f_2^{p-n}(1.)$	${\chi }^2/d.o.f$	${\chi }_{\mathrm{ext}}^2/d.o.f$	${\chi }_{0.268}^2/d.o.f.$	${\chi }_{\mathrm{all}}^2/d.o.f.$
$\overline{\mathrm{MS}}$ pQCD	0.66	6.04	$-0.103\pm 0.001\pm 0.011$	44.5	161.	$6.81\times {10}^4$	$\infty$
(F)APT	0.66	2.74	$-0.143\pm 0.001\pm 0.021$	23.3	85.6	63.0	970.
2$\delta$	0.66	0.397	$-0.077\pm 0.001\pm 0.024$	8.34	52.2	11.1	406.
3$\delta$	0.66	2.41	$-0.064\pm 0.001\pm 0.023$	6.85	56.1	6.17	130.
$\overline{\mathrm{MS}}$ pQCD	0.47	6.35	$-0.094\pm 0.000\pm 0.009$	66.1	177.	$5.71\times {10}^4$	$\infty$
(F)APT	0.47	3.03	$-0.137\pm 0.000\pm 0.015$	30.5	88.4	51.2	873.
2$\delta$	0.47	0.395	$-0.077\pm 0.000\pm 0.023$	7.91	45.2	10.9	402.
3$\delta$	0.47	2.58	$-0.064\pm 0.000\pm 0.023$	6.49	50.2	5.92	129.
$\overline{\mathrm{MS}}$ pQCD	0.268	16.	$-0.160\pm 0.001\pm 0.045$	58.7	145.	58.7	$1.20\times {10}^3$
(F)APT	0.268	3.71	$-0.134\pm 0.000\pm 0.102$	49.6	103.	49.6	827.
2$\delta$	0.268	0.338	$-0.074\pm 0.000\pm 0.022$	10.6	46.6	10.6	372.
3$\delta$	0.268	2.55	$-0.064\pm 0.000\pm 0.018$	5.91	44.5	5.91	129.

Table 2. Table 2: The values of the extracted fit parameters k = μ 2 / Q 2 𝑘 superscript 𝜇 2 superscript 𝑄 2 k=\mu^{2}/Q^{2} , f 2 p − n ​ ( 1 ​ GeV 2 ) superscript subscript 𝑓 2 𝑝 𝑛 1 superscript GeV 2 f_{2}^{p-n}(1\ {\rm GeV}^{2}) and μ 6 subscript 𝜇 6 \mu_{6} (in GeV 4 superscript GeV 4 {\rm GeV}^{4} ), obtained by fitting the OPE expression ( 10 ), truncated at D = 4 𝐷 4 D=4 ( μ 6 ≠ 0 subscript 𝜇 6 0 \mu_{6}\not=0 ). The notations are as in Table 1 . See the text in Sec.  III.1 for explanation of the various χ 2 / d . o . f . formulae-sequence superscript 𝜒 2 d o f \chi^{2}/{\rm d.o.f.} ’s.

QCD variant	$Q_{\min }^2\left(\mathrm{fit}\right)$	$k$	$f_2^{p-n}(1.)$	${\mu }_6$	${\chi }^2/d.o.f$	${\chi }_{\mathrm{ext}}^2/d.o.f$	${\chi }_{0.268}^2/d.o.f.$	${\chi }_{\mathrm{all}}^2/d.o.f.$
$\overline{\mathrm{MS}}$ pQCD	0.66	16.	$-0.219\pm 0.002\pm 0.111$	$0.032\pm 0.001\pm 0.055$	7.31	53.4	$2.92\times {10}^5$	$1.68\times {10}^6$
(F)APT	0.66	0.063	$-0.198\pm 0.002\pm 0.106$	$0.018\pm 0.001\pm 0.028$	7.85	44.2	17.7	$1.04\times {10}^4$
2$\delta$	0.66	0.999	$-0.116\pm 0.002\pm 0.094$	$0.014\pm 0.001\pm 0.023$	6.99	51.8	19.5	$5.92\times {10}^3$
3$\delta$	0.66	4.46	$-0.101\pm 0.002\pm 0.114$	$0.013\pm 0.001\pm 0.025$	6.45	57.5	28.0	$7.27\times {10}^3$
$\overline{\mathrm{MS}}$ pQCD	0.47	16.	$-0.216\pm 0.002\pm 0.120$	$0.031\pm 0.001\pm 0.045$	7.34	47.3	$2.85\times {10}^5$	$1.66\times {10}^6$
(F)APT	0.47	0.063	$-0.194\pm 0.002\pm 0.098$	$0.016\pm 0.001\pm 0.021$	7.65	35.6	10.0	$6.95\times {10}^3$
2$\delta$	0.47	0.772	$-0.098\pm 0.002\pm 0.078$	$0.008\pm 0.001\pm 0.017$	6.88	48.2	7.42	$1.40\times {10}^3$
3$\delta$	0.47	3.10	$-0.072\pm 0.002\pm 0.138$	$0.003\pm 0.001\pm 0.026$	6.64	53.5	6.56	218.
$\overline{\mathrm{MS}}$ pQCD	0.268	3.00	$-0.180\pm 0.001\pm 0.237$	$0.023\pm 0.000\pm 0.013$	6.52	38.7	6.52	$\infty$
(F)APT	0.268	5.27	$-0.181\pm 0.001\pm 0.194$	$0.013\pm 0.000\pm 0.011$	7.77	37.0	7.77	$3.91\times {10}^3$
2$\delta$	0.268	0.149	$-0.089\pm 0.001\pm 0.166$	$0.005\pm 0.000\pm 0.010$	6.27	41.9	6.27	293.
3$\delta$	0.268	2.72	$-0.007\pm 0.001\pm 0.262$	$0.001\pm 0.000\pm 0.015$	5.99	46.5	5.99	47.3

Table 3. Table 3: As in Table 2 , but now the OPE form is “massive”, Eq. ( 16 ). The fits were made only for Q min 2 ​ ( fit ) = 0.268 ​ GeV 2 subscript superscript 𝑄 2 min fit 0.268 superscript GeV 2 Q^{2}_{\rm min}({\rm fit})=0.268\ {\rm GeV}^{2} . The extracted values of the squared mass M 2 superscript 𝑀 2 M^{2} are in GeV 2 superscript GeV 2 {\rm GeV}^{2} .

QCD variant	$Q_{\min }^2\left(\mathrm{fit}\right)$	$k$	$f_2^{p-n}(1.)$	$M^2$	${\chi }^2/d.o.f$	${\chi }_{\mathrm{ext}}^2/d.o.f$	${\chi }_{\mathrm{all}}^2/d.o.f.$
$\overline{\mathrm{MS}}$ pQCD	0.268	16.0	$-0.286\pm 0.003\pm 1.035$	$0.623\pm 0.017\pm 0.471$	7.94	35.4	130.
(F)APT	0.268	16.0	$-0.209\pm 0.002\pm 0.744$	$0.439\pm 0.012\pm 0.463$	7.48	31.5	31.2
2$\delta$	0.268	0.689	$-0.097\pm 0.002\pm 0.155$	$0.336\pm 0.023\pm 0.371$	6.25	41.1	36.9
3$\delta$	0.268	2.71	$-0.065\pm 0.001\pm 0.462$	$0.036\pm 0.017\pm 0.377$	5.99	46.6	55.0

Table 4. Table 4: As in Table 3 , but now the squared mass in the higher-twist term in Eq. ( 16 ) is running according to Eq. ( 17 ).

QCD variant	$Q_{\min }^2\left(\mathrm{fit}\right)$	$k$	$f_2^{p-n}(1.)$	$M^2\left(1{\mathrm{GeV}}^2\right)$	${\chi }^2/d.o.f$	${\chi }_{\mathrm{ext}}^2/d.o.f$	${\chi }_{\mathrm{all}}^2/d.o.f.$
$\overline{\mathrm{MS}}$ pQCD	0.268	11.9	$-0.214\pm 0.002\pm 0.651$	$0.240\pm 0.010\pm 0.272$	5.91	48.1	252.
(F)APT	0.268	16.0	$-0.171\pm 0.002\pm 0.186$	$0.152\pm 0.007\pm 0.100$	6.65	43.2	38.8
2$\delta$	0.268	9.26	$-0.160\pm 0.002\pm 0.496$	$0.277\pm 0.011\pm 0.318$	5.88	49.0	48.0
3$\delta$	0.268	2.65	$-0.064\pm 0.001\pm 0.097$	$0.010\pm 0.009\pm 0.162$	5.99	46.5	48.3

Table 5. Table 5: As in Table 2 , but now the fit is performed on the sum which includes the elastic contribution, i.e., on Γ 1 p − n ( inel . + el . ) \Gamma_{1}^{p-n}({\rm inel.+el.}) .

QCD variant	$Q_{\min }^2\left(\mathrm{fit}\right)$	$k$	$f_2^{p-n}(1.)$	${\mu }_6$	${\chi }^2/d.o.f$	${\chi }_{\mathrm{ext}}^2/d.o.f$	${\chi }_{0.268}^2/d.o.f.$	${\chi }_{\mathrm{all}}^2/d.o.f.$
$\overline{\mathrm{MS}}$ pQCD	0.66	16.	$-0.253\pm 0.002\pm 0.125$	$0.101\pm 0.001\pm 0.064$	7.11	55.6	$3.65\times {10}^5$	$1.26\times {10}^6$
(F)APT	0.66	16.	$-0.227\pm 0.002\pm 0.105$	$0.086\pm 0.001\pm 0.028$	7.53	46.3	697.	$3.54\times {10}^5$
2$\delta$	0.66	0.487	$-0.121\pm 0.002\pm 0.057$	$0.073\pm 0.001\pm 0.011$	7.18	52.0	555.	$2.60\times {10}^5$
3$\delta$	0.66	2.01	$-0.095\pm 0.002\pm 0.145$	$0.067\pm 0.001\pm 0.033$	7.04	55.8	509.	$2.27\times {10}^5$
$\overline{\mathrm{MS}}$ pQCD	0.47	1.37	$-0.194\pm 0.002\pm 0.111$	$0.114\pm 0.001\pm 0.056$	8.82	39.0	$\infty$	$\infty$
(F)APT	0.47	1.67	$-0.191\pm 0.002\pm 0.086$	$0.072\pm 0.001\pm 0.021$	19.4	63.4	364.	$2.46\times {10}^5$
2$\delta$	0.47	0.194	$-0.089\pm 0.002\pm 0.052$	$0.059\pm 0.001\pm 0.093$	14.8	57.8	250.	$1.67\times {10}^5$
3$\delta$	0.47	0.775	$-0.076\pm 0.002\pm 0.135$	$0.053\pm 0.001\pm 0.029$	12.6	44.9	193.	$1.38\times {10}^5$
$\overline{\mathrm{MS}}$ pQCD	0.268	1.82	$-0.153\pm 0.001\pm 0.081$	$0.086\pm 0.000\pm 0.010$	28.0	82.3	28.0	$\infty$
(F)APT	0.268	1.98	$-0.096\pm 0.001\pm 0.054$	$0.044\pm 0.000\pm 0.011$	131.	254.	131.	$8.84\times {10}^4$
2$\delta$	0.268	0.063	$0.018\pm 0.001\pm 0.042$	$0.028\pm 0.000\pm 0.006$	37.9	113.	37.9	$3.91\times {10}^4$
3$\delta$	0.268	0.846	$0.000\pm 0.001\pm 0.091$	$0.030\pm 0.000\pm 0.016$	68.7	145.	68.7	$4.61\times {10}^4$

Table 6. Table 6: As in Table 3 , i.e., the OPE form is “massive”, Eq. ( 16 ), but now the elastic part of BSR is included in the fit.

QCD variant	$Q_{\min }^2\left(\mathrm{fit}\right)$	$k$	$f_2^{p-n}(1.)$	$M^2$	${\chi }^2/d.o.f$	${\chi }_{\mathrm{ext}}^2/d.o.f$	${\chi }_{\mathrm{all}}^2/d.o.f.$
$\overline{\mathrm{MS}}$ pQCD	0.268	1.47	$0.117\pm 0.000\pm 0.042$	$-0.247\pm 0.000\pm 0.017$	384.	795.	$\infty$
(F)APT	0.268	3.39	$0.041\pm 0.001\pm 0.032$	$-0.207\pm 0.002\pm 0.063$	395.	760.	$4.63\times {10}^3$
2$\delta$	0.268	0.063	$0.086\pm 0.001\pm 0.044$	$-0.156\pm 0.002\pm 0.032$	96.6	245.	$1.39\times {10}^4$
3$\delta$	0.268	1.13	$0.083\pm 0.001\pm 0.082$	$-0.151\pm 0.003\pm 0.071$	143.	320.	$1.05\times {10}^4$

Table 7. Table 7: As in Table 6 , but with Q 2 superscript 𝑄 2 Q^{2} -dependent mass Eq. ( 17 ) in the higher-twist part of Eq. ( 16 ).

QCD variant	$Q_{\min }^2\left(\mathrm{fit}\right)$	$k$	$f_2^{p-n}(1.)$	$M^2\left(1{\mathrm{GeV}}^2\right)$	${\chi }^2/d.o.f$	${\chi }_{\mathrm{ext}}^2/d.o.f$	${\chi }_{\mathrm{all}}^2/d.o.f.$
$\overline{\mathrm{MS}}$ pQCD	0.268	4.88	$0.113\pm 0.001\pm 0.053$	$-0.065\pm 0.001\pm 0.070$	461.	932.	$\infty$
(F)APT	0.268	3.22	$0.052\pm 0.001\pm 0.030$	$-0.072\pm 0.001\pm 0.034$	436.	825.	$3.25\times {10}^3$
2$\delta$	0.268	0.063	$0.102\pm 0.001\pm 0.034$	$-0.051\pm 0.001\pm 0.020$	119.	288.	$6.58\times {10}^5$
3$\delta$	0.268	1.14	$0.098\pm 0.001\pm 0.044$	$-0.048\pm 0.001\pm 0.023$	167.	366.	$5.49\times {10}^4$

Table 8. Table 8: Values of the parameters of 2 δ 𝛿 \delta and 3 δ 𝛿 \delta coupling used in the present work, for N f = 3 subscript 𝑁 𝑓 3 N_{f}=3 : the Lambert Λ L . subscript Λ L \Lambda_{{\rm L.}} scale (in GeV); and the dimensionless parameters s j ≡ M j 2 / Λ L . 2 subscript 𝑠 𝑗 superscript subscript 𝑀 𝑗 2 superscript subscript Λ L 2 s_{j}\equiv M_{j}^{2}/\Lambda_{{\rm L.}}^{2} and f j ≡ ℱ j / Λ L . 2 subscript 𝑓 𝑗 subscript ℱ 𝑗 superscript subscript Λ L 2 f_{j}\equiv{\cal F}_{j}/\Lambda_{{\rm L.}}^{2} . The “input” parameter choice is α s ​ ( M Z 2 ; MS ¯ ) = 0.1185 subscript 𝛼 𝑠 superscript subscript 𝑀 𝑍 2 ¯ MS 0.1185 \alpha_{s}(M_{Z}^{2};\overline{\rm MS})=0.1185 and r τ ( D = 0 ) = 0.201 superscript subscript 𝑟 𝜏 𝐷 0 0.201 r_{\tau}^{(D=0)}=0.201 . 2 δ 𝛿 \delta coupling is in the Lambert scheme with c 2 = − 4.9 subscript 𝑐 2 4.9 c_{2}=-4.9 ; 3 δ 𝛿 \delta coupling is in the four-loop LMM scheme. We refer for details to mathprg2 ; 4l3danQCD .

$𝒜$QCD	${\overline{\alpha }}_s\left(M_Z^2\right)$	$r_{\tau }^{\left(D=0\right)}$	$s_1$	$s_2$	$s_3$	$f_1$	$f_2$	$f_3$	$s_0$	${\Lambda }_{L.}$ [GeV]	$\pi {𝒜}_{\max }$
2$\delta$	0.1185	0.201	18.734	1.0361	-	0.2929	0.5747	-	25.610	0.2564	2.0713
3$\delta$	0.1185	0.201	3.970	18.495	474.20	-2.8603	11.801	5.2543	652.	0.11564	0.9156

Equations

$${\Gamma }_1^{p-n}(Q^2)={\int }_0^{1}dx\left[g_1^p(x,Q^2)-g_1^n(x,Q^2)\right].$$

$$\Gamma_{1}^{p-n,\rm{OPE}}(Q^{2})={\Big{|}}\frac{g_{A}}{g_{V}}{\Big{|}}\frac{1}{6}(1-{\cal D}_{\rm BS}(Q^{2}))+\sum_{i=2}^{\infty}\frac{\mu_{2i}(Q^{2})}{Q^{2i-2}}\ .$$

$${\mathcal{D}}_{\mathrm{BS}}(Q^2{)}_{\mathrm{pt}}$$

$${\mathcal{D}}_{\mathrm{BS}}(Q^2{)}_{\mathcal{A}\mathrm{QCD}}$$

$${\mu }_4=\frac{M_N^2}{9}(a_2^{p-n}+4d_2^{p-n}+4f_2^{p-n}(Q^2)),$$

$$d_2^{p-n}={\int }_0^{1}dx{x}^2(2g_1^{p-n}+3g_2^{p-n}).$$

$$f_2^{p-n}(Q^2)=f_2^{p-n}(Q_{\mathrm{in}}^2){\left(\frac{a(Q^2)}{a(Q_{\mathrm{in}}^2)}\right)}^{{\gamma }_0/8{\beta }_0},$$

$$f_2^{p-n}(Q^2)=f_2^{p-n}(Q_{\mathrm{in}}^2)\left(\frac{A_{{\gamma }_0/8{\beta }_0}(Q^2)}{{\mathcal{A}}_{{\gamma }_0/8{\beta }_0}(Q_{\mathrm{in}}^2)}\right),$$

$${\Gamma }_1^{p-n,\mathrm{OPE}[4]}(Q^2;k,f_2^{p-n}(1);{\mu }_6)$$

$${\Gamma }_1^{p-n}(Q^2{)}_{el.}=\frac{1}{2}{F}_1^p(Q^2)\left(F_1^p(Q^2)+F_2^p(Q^2)\right)-\frac{1}{2}{F}_1^n(Q^2)\left(F_1^n(Q^2)+F_2^n(Q^2)\right).$$

$${\mathcal{F}}_{\tau }\equiv {\left(1+\frac{Q^2}{M_{\rho }^2}\right)}^{-1}\times ...\times {\left(1+\frac{Q^2}{(2\tau -3)M_{\rho }^2}\right)}^{-1},(\tau =2,3,\dots ),$$

$$F_1^p(Q^2)$$

$$F_1^n(Q^2)$$

$${\Gamma }_1^{p-n}(Q^2{)}_{el.}=2.3368{\left(\frac{M_{\rho }^2}{Q^2}\right)}^4+13.8851{\left(\frac{M_{\rho }^2}{Q^2}\right)}^5+\mathcal{O}(1/Q^{12}).$$

$${\chi }^2(Q_{K+1}^2\le Q^2\le Q_{K+N}^2)/d.o.f.=\frac{1}{(N-p)}\sum _{j=K+1}^{K+N}\frac{1}{{\sigma }_{j,\mathrm{stat}}^2}{\left({\Gamma }_1^{p-n,OPE[4]}(Q_j^2;k,f_2^{p-n}(1);{\mu }_6=0)-{\Gamma }_1^{p-n}(Q_j^2{)}_{\exp }\right)}^2.$$

$${\Gamma }_1^{p-n,\mathrm{mOPE}[4]}(Q^2;k,f_2^{p-n}(1);M^2)$$

$$M^2(Q^2)=\frac{m_0^2}{1+(Q^2/\mathcal{M}{)}^{1+p}}(\mathcal{M}=0.5\mathrm{GeV},p=0.1),$$

$${\Gamma }_1^{p-n}(Q^2{)}_{inel.}=\frac{{\chi }_n^2-{\chi }_p^2}{8M_N^2}{Q}^2+A(Q^2{)}^2+B(Q^2{)}^3(Q^2\lesssim 0.5{\mathrm{GeV}}^2),$$

$$\Gamma_{1}^{p-n}(Q^{2})_{\rm inel.}={\Big{|}}\frac{g_{A}}{g_{V}}{\Big{|}}\frac{1}{6}\left[1-{\mathcal{A}}^{\rm(LFH)}(Q^{2})\right]={\Big{|}}\frac{g_{A}}{g_{V}}{\Big{|}}\frac{1}{6}\left[1-\exp\left(-\frac{Q^{2}}{2\kappa^{2}}\right)\right],\qquad(Q^{2}\lesssim 1\ {\rm GeV}^{2}).$$

$$\frac{\partial }{\partial \ln {\mu }^2}a({\mu }^2)=-{\beta }_{0}a({\mu }^2{)}^2-{\beta }_{1}a({\mu }^2{)}^3-{\beta }_{2}a({\mu }^2{)}^4-\dots ,$$

$$\frac{\partial a}{\partial c_2}=a^3+\mathcal{O}(a^5),\frac{\partial a}{\partial c_3}=\frac{1}{2}{a}^4+\mathcal{O}(a^5),\frac{\partial a}{\partial c_4}=\mathcal{O}(a^5),\dots$$

$$d_1(k)$$

$$d_2(k;c_2)$$

$$d_3(k;c_2,c_3)$$

$$\mathcal{A}(Q^2)-a(Q^2)\sim {\left(\frac{{\Lambda }^2}{Q^2}\right)}^N,$$

$${\tilde{a}}_n(Q^2)\mapsto {\tilde{\mathcal{A}}}_n(Q^2),$$

$${\tilde{a}}_n(Q^2)$$

$${\tilde{\mathcal{A}}}_n(Q^2)$$

$${\tilde{a}}_n(Q^2)=a(Q^2{)}^n+\mathcal{O}(a^{n+1}).$$

$${\tilde{a}}_2$$

$${\tilde{a}}_3$$

$$a^2$$

Full text

perturbation expansion in low-energy QCD; IR-safe QCD coupling; holomorphic behavior; spacelike quantities; QCD phenomenology

††preprint: USM-TH-357

I Introduction

The polarized Bjorken sum rule (BSR) $\Gamma_{1}^{p-n}(Q^{2})$ BjorkenSR is an important spacelike QCD observable for various reasons. It is a difference of the first moment of the spin-dependent structure functions of proton and neutron, therefore its isovector nature makes it easier to describe it theoretically, in pQCD in terms of OPE, than the separate integrals of the two nucleons. Further, high quality experimental results for this quantity, obtained in polarized deep inelastic scattering (DIS), are now available in a large range of spacelike squared momenta $-q^{2}\equiv Q^{2}$: $0.054\ {\rm GeV}^{2}\leq Q^{2}<5\ {\rm GeV}^{2}$ EG1a ; EG1b ; DVCS ; E143 ; COMP1 ; COMP2 ; E155 ; SMC ; HERMES . In particular, the newer experimental results DVCS with highly reduced statistical uncertainties, extracted mainly from the Jefferson Lab CLAS EG1-DVCS experiment EG1DVCS taken on polarized targets made of protons and deuterons, render BSR an attractive quantity to test on it various extensions of pQCD to low $Q^{2}\lesssim 1\ {\rm GeV}^{2}$.

Theoretically, pQCD (with OPE), in $\overline{\rm MS}$ scheme, has been the usual approach to describe such quantities, cf. EG1a ; EG1b ; DVCS . This approach, however, has the theoretical disadvantage that the running coupling $a(Q^{2})$ [$\equiv\alpha_{s}(Q^{2})/\pi$] possesses Landau singularities at low positive $Q^{2}\lesssim 0.1\ {\rm GeV}^{2}$, and this makes it inconvenient for evaluation of spacelike observables ${\cal D}(Q^{2})$ at low $Q^{2}$, such as BSR. In recent years, an extension of pQCD couplings to low $Q^{2}$, without Landau singularities, called (Fractional) Analytic Perturbation Theory [(F)APT)] ShS ; MS96 ; ShS98 ; Sh ; BMS05 ; BKS05 ; BMS06 ; BMS10 ; Bakulev has been applied in the fitting of the theoretical OPE expression to the experimental inelastic contributions to BSR PSTSK10 , with good results.

In this work we fit the theoretical OPE expressions to the experimental BSR results in pQCD, in (F)APT, and two additional extensions of QCD to low $Q^{2}$, namely the 2$\delta$ 2danQCD ; mathprg2 and 3$\delta$ 3l3danQCD ; 4l3danQCD ${\mathcal{A}}$QCD. The latter two extensions have the coupling ${\mathcal{A}}(Q^{2})$ [the analog of the pQCD coupling $a(Q^{2})$] which is free of Landau singularities and physically motivated in the entire relevant regime of $Q^{2}$ in the complex plane, $Q^{2}\in\mathbb{C}\backslash(-\infty,-M_{\rm thr}^{2}]$, where $M_{\rm thr}^{2}\lesssim 1\ {\rm GeV}^{2}$ is a positive threshold scale. The present work is an extension of our previous work on BSR prev ; in comparison with the latter work, we now vary the $Q^{2}$-range of the fit, include the consideration of the elastic contribution, and estimate the uncertainties of the values of the extracted fit parameters due to systematic (in addition to statistical) uncertainties of the experimental BSR data.

The mentioned QCD variants (F)APT, 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD were constructed directly by imposing certain physically motivated conditions on the QCD coupling. In this context, we point out that there exist several other approaches to the construction of the QCD coupling, among them those using Dyson-Schwinger equations which involve various versions of the dynamically generated gluon mass (the mentioned 2$\delta$ ${\mathcal{A}}$QCD gives a coupling with similar properties). For a recent review of various approaches, we refer to Brodrev .

In Sec. II we present the theoretical leading-twist and higher-twist OPE contributions for the considered quantity BSR, as well as the parametrization of the elastic contribution to BSR. In Sec. III we present the results of various fits of theoretical expressions to the experimental results. Finally, in Sec. IV we summarize our conclusions.

Most of the formal aspects of the calculations are relegated to Appendices: in Appendix A we present the form of the leading-twist perturbation coefficients for a general renormalization scale and scheme; in Appendix B we present construction of ${\mathcal{A}}_{n}$, the analogs of pQCD powers $a^{n}$, in extensions of pQCD without Landau singularities; in Appendix C we summarize such extensions [(F)APT, 2$\delta$ and 3$\delta$]; in Appendix D we explain how the statistical and systematic uncertainties of the experimental data are reflected in the corresponding uncertainties of the parameters extracted in the fits; and in Appendix E we estimate the effects of the finiteness of the charm quark mass in our evaluations.

II Bjorken sum rule: theoretical expressions

The polarized Bjorken sum rule (BSR), $\Gamma_{1}^{p-n}$, is defined as the difference between proton and neutron polarized structure functions $g_{1}$ integrated over the whole $x$-Bjorken interval

[TABLE]

Based on the various measurements of these and the related structure functions, the inelatic part of the above quantity, $\Gamma_{1}^{p-n}(Q^{2})_{\rm inel.}$, has been extracted at various values of squared momenta $Q_{j}^{2}$ ($0.054\ {\rm GeV}^{2}\leq Q_{j}^{2}<5\ {\rm GeV}^{2}$): EG1a ; EG1b ; DVCS (Jefferson Lab), E143 ; E155 (SLAC), COMP1 ; COMP2 (COMPASS at CERN). 111The index $j$ in $Q_{j}^{2}$ indicates from here on that these are the scales at which the experimental values are given.

Theoretically, this quantity can be written in the Operator Product Expansion (OPE) form BjorkenSR

[TABLE]

Here, $|g_{A}/g_{V}|$ is the ratio of the nucleon axial charge, $(1-{\cal D}_{\rm BS})$ is the perturbation expansion for the leading-twist contribution, and $\mu_{2i}/Q^{2i-2}$ are the higher-twist contributions. The value obtained from neutron $\beta$ decay is known to a high accuracy, $|g_{A}/g_{V}|=1.2723\pm 0.0023$ PDG2016 , and we will use the central value. In the higher-twist terms, we will take only the terms $\sim 1/Q^{2}$ and $1/(Q^{2})^{2}$.

II.1 Perturbation expansion of the leading-twist

The leading-twist term has the canonical part ${\cal D}_{\rm BS}(Q^{2})$ whose perturbation expansion in $a\equiv\alpha_{s}/\pi$ is known up to N3LO ($\sim a^{4}$)

[TABLE]

where the bar indicates that the expansion is in the $\overline{\rm MS}$ scheme, and the renormalization scale $\mu^{2}$ is implicitly understood to be equal to the physical scale $Q^{2}$. The NLO, N2LO and N3LO coefficients ${\bar{d}}_{j}$ ($j=1,2,3$) were obtained in nloBSR ; nnloBSR ; nnnloBSR , respectively. In the considered range of momentum transfer $0<Q^{2}<5\ {\rm GeV}^{2}$, we will assume that the effective number of active quark flavors is $N_{f}=3$, and therefore only the nonsinglet (NS) contributions appear.222The singlet contribution apeears for the first time at $\sim a^{4}$, and only if $N_{f}\not=3$ Baikov:2015tea ; Larin:2013yba .

When the renormalization scale is changed from $\mu^{2}=Q^{2}$ to a general value $\mu^{2}=kQ^{2}$ ($0<k\sim 1$), and when the renormalization scheme parameters are changed from the $\overline{\rm MS}$ values ${\bar{c}}_{j}\equiv{\bar{\beta}}_{j}/\beta_{0}$ to general scheme parameter values $c_{j}$ ($j\geq 2$), the perturbation expansion changes accordingly

[TABLE]

The expressions for the new coefficients $d_{1}(k)$, $d_{2}(k;c_{2})$ and $d_{3}(k;c_{2},c_{3})$ are obtained on the basis of independence of the observable ${\cal D}_{\rm BS}(Q^{2})_{\rm pt}$ from $k$ and $c_{j}$ ($j\geq 2$), and are given in Appendix A.

In those versions of ${\mathcal{A}}$QCD where the coupling is a holomorphic function ${\mathcal{A}}(Q^{2})$ [instead of the nonholomorphic $a(Q^{2})$] with nonperturbative contributions, the power expansion (4) becomes a nonpower expansion where $a^{n}$ get replaced by ${\mathcal{A}}_{n}$ ($\not={\mathcal{A}}^{n}$)

[TABLE]

The construction of the power analogs ${\mathcal{A}}_{n}$ of $a^{n}$ were obtained in Ref. CV1 ; CV2 for integer $n$ and in Ref. GCAK for general real $n>-1$. A brief description for the construction of ${\mathcal{A}}_{n}$ is given in Appendix B. These expressions are based on the renormalization group equation (RGE), in close analogy with the RGE in the perturbation theory. The couplings ${\mathcal{A}}_{n}(Q^{2})$ can be obtained once the coupling ${\mathcal{A}}(Q^{2})$ is known. The construction of ${\mathcal{A}}(Q^{2})$ coupling is summarized in Appendix C for various variants of QCD with holomorphic coupling: (F)APT, 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD, and we refer to that Appendix for more details.

II.2 Higher-twist

In the theoretical OPE expression (2), the term with the dimension $D=2$ (i.e., $\propto 1/Q^{2}$) has the coefficient

[TABLE]

where $M_{N}\approx 0.94$ GeV is the nucleon mass, $a_{2}^{p-n}$ is the (twist-2) target mass correction, and $d_{2}^{p-n}$ is a twist-3 matrix element

[TABLE]

At $Q^{2}=1\ {\rm GeV}^{2}$, we have $a_{2}^{p-n}=0.031\pm 0.010$ and $d_{2}^{p-n}=0.008\pm 0.0036$ DVCS . We will neglect $Q^{2}$-dependence of these two quantities [as was done also in Ref. DVCS ], and will take the central values, i.e., $a_{2}^{p-n}+4d_{2}^{p-n}=0.063$. On the other hand, the coefficient $f_{2}^{p-n}(Q^{2})$ will be a parameter of the fit, and its $Q^{2}$-dependence will not be neglected, its evolution is known ShuVa ; KUY in pQCD

[TABLE]

where $\nu\equiv\gamma_{0}/(8\beta_{0})=32/81$ when $N_{f}=3$, and the reference scale will be taken $Q_{\rm in}^{2}=1\ {\rm GeV}^{2}$. In QCD with holomorphic coupling ${\mathcal{A}}(Q^{2})$, the power $a^{\gamma_{0}/8\beta_{0}}$ gets replaced by ${\mathcal{A}}_{\gamma_{0}/8\beta_{0}}$ (which is in general not equal to the simple power ${\mathcal{A}}^{\gamma_{0}/8\beta_{0}}$, as mentioned above, cf. also Appendix B,

[TABLE]

In addition, in some of the fits we will also include the $D=4$ term in the theoretical OPE expression (2) $\mu_{6}/(Q^{2})^{2}$, where we will consider the coefficient $\mu_{6}$ as $Q^{2}$-independent. Thus, our theoretical expression for BSR will be the truncated (at $D=4$, or $D=2$) OPE expression

[TABLE]

where the leading-twist truncated expressions are given in Eqs. (4) and (5), and the $Q^{2}$-dependent part of the $D=2$ term (twist-4) in Eqs. (8) and (9), for the pQCD and ${\mathcal{A}}$QCD version of QCD, respectively. We will regard the renormalization scale parameter $k\equiv\mu^{2}/Q^{2}$ and the higher-twist coefficients $f_{2}^{p-n}(Q^{2}_{\rm in})$ (with $Q^{2}_{\rm in}=1\ {\rm GeV}^{2}$) and $\mu_{6}$ as the free parameters to be determined in the fitting procedure.

In our approach, we will include in the leading-twist part of the OPE all the known terms (i.e., up to order $\sim a^{4}\sim{\mathcal{A}}_{4}$). The order of the leading-twist terms in general affects the fitted higher-twist contributions (cf. Kataev:1997nc ; Narison:2009ag ; Parente:1994bf for the fit of truncated OPE to structure functions).

II.3 Elastic contribution

The OPE evaluation is, in principle, for inclusive observables; in the case of BSR, this means that the OPE fit should be applied to the experimental values of the sum of the inelastic and elastic contribution.

The elastic contribution to BSR can be expressed Osip ; Merg in terms of the proton and neutron electromagnetic form factors $F_{1}$ and $F_{2}$

[TABLE]

The most recent parametrization of these form factors was performed in Sabbiretal from light-front holographic QCD (LFH QCD). Namely, these form factors can be expressed in terms of the inverse power expressions

[TABLE]

which are products of $\tau-1$ poles along the vector meson Regge radial trajectory in terms of the $\rho$-vector mass $M_{\rho}=0.755$ GeV and its radial excitations. We have

[TABLE]

Here: $\chi_{p}=\mu_{p}-1=1.793$ and $\chi_{n}=\mu_{n}=-1.913$ are the anomalous moments of $p$ and $n$, respectively; $\gamma_{n,p}$ are the higher Fock probabilities for the $L=0\to L=1$ spin-flip electromagnetic form factors, and $r$ is a phenomenological factor. The values of these three parameters are obtained by fitting to polarization data for the form factors Sabbiretal : $\gamma_{p}=0.27$; $\gamma_{n}=0.38$; $r=2.08$.

As a consequence, the elastic contribution (11) can be represented as a combination of inverse powers $(1+Q^{2}/M_{\rho}^{2}/(2m-1))$, and for large $Q^{2}$ this can be expanded in inverse powers of $Q^{2}$

[TABLE]

This means that at high $Q^{2}$ the elastic contribution starts with the dimension $D=8$ term $\sim(1/Q^{2})^{4}$. Theoretically, it is not included in the truncated OPE expression (10) in pQCD.333However, the last coefficient in the truncated OPE, e.g. $\mu_{6}$, is sometimes in the literature regarded to include, in a certain effective way, the contributions from the higher dimensional terms $D\geq 6$, cf. DVCS for a discussion of these aspects.*,*444The elastic contributions parametrized in a way different from that of Ref. Sabbiretal could in general contain terms of dimension $D<8$. Nonetheless, the fitted expressions of Ref. Merg give for the elastic contributions an expansion similar to Eq. (14), where the first nonzero term is a very small $D=6$ contribution which is negligible in comparison with $D=8$ term (the coefficient at $D=6$ term is by about a factor of $10^{-4}$ smaller than at $D=8$ term). All these terms there are divided by a large weakly $Q^{2}$-dependent factor $[\ln(27.81+1.72\;Q^{2}/M_{\rho}^{2})]^{4.296}$. Further, in 3$\delta$ and 2$\delta$ ${\mathcal{A}}$QCD, the higher dimensional terms up to (and including) the $D=8$ term are not included in the leading-twist contribution (5) because in these approaches ${\mathcal{A}}(Q^{2})-a(Q^{2})\sim(\Lambda^{2}/Q^{2})^{5}$ as explained in Appendix C.2 Eq. (41). Therefore, the truncated OPE expression (10) does not contain the dimension $D=8$ term $\sim(1/Q^{2})^{4}$ also in 3$\delta$ and 2$\delta$ ${\mathcal{A}}$QCD approaches. The only exception is the (F)APT where ${\mathcal{A}}(Q^{2})-a(Q^{2})\sim(\Lambda^{2}/Q^{2})^{1}$, cf. Appendix C.1, and where the truncated OPE (10) could contain, in principle, all the elastic contributions, including the $\sim(1/Q^{2})^{4}$ term.

Therefore, it looks more natural to first fit the truncated OPE expression (10), with pQCD and (3$\delta$ and 2$\delta$) ${\mathcal{A}}$QCD, to the BSR experimental results EG1a ; EG1b ; DVCS ; E143 ; COMP2 which are obtained for the inelastic contribution; and after such a fit, add the elastic contribution (11) [with the parametrization (12)-(13b)]. We note that the theoretical expression obtained for the total BSR in this way is the sum of the expressions (10) and (11), which is again an OPE expression as it should be for such an inclusive spacelike observable as the total BSR.

The other approach would be to fit the truncated OPE expression (10) with the experimental points for the total BSR, i.e., fit with the experimental points for inelastic BSR with the elastic contribution added to them. Such an approach seems less natural, at least in pQCD and 3$\delta$ and 2$\delta$ ${\mathcal{A}}$QCD, because the truncated expression (10) in these cases in principle does not contain the leading elastic contribution $\sim 1/(Q^{2})^{4}$, cf. Eq. (41) in Appendix C.2.

Nonetheless, below we will apply, for completeness, both approaches in our numerical fitting procedures.

III Numerical fits

In the numerical fits, we will consider the following parameters to be fitted in the truncated OPE expression (10): (i) the renormalization scale parameter $k\equiv\mu^{2}/Q^{2}$ of the (truncated) leading-twist contribution (4) [for pQCD] or (5) [for ${\mathcal{A}}$QCD]; (ii) the parameter $f_{2}^{p-n}(1\ {\rm GeV}^{2})$ appearing in the $D=2$ term of Eq. (10) [cf. Eqs. (8) and (9)]; (iii) and in some fits the parameter $\mu_{6}$ in the $D=4$ term of Eq. (10) will be included. We will take for the experimental data points for the inelastic contribution to BSR the data from EG1a ; EG1b ; DVCS ; E143 ; COMP2 ;555 We will not take into account the SLAC E155, SMC and HERMES points E155 ; SMC ; HERMES , at $Q^{2}=5\ {\rm GeV}^{2}$, $10\ {\rm GeV}^{2}$, and $2.5$ (and $5$) ${\rm GeV}^{2}$ , respectively, as they were obtained by NLO pQCD $Q^{2}$-evolution from data distributed across a wide range of $Q^{2}$. Among the points of SLAC E143 we exclude the point at $Q^{2}=5\ {\rm GeV}^{2}$ for the same reason. they are in the momentum interval $0.054\ {\rm GeV}^{2}\leq Q_{j}^{2}<5\ {\rm GeV}^{2}$. Among these data points, we exclude from the fit the four points with $Q_{j}^{2}\geq 3\ {\rm GeV}^{2}$ (three from DVCS , at $Q_{j}^{2}=3.223$, $3.871$ and $4.739\ {\rm GeV}^{2}$; and one from COMP2 at $Q^{2}=3\ {\rm GeV}^{2}$) because they tend to decrease the quality parameter $\chi^{2}/{\rm(d.o.f.)}$ of the fit significantly. Nonetheless, we will show both quality parameters $\chi^{2}/{\rm(d.o.f.)}$ obtained from the mentioned fits, i.e., the one with the four points excluded, and included.666If we perform the fit with the four points included, the results differ somewhat, but not significantly. The number of active flavors will be kept all the time at $N_{f}=3$. The fits will be performed by the least squares method, taking into account the statistical uncertainties $\sigma_{j,{\rm stat}}$ in the points $Q_{j}^{2}$, we will consider them independent of each other. These uncertainties $\sigma_{j,{\rm stat}}$ are in general significantly smaller than the systematic uncertainties $\sigma_{j,{\rm sys}}$. The latter are strongly correlated, and we will consider them as completely correlated. The uncertainties in the values of the extracted parameters $k\equiv\mu^{2}/Q^{2}$, $f_{2}^{p-n}(Q^{2}_{\rm in})$ and $\mu_{6}$ are then due to statistical (small) and systematic (larger) uncertainties of the data. We refer to Appendix D on how we obtained the uncertainties of the extracted values of the fit parameters.

Using specific values of the renormalization scale parameter $k\equiv\mu^{2}/Q^{2}$ may allow us to incorporate in our evaluations (5) at least a part of the contributions of higher orders of the series. It is expected that $k\sim 1$, and usually it is taken in the literature in the range $1/2<k<2$, sometimes $1/4<k<4$. In the considered work, after replacing the powers of the pQCD coupling by their analytic counterparts, cf. Eqs. (4)-(5), a spacelike observable depends usually weakly on the contributions of higher orders. Still, in order not to miss the possibly relevant influence of higher orders, we decided to increase the range of possible $k$ values to: $1/16<k<16$. This choice still avoids very large or very small renormalization scales where the corresponding coefficients $d_{n}(k)$ at the powers $a(kQ^{2})^{n+1}$, or at their analytic counterparts ${\mathcal{A}}_{n+1}(kQ^{2})$, contain powers of the large terms $\sim\ln k$ [cf. Eqs. (22)] which may destroy the convergence of the series already at low $n$. We will see that the results of the fitting will rarely give the extreme values $k=16$ or $k=1/16\approx 0.063$, mostly in the cases of ($\overline{\rm MS}$) pQCD and (F)APT.

We present below the results of the fits of the inelastic contributions to the truncated OPE expression (10), first in Sec. III.1 with $\mu_{6}=0$, and then in Sec. III.2 for $\mu_{6}\not=0$ included as a fit parameter; the elastic contribution is added afterwards. In Sec. III.3 we present the corresponding fits for the case when the higher-twist term in the OPE has a mass parameter. In Sec. III.4 we include the fits for two different ansätze for BSR at very low $Q^{2}$. Finally, in Sec. III.5 we present the fits of the truncated OPE expression (10), with $\mu_{6}\not=0$ included as a fit parameter, applied to the total BSR, i.e., to the sum of the inelastic data points (at $Q_{j}^{2}$) and the corresponding elastic contribution $\Gamma_{1}^{p-n}(Q_{j}^{2})_{\rm el.}$.

In each case, the fits are performed for four variants of QCD: in pQCD (in $\overline{\rm MS}$), in (F)APT (in $\overline{\rm MS}$), and in 2$\delta$ mathprg2 and 3$\delta$ ${\mathcal{A}}$QCD 4l3danQCD . Each of these fits is performed by excluding among the data points those with $Q^{2}<Q^{2}_{\rm min}$, where $Q^{2}_{\rm min}=0.268$ or $0.66\ {\rm GeV}^{2}$ (and sometimes also: $Q^{2}_{\rm min}=0.47\ {\rm GeV}^{2}$). As mentioned earlier, the four points at high $Q^{2}\geq 3\ {\rm GeV}^{2}$ are excluded from the fit as well.

It is reasonable to assume that the number of active quark flavors in the considered interval $0.054\ {\rm GeV}^{2}<Q^{2}<5\ {\rm GeV}^{2}$ is $N_{f}=3$, cf. Appendix E. All the fits will be performed in $\overline{\rm MS}$ pQCD approach, in the Fractional Analytic Perturbation Theory [(F)APT], and in 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD approach, cf. Appendix C where also the values of the parameters of these three QCD variants are presented. In $\overline{\rm MS}$ pQCD and in 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD the (underlying) pQCD running coupling $a(Q^{2})$ is determined by the requirement $\pi a(M_{Z}^{2};\overline{\rm MS})=0.1185$ PDG2014 ; PDG2016 , and in (F)APT we use ${\overline{\Lambda}}_{3}=0.45$ GeV (at the end of Sec. III.2 we also comment on the case ${\overline{\Lambda}}_{3}=0.40$ GeV); we refer to Appendix C for more details.

III.1 Fits with $\mu_{6}=0$

In Figs. 1(a),(b), we present the curves in the mentioned four QCD variants, for $Q^{2}_{\rm min}=0.66$ and $0.268\ {\rm GeV}^{2}$, respectively.

The corresponding results and four fit quality parameters $\chi^{2}/{\rm d.o.f}$ are given in Table 1, for $Q^{2}_{\rm min}=0.66$, $0.47$ and $0.268\ {\rm GeV}^{2}$ .

As mentioned earlier, the fit is performed for the data in the momentum interval $Q^{2}_{\rm min}\leq Q_{j}^{2}\leq 3\ {\rm GeV}^{2}$ (with the COMPASS data point COMP2 excluded), and the corresponding quality parameter for that interval is denoted as $\chi^{2}/{\rm d.o.f}$. In addition, $\chi^{2}_{\rm ext}/{\rm d.o.f.}$ is the quality parameter for the wider interval $Q^{2}_{\rm min}\leq Q_{j}^{2}<5\ {\rm GeV}^{2}$; $\chi^{2}_{0.268}/{\rm d.o.f.}$ is the parameter when $0.268\ {\rm GeV}^{2}\leq Q_{j}^{2}\leq 3\ {\rm GeV}^{2}$; and $\chi^{2}_{\rm all}/{\rm d.o.f.}$ is the parameter when all the experimental points EG1a ; EG1b ; DVCS ; E143 ; COMP2 are included, $0.054\ {\rm GeV}^{2}\leq Q_{j}^{2}<5\ {\rm GeV}^{2}$. In these quantities $\chi^{2}$, the central values of $k$ and $f_{2}^{p-n}(Q^{2}_{\rm in})$ ($Q^{2}_{\rm in}=1\ {\rm GeV}^{2}$) parameters obtained from the mentioned fit interval $Q^{2}_{\rm min}\leq Q_{j}^{2}\leq 3\ {\rm GeV}^{2}$ are used; we take here “d.o.f” as $(N-p)$ where $N$ is the number of fitted data points entering the considered $\chi^{2}$ and $p$ is the number of parameters of the fit ($p=2$ here)

[TABLE]

The values of these quantities are always large, in the best cases between 1 and 10. This is so because the statistical uncertainties of the newer JLAB data DVCS are very small, $\sigma_{j,{\rm stat}}\lesssim 10^{-3}$, and simultaneously, our fit function (truncated OPE) is different from the ideal function which we do not know. Nonetheless, they decrease when analytic QCD variants are employed, especially 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD.

We wish to point out that the approach of ($\overline{\rm MS}$) QCD in the case of $Q_{\rm min}^{2}=0.268\ {\rm GeV}^{2}$ is, in principle, not applicable. This is so because the corresponding coupling $a(Q^{2})$ has a Landau branching point at $Q^{2}_{\rm branch}=0.371\ {\rm GeV}^{2}$, which makes the $D=2$ running coefficient $f_{2}^{p-n}(Q^{2})$, Eq. (8), undefined at $Q^{2}\leq 0.371\ {\rm GeV}^{2}$. Nonetheless, in order to be able to present a curve, we applied in the fitting case $Q_{\rm min}^{2}=0.268\ {\rm GeV}^{2}$ in the $\overline{\rm MS}$ pQCD approach a restriction on the leading-twist renormalization scale $\mu^{2}=kQ^{2}$, namely $k>1.383$; and this not just in the leading-twist contribution, but we also made an ad hoc replacement in the $D=2$ running coefficient $f_{2}^{p-n}(Q^{2})$, Eq. (8): $f_{2}^{p-n}(Q^{2})\mapsto f_{2}^{p-n}(kQ^{2})$. We applied this also in the case when $\mu_{6}\not=0$ (Sec. III.2, for $\overline{\rm MS}$ pQCD approach with $Q_{\rm min}^{2}=0.268\ {\rm GeV}^{2}$). In other approaches (APT and ${\mathcal{A}}$QCD’s) this is not necessary, as there are no Landau singularities. Further, we note that in the last column in Table 1 we have $\chi^{2}/{\rm d.o.f.}=\infty$ in the case of $\overline{\rm MS}$ pQCD for $Q^{2}_{\rm min}=0.66$ and $0.47\ {\rm GeV}^{2}$. This is so because at the lowest available experimental point, $Q_{j=1}^{2}=0.054\ {\rm GeV}^{2}$, the coupling is $a(kQ_{1}^{2};\overline{\rm MS})=\infty$ because $kQ_{1}^{2}<Q^{2}_{\rm branch}=0.371\ {\rm GeV}^{2}$ and thus we hit Landau singularities there.

We can deduce from Figs. 1 and Table 1: (a) the best results in the considered approach ($\mu_{6}=0$) are obtained in 3$\delta$ ${\mathcal{A}}$QCD; (b) the quality of extrapolation of the obtained fitted curves from the fitting interval, $Q^{2}_{\rm min}\leq Q^{2}\leq 3\ {\rm GeV}^{2}$, to the entire interval, $0.054\ {\rm GeV}^{2}\leq Q^{2}<5\ {\rm GeV}^{2}$, does not improve significantly when the fitted interval is extended (i.e., when $Q^{2}_{\rm min}$ is lowered), cf. also the last column in Table 1. This indicates that the behavior of the curves in the extrapolated regions, $0.054\ {\rm GeV}^{2}<Q^{2}<0.268\ {\rm GeV}^{2}$ and $3\ {\rm GeV}^{2}<Q^{2}\leq 5\ {\rm GeV}^{2}$, is of similar quality in the cases of different values of $Q^{2}_{\rm min}({\rm fit})$, and gives the dominant part of $\chi^{2}_{\rm all}/{\rm d.o.f.}$. The same can be observed for the values of $f_{2}^{p-n}(1)$: they do not depend much on the value of $Q^{2}_{\rm min}({\rm fit})$, but only on the QCD variant used in the fit. This appears to be related with the fact that only one parameter beyond the leading-twist contribution is used here ($\mu_{4}\leftrightarrow f_{2}^{p-n}$), representing a truncated OPE ansatz with a significantly restricted freedom.

III.2 Fits with $\mu_{6}\not=0$

When we include $\mu_{6}$ in the fit as the third parameter, the resulting curves are presented in Figs. 2(a),(b), for $Q^{2}_{\rm min}=0.66$ and $0.47\ {\rm GeV}^{2}$, respectively. Further, when $Q^{2}_{\rm min}=0.268\ {\rm GeV}^{2}$, the results are shown in Figs. 3(a),(b), at the higher $Q^{2}$ and the lower $Q^{2}<1\ {\rm GeV}^{2}$ momenta, respectively. The corresponding results are given in Table 2, for $Q^{2}_{\rm min}=0.66$, $0.47$ and $0.268\ {\rm GeV}^{2}$, with the same notations as in Table 1.

The various versions of $\chi^{2}/{\rm d.o.f.}$ are those as explained in the previous Sec. III.1, except that now in the relation (15) the factor in front of the sum is $1/(N-3)$ ($p=3$, d.o.f. is $N-3$).

Comparing Table 2 and Figs. 2-3 with Table 1 and the Figs. 1 of the previous Sec. III.1 where $\mu_{6}=0$ was kept, we can see that in the cases of $Q^{2}({\rm fit})\geq 0.66\ {\rm GeV}^{2}$ and $Q^{2}({\rm fit})\geq 0.47\ {\rm GeV}^{2}$ the $\mu_{6}=0$ fits give in general better $\chi^{2}_{\rm all}/{\rm d.o.f.}$ (the last columns of Tables 1 and 2). This means that the extrapolation down to the lowest experimental point $Q^{2}=Q_{j=1}^{2}=0.054\ {\rm GeV^{2}}$ is better when $\mu_{6}=0$ (with the exception of $\overline{\rm MS}$ pQCD case where problems with Landau singularities appear). This occurs because the inclusion of the $\mu_{6}/(Q^{2})^{2}$ term in the truncated OPE makes this expression less stable at very low $Q^{2}$.

On the other hand, when $Q^{2}({\rm fit})\geq 0.268\ {\rm GeV}^{2}$, some of the fits (2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD) give better extrapolation when $\mu_{6}\not=0$. Table 2 and Figs. 3 also show that when $\mu_{6}\not=0$ and $Q^{2}({\rm fit})\geq 0.268\ {\rm GeV}^{2}$, the best extrapolation to low $Q^{2}$ is obtained in 3$\delta$ ${\mathcal{A}}$QCD, followed by 2$\delta$ ${\mathcal{A}}$QCD. As Fig. 3(b) suggests, the fitted curve in pQCD $\overline{\rm MS}$ extrapolated to low $Q^{2}$ appears to be almost as good; in this case, however, we should keep in mind that the renormalization scale is $\mu^{2}=kQ^{2}$ ($k=3.00$) and that this scale was used also in $f^{p-n}_{2}(Q^{2})$, i.e., the ad hoc replacement $f^{p-n}_{2}(Q^{2})\mapsto f^{p-n}_{2}(kQ^{2})$ was performed in order to avoid the Landau singularities in the $D=2$ term at $Q^{2}\geq 0.268\ {\rm GeV}^{2}$ (cf. also the discussion about that point in Sec. III.1). Despite this replacement, in $\chi^{2}_{\rm all}/{\rm d.o.f.}$ Landau singularities are hit, because the three lowest experimental points give $kQ_{j}^{2}<Q^{2}_{\rm branch}=0.371\ {\rm GeV}^{2}$ ($j=1,2,3$) and are thus in the Landau singularity region.777We have $Q_{1}^{2}=0.054\ {\rm GeV}^{2}$, $Q_{2}^{2}=0.078\ {\rm GeV}^{2}$, and $Q_{3}^{2}=0.101\ {\rm GeV}^{2}$.

In the case of (F)APT, in contrast to 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD, the coupling ${\mathcal{A}}(Q^{2})$ differs from the underlying pQCD coupling $a(Q^{2})$ nonnegligibly at high $|Q^{2}|>1\ {\rm GeV}^{2}$ [the index in Eq. (23) is $N=1$ in (F)APT; $N=5$ in 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD]. This implies that (F)APT has a certain ambiguity when “normalizing” the strength of ${\mathcal{A}}(Q^{2})$. As mentioned, we fixed the strength of the coupling ${\mathcal{A}}(Q^{2})$ in (F)APT to the value ${\overline{\Lambda}}_{N_{f}=3}=0.45$ GeV, since in such a case (F)APT reproduces approximately the QCD phenomenology at high energies (cf. also Sh ). A question appears whether the results of our fits in (F)APT depend significantly on this value. We repeated the analysis in (F)APT with the value ${\overline{\Lambda}}_{N_{f}=3}=0.40$ GeV, and it turned out that the results of the fits did not change significantly. For example, when $Q^{2}_{\rm min}=0.268\ {\rm GeV}^{2}$ and $\mu_{6}=0$, we obtained $k=2.86$ and the central value $f_{2}^{p-n}(1.)=-0.137$, and for the three quality parameters $\chi^{2}_{\rm ext}/{\rm d.o.f.}$, $\chi^{2}_{0.268}/{\rm d.o.f.}$ and $\chi^{2}_{\rm all}/{\rm d.o.f}$ the values $120.$, $54.7$ and $878.$, respectively (to be compared with the corresponding values in Table 1, the third line from below). When $\mu_{6}\not=0$ in the fit, we obtained $k=2.62$, $f_{2}^{p-n}(1.)=-0.187$, $\mu_{6}=0.013\ {\rm GeV}^{4}$ and for the mentioned three quality parameters the values $42.3$, $7.94$ and $4.43\times 10^{3}$, respectively (to be compared with the corresponding values in Table 2, the third line from below).

When we add the parametrized elastic contribution of BSR, Eqs. (11)-(13), to the experimental points and to the theoretical curves of Figs. 3, we obtain the results presented in Figs. 4(a),(b). In comparison with Figs. 3, the values of BSR are shifted to significantly higher values at low $Q^{2}$.

With this approach, the quality of fits ($\chi^{2}/{\rm d.o.f}$) does not change, as we consider the elastic contribution as known (and parametrized) and added here simultaneously to the (inelastic) theoretical fitting curves and to the data points. This means that the results of Table 2 remain unchanged under this subsequent addition of $\Gamma_{1}^{p-n}(Q^{2})_{\rm el.}$.

We can observe in the results of Tables 2 and 1 that the values of the higher-twist parameters, $f_{2}^{p-n}(1)$ and $\mu_{6}$, are in the analytic variants of QCD smaller than in pQCD, this reduction being especially strong in the 3$\delta$ QCD variant. It has been noted in the literature that in pQCD OPE there is a duality between the order of truncation of the leading-twist series and the higher-twist contribution PSTSK10 ; Kataev:1997nc ; Narison:2009ag ; Parente:1994bf : higher-twist contribution often significantly decreases with the inclusion of higher orders in the leading-twist part. This effect and ambiguity become stronger in the ranges where the perturbation theory becomes questionable (for example, at the large and low values of the Bjorken variable $x$, as it was shown in Refs. Kotikov:1992ht ; Krivokhizhin:2005pt , respectively). It has been observed that the higher-twist contribution is smaller, but also more stable (under the inclusion of more terms in the leading-twist), in QCD variants with infrared modifications of the coupling (various modifications lead to quite similar results Kotikov:2004uf ). The latter probably incorporate a part of the higher-twist contributions (which are rather cumbersome Illarionov:2004nw ) into (formally) the leading-twist contribution for small $x$ range at moderately small $Q^{2}$ values ($\lesssim 1\ {\rm GeV}^{2}$) (see Ref. Cvetic:2009kw and more recent studies Kotikov:2012sm of the precise combined H1 and ZEUS data Aaron:2009aa for the DIS structure function $F_{2}$).

Following the above observations in Tables 2 and 1, we can conclude that the applications of the 2$\delta$ and especially 3$\delta$ ${\mathcal{A}}$QCD are very appropriate frameworks for the BSR studies because they appear to resum effectively a large part of the perturbative contribution into the leading-twist part (5). In this context, we wish to recall that 3$\delta$ ${\mathcal{A}}$QCD is significantly different from the other two ${\mathcal{A}}$QCD variants [(F)APT and 2$\delta$ ${\mathcal{A}}$QCD] in the infrared region, because its coupling is not just finite there but goes to zero, ${\mathcal{A}}(Q^{2})\sim Q^{2}\to 0$, as motivated by large-volume lattice calculations, cf. LattcoupNf0 ; LattcoupNf0b ; LattcoupNf2 ; LattcoupNf4 and Appendix C.2. Further, we wish to point out that the 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD couplings ${\mathcal{A}}$ (and thus ${\widetilde{\mathcal{A}}}_{n}$ and ${\mathcal{A}}_{n}$) are at large $Q^{2}$ indistinguishable from their underlying pQCD couplings $a$ (${\widetilde{a}}_{n}$, $a^{n}$), cf. Eq. (41), in contrast to (F)APT which satisfies the relation Eq. (23) with $N=1$. Therefore, theoretically, neither the higher-twist contribution of order $1/(Q^{2})^{N}$ with $N\leq 4$ ($D\leq 8$), nor a part of it, is incorporated in the leading-twist contribution (5) in the 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD, in contrast to (F)APT. This indicates that the higher-twist terms extracted here with 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD (in truncated OPE) represent an effective form for the true higher-twist contribution with dimension $D\leq 8$ [and a part of the other ($D\geq 10$) presumably small contribution]. In pQCD this is definitely not so, because of the mentioned duality there between the order of truncation of the leading-twist series and the extracted higher-twist contribution. As a consequence, the extracted effective higher-twist contribution in pQCD represents a sum of the true higher-twist contribution and a significant part of the perturbative (leading-twist) contribution; this effective higher-twist contribution appears to be in general larger than the true higher-twist contribution.

III.3 Fits with “massive” OPE

For comparison, we performed a similar fit, but now with a “massive” higher-twist term instead of the truncated OPE expression (10)

[TABLE]

where the squared mass $M^{2}$ in the denominator of the higher-twist part888Similar higher-twist expressions were used in the analyses of BSR in KTG1 ; KTG2 where the leading-twist contribution was evaluated with the “Massive” Perturbation Theory (MPT) MPT . MPT is an extension of APT ShS ; MS96 ; ShS98 ; Sh where, in contrast to APT, the coupling ${\mathcal{A}}(Q^{2})$ is analytic in the point $Q^{2}=0$; nonetheless, the index $N$ of Eq. (23) (Appendix B) remains in MPT at the minimal value as in APT, i.e., $N=1$. is taken to be constant (not running), and is expected to be $0<M^{2}\lesssim 1\ {\rm GeV}^{2}$. Now, instead of $f^{p-n}_{2}(1)$ and $\mu_{6}$, the fit parameters are $f_{2}^{p-n}(1)$ and $M^{2}$. The resulting curves, for $Q^{2}_{\rm min}=0.268\ {\rm GeV}^{2}$, are given in Figs. 5(a),(b), at the higher $Q^{2}$ and the lower $Q^{2}<1\ {\rm GeV}^{2}$ momenta, respectively.

The corresponding results are given in Table 3. These curves are analogous to those in the previous Figs. 3(a),(b) in which the truncated OPE (10) was used. Numerically, the behavior at low $Q^{2}$ in the “massive” case is significantly influenced by the $Q^{2}$-dependence of $f_{2}^{p-n}(Q^{2})$. In the $\overline{\rm MS}$ pQCD case, as in the $\overline{\rm MS}$ pQCD cases of the analyses in all the Sections, we replaced in the higher-twist running parameter $f_{2}^{p-n}(Q^{2})$ [cf. Eq. (8)] the scale $Q^{2}$ in an ad hoc way by the renormalization scale $kQ^{2}$ used in the leading-twist part ($k=16$ resulted here), in order to artificially avoid the problem of Landau singularities in the pQCD coupling $a(Q^{2})$. Comparing Figs. 5 and Table 3 with the corresponding “nonmassive” case Figs. 3 and Table 2, we see that the results and extrapolations in the case of 3$\delta$ ${\mathcal{A}}$QCD are now comparably good in the “massive” and the $\mu_{4}\&\mu_{6}$ approaches. Stated differently, the corresponding $\chi^{2}_{\{...\}}/{\rm d.o.f.}$ values in Tables 3 and 2 (with $Q^{2}_{\rm fit}=0.268\ {\rm GeV}^{2}$) are very similar. In the case of 2$\delta$ ${\mathcal{A}}$QCD and (F)APT, the extrapolations are better in the “massive” than in the $\mu_{4}\&\mu_{6}$ approach, i.e., $\chi^{2}_{\rm all}/{\rm d.o.f.}$ is significantly reduced in the “massive” case. One reason for this lies perhaps in the fact that the massive higher-twist term is under control at very low $Q^{2}$, unlike the separate $\mu_{4}(Q^{2})/Q^{2}$ and $\mu_{6}/(Q^{2})^{2}$ terms. Further, the extracted values of $f_{2}^{p-n}(1)$ are in general similar in the $\mu_{6}=0$, $\mu_{4}\&\mu_{6}$ and the “massive” approaches, although the uncertainties of the extracted parameters are quite high in the “massive” approach.999The systematic uncertainties of the extracted parameters in Table 3 are large and should therefore be regarded as crude estimates only, cf. comments in Appendix D.3.

The results of Table 3 show that the QCD variants with infrared-finite analytic coupling, and especially 3$\delta$ ${\mathcal{A}}$QCD, give smaller values of higher-twist parameters $f_{2}^{p-n}(1)$ and $M^{2}$ than pQCD. On the one hand, at low $Q^{2}$, the smaller values of $M^{2}$ compensate partially the decreased value of $f_{2}^{p-n}$ in the higher-twist contribution. On the other hand, smaller values of $M^{2}$ and $f_{2}^{p-n}(1)$ mean that at higher values of $Q^{2}$ the higher-twist contribution is significantly reduced; this can be seen also by expanding the massive higher-twist term of Eq. (16) in powers of $M^{2}/Q^{2}$. Such effect is in full agreement with one observed in Ref.PSTSK10 ; the effect can be considered as a stabilization of the higher-twist contribution. We also recall that 3$\delta$ ${\mathcal{A}}$QCD is significantly different from (F)APT and 2$\delta$ ${\mathcal{A}}$QCD in the infrared region, since its coupling goes to zero there, ${\mathcal{A}}(Q^{2})\sim Q^{2}\to 0$.

It is possible to choose for the higher-twist term a massive form with a running mass, in the spirit of a dynamical effective gluon mass of the gluon propagator at low $Q^{2}$ Cornwall . Such masses appear in the literature often in definitions of QCD couplings at low $Q^{2}$, and are responsible for the freezing (finiteness) of the coupling at $Q^{2}\to 0$. The couplings in (F)APT and 2$\delta$ ${\mathcal{A}}$QCD have a freezing which could be described also via a running effective gluon mass. Our view is that the OPE higher-twist terms represent a new contribution not contained in the QCD coupling itself. The squared mass $M^{2}$ in such terms, Eq. (16), is considered constant, in the spirit of the approach of Ref. KTG1 (cf. also KTG2 ), where the basic component (delta function) in the spectral function of the higher-twist contribution gives such a mass term. Nonetheless, we repeated the aforementioned analysis for the case of a running squared mass $M^{2}(Q^{2})$ representative of a dynamical effective gluon mass, chosen with a simple parametrization of Ref. ABP

[TABLE]

where we chose for the parameters ${\cal M}$ and $p$ values within the expected regions ABP . The adjustable squared mass scale was taken (instead of $m_{0}^{2}$) to be $M^{2}(1\ {\rm GeV}^{2})$. The same analysis then gave the results presented in Table 4. We can see that the results are qualitatively similar to those with the constant squared mass $M^{2}$, Table 3, except that the values of $M^{2}(1\ {\rm GeV}^{2})$ are now lower. Nonetheless, a reasonable definition of the average squared mass value in the present analysis, for the considered fit interval $0.268\ {\rm GeV}^{2}<Q^{2}<3\ {\rm GeV}^{2}$ may be $\langle M^{2}(Q^{2})\rangle=(1/2)\times(M^{2}(0.268)+M^{2}(3))\approx 1.52\times M^{2}(1\ {\rm GeV}^{2})$. Or, if we regard that the squared mass is most relevant only in the low-$Q^{2}$ part $0.268\ {\rm GeV}^{2}<Q^{2}<1\ {\rm GeV}^{2}$ of the fit interval, a reasonable definition of the average squared mass would be $\langle M^{2}(Q^{2})\rangle=(1/2)\times(M^{2}(0.268)+M^{2}(1))\approx 1.85\times M^{2}(1\ {\rm GeV}^{2})$.

III.4 Testing low-$Q^{2}$ regime ansätze

At low $Q^{2}$, the inelastic contribution to BSR behaves as $\sim Q^{2}$, according to Gerasimov-Drell-Hearn sum rule GDHsr as pointed out and used in Ansel ; PSTSK10 ; GDHlow . Based on this, an expansion EG1b motivated by chiral perturbation theory ($\chi$PT) can be constructed

[TABLE]

where, according to Gerasimov-Drell-Hearn sum rule GDHsr , $\chi_{n}$ and $\chi_{p}$ are anomalous magnetic moments of nucleons [which appear also in the elastic BSR contributions, cf. Eqs. (13)]; the parameters $A$ and $B$ are determined in the fit. When we fit with this expression the inelastic BSR data EG1a ; EG1b ; DVCS for $Q^{2}\leq Q^{2}_{\rm max}=0.2\ {\rm GeV}^{2}$, we obtain $A=0.765\ {\rm GeV}^{-2}$ and $B=0.678\ {\rm GeV}^{-4}$, with $\chi^{2}(Q^{2}\leq Q^{2}_{\rm max})/{\rm d.o.f.}=0.720$. On the other hand, if we take $Q^{2}_{\rm max}=0.5\ {\rm GeV}^{2}$ in the fitting, we obtain $A=0.744\ {\rm GeV}^{-2}$ and $B=-1.033\ {\rm GeV}^{-4}$, $\chi^{2}(Q^{2}\leq Q^{2}_{\rm max})/{\rm d.o.f.}=1.313$.

Another possible ansatz for the inelastic contribution to BSR at low $Q^{2}$ is the form of the light-front holographic (LFH) effective charge ${\mathcal{A}}^{\rm(LFH)}$ in the BSR ($g_{1}$) scheme LFH ; LFHBSR [${\mathcal{A}}(0)_{g_{1}}=1$]

[TABLE]

When we fit with this expression the inelastic BSR data EG1a ; EG1b ; DVCS for $Q^{2}<Q^{2}_{\rm max}=0.5\ {\rm GeV}^{2}$, we obtain $\kappa=0.479$ GeV.101010This is not far from the universal nonperturbative scale $\kappa=M_{\rho}/\sqrt{2}=0.548$ GeV, cf. LFHBSR .

In Figs. 6(a),(b), we present these low-$Q^{2}$ expressions. The $\overline{\rm MS}$ pQCD and 3$\delta$ ${\mathcal{A}}$QCD curves (obtained from fit with $Q^{2}\geq Q^{2}_{\rm min}=0.268\ {\rm GeV}^{2}$) are also included for comparison.

In Fig. 6(b) we can see that the theoretical fitted curves of the presented QCD variants connect smoothly with the nonperturbative low-$Q^{2}$ curves ($\chi$PT and LFH, both fitted to BSR up to $Q^{2}=0.5\ {\rm GeV}^{2}$) at the values of $Q^{2}$ around $0.3\ {\rm GeV}^{2}$. Nonetheless, we recall that the apparent success of the ($\overline{\rm MS}$) pQCD curve, down to about $0.3\ {\rm GeV}^{2}$, was achieved due to the ad hoc change of scale in $f^{p-n}_{2}(Q^{2})$ to $kQ^{2}$ (with $k=3$), to avoid the Landau singularities in the $D=2$ term at $Q^{2}\geq 0.268\ {\rm GeV}^{2}$ (cf. also the discussion about that point in Sec. III.2). We further notice that the 3$\delta$ QCD curve agrees well with the mentioned $\chi$PT curve in a broader interval, $0.17\ {\rm GeV}^{2}<Q^{2}<0.3\ {\rm GeV}^{2}$. Similar analyses with the goal of connecting the curves of pQCD (or of specific QCD variants) with nonperturbative curves at low $Q^{2}$ were performed in some of the references PSTSK10 , in LFHmatch , and was discussed also in Burkert .

In Refs. Bernard ; Lensky , BSR at low $Q^{2}<0.3\ {\rm GeV}^{2}$ was calculated with baryon chiral perturbation theory (B$\chi$PT) up to NLO. The authors of Bernard ; Lensky did not give the values of their parameters for BSR. However, careful visual comparison of their obtained $Q^{2}$-dependence [cf. curves and bands in their Fig. 6(b) in Ref. Lensky ] with our $\chi$PT curve (18) with $Q^{2}_{\rm fit}\leq 0.5\ {\rm GeV}^{2}$ [the long-dashed curve in our Fig. 6(b)] shows very good agreement between them.

III.5 Fitting to total (inelastic $+$ elastic) BSR data

We perform also the fitting of the truncated (at $D=4$) OPE expression (10) to the data for the total BSR $\Gamma_{1}^{p-n}(Q^{2})_{\rm inel.+el.}$. As argued at the end of Sec. II.3 [after Eq. (14)], such a fit has problematic aspects. These data are obtained by adding to the experimental data $\Gamma_{1}^{p-n}(Q_{j}^{2})_{\rm inel.}$ of Refs. EG1a ; EG1b ; DVCS the parametrized elastic contribution Eqs. (11)-(13) of Ref. Sabbiretal . The uncertainties of the parametrized elastic form factors are considered to be less than $10\%$ Sabbiretal . Since the elastic contribution (11) is quadratic in the form factors (and is numerically dominated by the proton contribution), the relative uncertainties of the elastic contribution are less than $5\%$. Since we do not have more information about these uncertainties, we will neglect them in this analysis. This means that in the total BSR, we will consider that the statistical and the systematic uncertainties $\sigma_{j,{\rm stat}}$ and $\sigma_{j,{\rm sys}}$ are those of the inelastic contribution. Otherwise, the fitting is performed as in the previous Sec. III.2.

The resulting curves are presented in Figs. 7(a),(b), for $Q^{2}\geq Q^{2}_{\rm min}$ with $Q^{2}_{\rm min}=0.268\ {\rm GeV}^{2}$. The obtained results are given in Table 5, for fits with various $Q^{2}_{\rm min}=0.66$, $047$ and $0.268\ {\rm GeV}^{2}$.

Comparing the obtained results in Figs. 7 and Table 5 with the corresponding results in Figs. 4 and Table 2 (where the elastic part was not included in the fit procedure), we see that the inclusion of the elastic contribution in the fit procedure significantly deteriorates (increases) the values of the various fit quality parameters $\chi^{2}$, especially when the fit is performed in the larger interval $(Q^{2}_{\rm min}=)0.268\ {\rm GeV}^{2}\leq Q^{2}<3{\rm GeV}^{2}$. In particular, the differences in quality are clearly visible when comparing Fig. 7(b) with Fig. 4(b) at low $Q^{2}$. Only when $Q^{2}_{\rm min}$ is relatively high, $Q^{2}_{\rm min}=0.66\ {\rm GeV}^{2}$, are some of the $\chi^{2}$ parameters comparable in the two cases (but not the extrapolation quality parameters $\chi^{2}_{0.268}$ and $\chi^{2}_{\rm all}$). Further, the inclusion of the elastic contribution in the fit in general does not change significantly the extracted (mostly negative) values of $f_{2}^{p-n}(1{\rm GeV}^{2})$, but increases significantly the (positive) values of $\mu_{6}$.

For the fit analysis with the “massive” truncated OPE, Eq. (16), the results are presented in Table 6 for $Q^{2}_{\rm min}=0.268\ {\rm GeV}^{2}$.

In general, the results are worse than in the corresponding “massive” OPE case fitted to the inelastic contribution (cf. Table 3 in Sec. III.2). Namely, the extracted values of the parameter $M^{2}$ become negative ($M^{2}\sim-0.1\ {\rm GeV}^{2}$) and the fit quality values $\chi^{2}/{\rm d.o.f.}$ are in general significantly higher than those in Table 3. The uncertainties of the extracted values of $f_{2}^{p-n}(1)$ and $M^{2}$ become reduced (in comparison to the case when the elastic part is not included in the fit); this is so because the elastic parts increase significantly BSR $\Gamma_{1}^{p-n}(Q^{2})$ at low $Q^{2}$ where they are represented mostly by the “massive” higher-twist term, and hence the relative uncertainties of BSR at low $Q^{2}$ become smaller.

In addition to the massive case with constant squared mass $M^{2}$, we include also the case of $Q^{2}$-dependent squared mass Eq. (17) in the analysis with the elastic contribution included in the fit. The results for this case are presented in Table 7. These results are qualitatively similar to those with constant squared mass, Table 6; and the comparison of Table 7 with its counterpart Table 4 without the elastic contribution is similar to the above comparison when the squared mass is constant. We recall that an average mass in the $Q^{2}$-independent case can be regarded to be $\langle M^{2}(Q^{2})\rangle=(1/2)\times(M^{2}(0.268)+M^{2}(1))\approx 1.85\times M^{2}(1\ {\rm GeV}^{2})$.

IV Summary

Experimental results for the polarized Bjorken sum rule (BSR) $\Gamma_{1}^{p-n}(Q^{2})$ were fitted, for various ranges of $Q^{2}$, with OPE theoretical expressions using QCD couplings obtained in four different approaches: perturbative QCD (pQCD) in $\overline{\rm MS}$ scheme; (Fractional) Analytic Perturbation Theory [(F)APT]; Two-delta ${\mathcal{A}}$QCD (2$\delta$); and Three-delta lattice-motivated ${\mathcal{A}}$QCD (3$\delta$). The QCD running coupling ${\mathcal{A}}(Q^{2})$ in the latter three QCD variants does not have Landau singularities, in contrast to the pQCD coupling $a(Q^{2})$ [$\equiv\alpha_{s}(Q^{2})/\pi$].

In the fit of the inelastic experimental BSR results, up to two higher-twist terms [$\sim 1/Q^{2},1/(Q^{2})^{2}$] were added to the theoretical leading-twist contribution. The elastic contributions, which are $\sim 1/(Q^{2})^{n}$ with typically $n\geq 4$ at $Q^{2}>1\ {\rm GeV}^{2}$, were then added by using the parametrization obtained from the literature. The fits were performed for the ranges $Q^{2}_{\rm min}\leq Q^{2}\leq 3\ {\rm GeV}^{2}$, where $Q^{2}_{\rm min}=0.66$, $0.47$ and $0.268\ {\rm GeV}^{2}$. In general, the best curves were obtained when 2$\delta$ or 3$\delta$-couplings were used. When only $D=2$ ($\sim 1/Q^{2}$) higher-twist term was included in the fit, the quality of the fitted curves, in the range of the fit and in the extrapolated ranges of $Q^{2}$, in general did not depend significantly on $Q^{2}_{\rm min}$ of the fit. On the other hand, when both $D=2$ and $D=4$ terms were included, the quality in the extrapolated ranges of $Q^{2}$ was in general better for the lowest $Q^{2}_{\rm min}$ value ($0.268\ {\rm GeV}^{2}$), i.e., when the $Q^{2}$-range of the fit was the largest. Comparably good results were obtained when “massive” higher-twist term was used in the OPE and the QCD coupling was either from (F)APT or 2$\delta$ or 3$\delta$ ${\mathcal{A}}$QCD.

When the range of fit had $Q^{2}_{\rm min}=0.268\ {\rm GeV}^{2}$, the pQCD $\overline{\rm MS}$ coupling approach worked and gave acceptable results only if the renormalization scale of the coupling was maintained everywhere at sufficiently high values, and the coefficient $f_{2}^{p-n}(Q^{2})$ [$\sim a(Q^{2})^{\gamma_{0}/8\beta_{0}}$] at the $D=2$ term had an (ad hoc) increased scale $Q^{2}\mapsto kQ^{2}$, in order to avoid the problem of the Landau singularities.

When the fit procedure was performed by fitting the theoretical OPE, truncated at $D=4$ ($\sim 1/(Q^{2})^{2}$) terms, to the sum of (experimental) inelastic and (parametrized) elastic BSR, the quality of the results turned out to be significantly worse in all the cases of the theoretical curves, something expected by the arguments presented at the end of Sec. II.3. Namely, the elastic contribution is dominated by terms which behave at high values of $Q^{2}$ as $\sim 1/(Q^{2})^{(D/2)}$ where usually $D\geq 8$, and these terms are not contained in the theoretical expressions for BSR which are usually OPE series truncated at $1/(Q^{2})^{2}$.

The results of this work can be interpreted as an additional indication of the following important property: the evaluation of the (truncated) leading-twist contribution of spacelike low-$Q^{2}$ QCD observables such as inelastic BSR, in QCD variants 2$\delta$ and in particular 3$\delta$ ${\mathcal{A}}$QCD [both have infrared finite and holomorphic coupling ${\mathcal{A}}(Q^{2})$], appear to resum effectively a large part of the perturbative contribution of the observables, and leads to reduced extracted values of the higher-twist terms ($D=2,4$) in the truncated OPE. This property was noted earlier, for different observables, in Refs. anOPE ; 3l3danQCD ; 4l3danQCD . In this context, it appears to be important that in 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD the coupling practically merges with the underlying pQCD coupling $a(Q^{2})$ at higher values of $Q^{2}\gg\Lambda^{2}_{\rm QCD}$. This property is not shared by the (F)APT holomorphic coupling where the leading-twist series contains parts of the higher-twist contribution of as low dimensionality as $D=2$. The extracted parameters in the higher-twist contribution, including those in the “massive” OPE, are especially reduced in 3$\delta$ ${\mathcal{A}}$QCD. This suggests the possibility that the true higher-twist contribution is small, including the (sum of) terms of high dimension; and that the (truncated) OPE with 3$\delta$ ${\mathcal{A}}$QCD leading-twist gives, through fitting, an extracted value which is a good approximation to this true value of the higher-twist contribution. Numerically, the significantly reduced extracted value in (truncated) OPE with 3$\delta$ ${\mathcal{A}}$QCD is probably partly related with the fact that 3$\delta$ ${\mathcal{A}}$QCD differs from both 2$\delta$ and (F)APT ${\mathcal{A}}$QCD variants in that its coupling goes to zero in the deep infrared regime, ${\mathcal{A}}^{(3\delta)}(Q^{2})\sim Q^{2}\to 0$. The latter property, we recall, is suggested by the large-volume lattice calculations of the dressing functions of the Landau-gauge gluon and ghost propagators at low $Q^{2}$ values.

Acknowledgements.

This work was supported by FONDECYT Postdoctoral Grant No. 3170116 (C.A.), by FONDECYT Regular Grant No. 1180344 (G.C. and C.A.), and by the RFBR Foundation through Grant No. 16-02-00790-a (A.V.K. and B.G.S.). We thank A. Deur, J. Blümlein and A. L. Kataev for helpful suggestions.

Appendix A Renormalization scale and scheme dependence of the expansion coefficients

The dependence of the coupling $a(\mu^{2})\equiv\alpha_{s}(\mu^{2})/\pi$ (where $\mu^{2}>0$ means the spacelike region) is governed by the perturbative renormalization group equation (RGE)

[TABLE]

where in the mass-independent schemes the coefficients $\beta_{0}=(11-2N_{f}/3)/4$ and $\beta_{1}=(102-38N_{f}/3)/16$ are universal (scheme-independent), while the coefficients $\beta_{j}$ (or equivalently, $c_{j}\equiv\beta_{j}/\beta_{0}$) for $j\geq 2$ are the (arbitrary) parameters which characterize the renormalization scheme.111111There is another scheme parameter, the scale $\Lambda^{2}$, such that $a(\mu^{2})=f(\mu^{2}/\Lambda^{2})$. However, the change of $\Lambda^{2}$ can be regarded as the change in the definition of the renormalization scale $\mu^{2}$. The running of the coupling $a(\mu^{2};c_{2},c_{3},\ldots)$ with these scheme parameters is governed by the following relations (cf. App. A of Ref. Stevenson , and App. A of Ref. GCRK63 ):

[TABLE]

When we use the relations (20)-(21) in the perturbation expansion (4) and account for the fact that ${\cal D}_{\rm BSR}(Q^{2})$ is a (spacelike) observable and thus independent of the scale $\mu^{2}$ (i.e., independent of $k\equiv\mu^{2}/Q^{2}$) and of the scheme parameters $c_{j}$ ($j\geq 2$), we obtain the following expressions for the perturbation coefficients $d_{j}$ in terms of the general renormalization scale and scheme parameters ($k,c_{2},c_{3}$):

[TABLE]

Here we used the bar symbol to denote the choice of the scheme $\overline{\rm MS}$ with the renormalization scale $\mu^{2}=Q^{2}$ ($k=1$). We recall that $k\equiv\mu^{2}/Q^{2}$ is the renormalization scale parameter ($0<k\sim 1$), and that ${\bar{c}}_{j}\equiv{\bar{\beta}}_{j}/\beta_{0}$ ($j\geq 2$) are the $\overline{\rm MS}$ scheme parameters.

Appendix B Power analogs ${\mathcal{A}}_{n}$ in ${\mathcal{A}}$QCD

A QCD running coupling ${\mathcal{A}}(Q^{2})$ which is holomorphic (analytic) in the non-timelike $Q^{2}$ complex plane sector $(-q^{2}\equiv)Q^{2}\in\mathbb{C}\backslash(-\infty,-M_{\rm thr}^{2}]$ (where $0\leq M_{\rm thr}^{2}\lesssim 1\ {\rm GeV}^{2}$), has in general nonperturbative (NP) contributions $\sim 1/(Q^{2})^{n}$ appreciable at small $|Q^{2}|$. It differs from the underlying pQCD coupling $a(Q^{2})$ by

[TABLE]

for $|Q^{2}|>\Lambda^{2}$ ($\gtrsim 0.1\ {\rm GeV}^{2}$). Here, $N=1$ in the case of (F)APT, and $N=5$ in 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD. The analytization procedure can be presented schematically as $a(Q^{2})\mapsto{\mathcal{A}}(Q^{2})$. This procedure involves $a$ and ${\mathcal{A}}$ linearly (not as powers). Namely, when $Q^{2}$ is varied, $Q^{2}\mapsto Q^{2}+\Delta Q^{2}$, we have $a(Q^{2}+\Delta Q^{2})\mapsto{\mathcal{A}}(Q^{2}+\Delta Q^{2})$. Therefore,

[TABLE]

where we denoted the (logarithmic) derivatives

[TABLE]

In this notation, ${{\widetilde{a}}}_{1}=a$ and ${{\widetilde{\mathcal{A}}}}_{1}\equiv{\mathcal{A}}$. We note that by pQCD RGE (20) we have

[TABLE]

More specifically, we have

[TABLE]

where, as mentioned in Appendix A, $c_{j}\equiv\beta_{j}/\beta_{0}$. When we invert these relations, we obtain

[TABLE]

The linearity of analytization, Eq. (24), then gives us the analogs ${\mathcal{A}}_{n}$ of the powers $a^{n}$

[TABLE]

We note that in general ${\mathcal{A}}_{n}(Q^{2})\not={\mathcal{A}}(Q^{2})^{n}$. The described construction (for $n=1,2,3,\ldots$) was performed in CV1 ; CV2 .

The above approach was extended in Ref. GCAK to the case of general real index $n=\nu$

[TABLE]

where ${\rm Li}_{-\nu+1}(z)$ is the polylogarithm function of order $-\nu+1$, and $\rho_{{\mathcal{A}}}(\sigma)={\rm Im}{\mathcal{A}}(Q^{2}=-\sigma-i\epsilon)$ is the cut discontinuity (spectral) function of ${\mathcal{A}}$. The coupling ${{\widetilde{\mathcal{A}}}}_{\nu}$ can also be presented in an alternative form applicable in an extended region $-1<\nu$ (cf. GCAK for details). The expression ${\mathcal{A}}_{\nu}$, the analog of the power $a^{\nu}$, was then obtained in the form

[TABLE]

with the coefficients ${\widetilde{k}}_{m}(\nu)$ given in Appendix A of Ref. GCAK . Eqs. (29) are a special case of Eq. (31).

The perturbation expansion of the type Eq. (4) in Sec. II.1, for any spacelike obervable ${\cal D}(Q^{2})$ in pQCD, can be reexpressed in terms of derivatives ${{\widetilde{a}}}_{n}$ of Eq. (25a)

[TABLE]

where we denoted $a\equiv a(kQ^{2};c_{2},\ldots)$, ${{\widetilde{a}}}_{n}\equiv{{\widetilde{a}}}_{n}(kQ^{2};c_{2},\ldots)$, and the coefficients ${\widetilde{d}}_{n}\equiv{\widetilde{d}}_{n}(k;c_{2},\ldots,c_{n})$ of this “modified” perturbation expansion (mpt) can be obtained by using the RGE-relations (28)

[TABLE]

and the coefficients $d_{j}\equiv d_{j}(k;c_{2},\ldots,c_{j})$ are those of Eqs. (22). The expressions in ${\mathcal{A}}$QCD, corresponding to the perturbation expansions (32) and (4), are then

[TABLE]

The expansion (34b) is written again, in a more detailed form, in Eq. (5) in Sec. II.1. Both expressions (34) are equivalent, but in practice it is more economical to do numerical evaluations using the expression (34a).

We have ${\mathcal{A}}_{n}(Q^{2})={\mathcal{A}}(Q^{2})^{n}$ only when ${\mathcal{A}}(Q^{2})$ is a perturbative coupling, i.e., when it has no NP terms ($\sim 1/(Q^{2})^{m}$). It is important not to use the power expansion in ${\mathcal{A}}$ for the evaluation of spacelike observables. Namely, if we used power expansion in ${\mathcal{A}}$, the truncated series for ${\cal D}(Q^{2})$ would have increasingly large (out of control) NP contributions when more power terms were included, and renormalization scale invariance would be increasingly violated, as emphasized in Techn . It turns out that in practice the sequence ${{\widetilde{\mathcal{A}}}}_{n}(Q^{2})$ ($n=0,1,2,\ldots$) is, in a general holomorphic ${\mathcal{A}}$QCD, a sequence with decreasing absolute values, at any finite $Q^{2}$: $|{{\widetilde{\mathcal{A}}}}_{n}(Q^{2})|>|{{\widetilde{\mathcal{A}}}}_{n+1}(Q^{2})|>\ldots$. In pQCD (${{\widetilde{a}}}_{n}(Q^{2})$) this is in general not valid at low values $|Q^{2}|\lesssim 1\ {\rm GeV}^{2}$.

Appendix C QCD variants with holomorphic coupling ${\mathcal{A}}(Q^{2})$

C.1 (Fractional) Analytic Perturbation Theory [(F)APT]

The pQCD running coupling $a(Q^{2})$, in a given renormalization scheme (usually $\overline{\rm MS}$), has in the complex $Q^{2}$-plane cut along the real axis, $(-\infty,\Lambda_{\rm Lan.}^{2})$, where $0<\Lambda_{\rm Lan.}^{2}\sim 0.1\ {\rm GeV}^{2}$ is the branching point of the interval of the Landau singularities $(0,\Lambda_{\rm Lan.}^{2})$ in the plane. Spacelike QCD observables ${\cal D}(Q^{2})$ are holomorphic (analytic) functions of complex $Q^{2}$, with the exception of the negative (timelike) semiaxis $(-\infty,-M^{2}_{\rm thr.})$ where $0\leq M^{2}_{\rm thr.}\lesssim 1\ {\rm GeV}^{2}$ is a threshold scale. The pQCD coupling $a(Q^{2})$ does not reflect these properties, because of the mentioned cut interval $(0,\Lambda_{\rm Lan.}^{2})$, called Landau singularities, on the positive semiaxis. Application of the Cauchy theorem to the integrand $a(Q^{{}^{\prime}2})/(Q^{{}^{\prime}2}-Q^{2})$ in the complex $Q^{{}^{\prime}2}$-plane, with the use of the asymptotic freedom of QCD ($|a(Q^{{}^{\prime}2})|\to 0$ when $|Q^{{}^{\prime}2}|\to\infty$), then gives the following dispersion integral for the value of the pQCD coupling $a(Q^{2})$:

[TABLE]

where $\rho_{1}^{\rm(pt)}(\sigma)={\rm Im}a(Q^{{}^{\prime}2}=-\sigma-i\epsilon)$ is the cut discontinuity (spectral) function of $a(Q^{{}^{\prime}2})$. The elimination of the Landau cut contribution in this integral, while keeping the spectral function unchanged at other $\sigma$, then gives us the APT coupling ShS

[TABLE]

which has $M^{2}_{\rm thr.}=0$. It is straightforward to check that the difference ${\mathcal{A}}^{\rm(APT)}(Q^{2})-a(Q^{2})$ at large $|Q^{2}|>\Lambda^{2}_{\rm Lan.}$ remains appreciable, $\sim(\Lambda_{\rm Lan.}^{2}/Q^{2})$, i.e., the index $N$ in eq. (23) is $N=1$.

The analog ${\mathcal{A}}_{\nu}^{\rm(APT)}(Q^{2})$ of the power $a(Q^{2})^{\nu}$ is then constructed in complete analogy, by replacing $\rho_{1}^{\rm(pt)}(\sigma)$ by $\rho_{\nu}^{\rm(pt)}(\sigma)={\rm Im}a(-\sigma-i\epsilon)^{\nu}$ MS96 ; Sh

[TABLE]

We use in this work this form of (F)APT couplings ${\mathcal{A}}_{\nu}$ for $\nu=1,2,\ldots$, in $\overline{\rm MS}$ scheme. Specifically, we apply the underlying $\overline{\rm MS}$ pQCD coupling $a(Q^{2};\overline{\rm MS})$ with $N_{f}=3$ to evaluate $\rho_{n}^{\rm(pt)}$ and thus ${\mathcal{A}}_{n}^{\rm(APT)}(Q^{2})$.

The authors of BMS05 ; BKS05 obtained explicit form of ${\mathcal{A}}_{\nu}^{\rm(APT)}(Q^{2})$ at the one-loop level of the underlying pQCD coupling, and extended it to higher loop level by a perturbative approach BMS06 ; BMS10 . This theory has the name Fractional Analytic Perturbation Theory (FAPT).

Numerical programs were constructed to calculate $A_{\nu}^{\rm(APT)}(Q^{2})$ up to four-loop level of the underlying $\overline{\rm MS}$ pQCD coupling $a$, in NeSi ; mathprg1 in Maple and/or Fortran, and in BK ; mathprg2 in Mathematica.

We note that the general approach presented in Appendix B for calculation of ${\mathcal{A}}_{\nu}$ gives in the APT case (where $\rho_{{\mathcal{A}}}(\sigma)=\rho_{1}^{\rm(pt)}(\sigma)$, for $\sigma>0$) approximately the same numerical results as the approach of Eq. (37). If in the approach described in Appendix B we take into account in Eq. (31) [or (29) when $\nu=n$ is integer] a large number of terms ${{\widetilde{\mathcal{A}}}}_{\nu+m}$,121212This is in principle not necessary; if, for example, a physical quantity ${\cal D}(Q^{2})$ is calculated up to ${\cal O}(a^{4})\mapsto{\cal O}({\mathcal{A}}_{4})$, then only terms up to ${\widetilde{\mathcal{A}}}_{4}$ are in principle needed on the right-hand sides of Eqs. (29). it turns out that the obtained ${\mathcal{A}}_{\nu}(Q^{2})$ numerically converges to that in Eq. (37). We wish to stress that the approach presented in Appendix B works for any ${\mathcal{A}}$QCD, i.e., QCD with any holomorphic ${\mathcal{A}}(Q^{2})$, while the approach Eq. (37) is applicable only in the (F)APT case, i.e., when $\rho_{{\mathcal{A}}}(\sigma)=\rho_{1}^{\rm(pt)}(\sigma)$.

The only adjustable parameter in (F)APT (with $N_{f}=3$) is the scale ${\overline{\Lambda}}_{N_{f}=3}$ of the underlying QCD coupling $a(Q^{2})=f(Q^{2}/{\overline{\Lambda}}_{3}^{2})$; the high energy QCD is approximately reproduced with $N_{f}=3$ (F)APT when ${\overline{\Lambda}}_{3}=0.45$ GeV (cf. also Ref. Sh ); we use this value in our analysis, but we also comment on the case ${\overline{\Lambda}}_{3}=0.40$ GeV.

C.2 2$\delta$ and 3$\delta$ ${\mathcal{A}}$QCD

This type of QCD variants 2danQCD ; 3l3danQCD ; 4l3danQCD with coupling ${\mathcal{A}}(Q^{2})$ holomorphic in $Q^{2}\in\mathbb{C}\backslash(-\infty,-M_{\rm thr}^{2}]$ are constructed on the idea that: (a) the coupling ${\mathcal{A}}(Q^{2})$ at high $|Q^{2}|>1\ {\rm GeV}^{2}$ should practically coincide with the underlying pQCD131313This means that in Eq. (23) the index $N$ is large; as a consequence, the well demostrated success of pQCD at high momenta is reproduced in ${\mathcal{A}}$QCD. coupling $a(Q^{2})$; (b) and at moderate and low $|Q^{2}|\lesssim 1{\rm GeV}^{2}$ the coupling should reproduce the well measured semihadronic $\tau$-lepton decay physics and possibly some other experimental indicators. We note that (F)APT does not fulfill these requirements.

The condition (a) then implies that the spectral function $\rho_{{\mathcal{A}}}(\sigma)\equiv{\rm Im}{\mathcal{A}}(-\sigma-i\epsilon)$ is at large $\sigma>1\ {\rm GeV}^{2}$ (approximately) equal to the spectral function of the underlying pQCD, $\rho_{a}(\sigma)\equiv\rho_{1}^{\rm(pt)}(\sigma)\equiv{\rm Im}a(-\sigma-i\epsilon)$. At low positive $\sigma$, the unknown behavior of the spectral function $\rho_{{\mathcal{A}}}(\sigma)$ is parametrized as a sum of delta functions

[TABLE]

Implicitly, we expect $M_{1}^{2}<M_{2}^{2}<\cdots<M_{n}^{2}<M_{0}^{2}$, where $M_{1}^{2}=M^{2}_{\rm thr.}$ is the mentioned threshold scale (expected to be $\sim m_{\pi}^{2}\sim 10^{-2}\ {\rm GeV}^{2}$), and $M_{0}^{2}$ ($\sim 1\ {\rm GeV}^{2}$) can be called the pQCD-onset scale. Application of the Cauchy theorem then gives for the coupling the expression

[TABLE]

The parametrization with $n$ delta functions in Eq. (38) means that the part $\Delta{\mathcal{A}}$ of the QCD coupling originating from the unknown low-$\sigma$ part of the spectrum is parametrized by a near-diagonal Padé $[n-1/n](Q^{2})$ approximant

[TABLE]

where $P_{n-1}$ and $Q_{n}$ are polynomials of degree $n-1$ and $n$ in $Q^{2}$, respectively. If $\Delta{\mathcal{A}}$ is a Stieltjes function141414If $\rho_{{\mathcal{A}}}(\sigma)\geq 0$ for $\sigma>0$, then ${\mathcal{A}}(Q^{2})$ and $\Delta{\mathcal{A}}(Q^{2})$ are Stieltjes functions., a theorem in the Padé theory BakerMorris (cf. also Peris ) states that there exists a sequence of Padé approximants $[n-1/n](Q^{2})$ which converges to $\Delta{\mathcal{A}}(Q^{2})$ for any $Q^{2}\in\mathbb{C}\backslash(-\infty,-M_{\rm thr}^{2}]$ when $n$ increases (if $\Delta{\mathcal{A}}$ is not Stieltjes, it is not known whether such a convergence is guaranteed).

The underlying pQCD coupling $a(Q^{2})$ is determined, to a rather high degree of accuracy, by the world average value $a(Q^{2}=M_{Z}^{2};\overline{\rm MS})=0.1185/\pi$ PDG2014 . We use this value, and we RGE-evolve $a(\overline{\rm MS})$ to lower values of $|Q^{2}|$, by using the four-loop RGE (20) and three-loop quark mass threshold conditions CKS in $\overline{\rm MS}$, into the region of interest where $N_{f}=3$; subsequently, we change the coupling $a(Q^{2};\overline{\rm MS})$ to the considered renormalization scheme to obtain the underlying coupling $a(Q^{2})$; cf. 4l3danQCD for more details. Thus we also obtain the perturbative spectral function $\rho_{1}^{\rm(pt)}(\sigma)={\rm Im}a(Q^{2}=-\sigma-i\epsilon)$, in the considered scheme and for $N_{f}=3$.

At that point, the considered ${\mathcal{A}}$QCD, Eqs. (38)-(39) has altogether $(2n+1)$ parameters, namely ${\cal F}_{j}$ and $M^{2}_{j}$ ($j=1,\ldots,n$) and the pQCD-onset scale $M_{0}^{2}$. These are to be fixed by the conditions (a) and (b) for ${\mathcal{A}}(Q^{2})$ at high and low $|Q^{2}|$ mentioned at the beginning of this Section C.2. The condition (a) is implemented by requiring a large index value for the difference Eq. (23); in our considered cases of 2$\delta$ and 3$\delta$ coupling we took $N=5$

[TABLE]

This represents altogether four conditions (for each increase of $N$ by one, from $N=1$, there is one condition).

When we take $n=2$, i.e., two delta functions in the spectral function (38), we have five parameters to determine; therefore one additional condition is needed. This condition will be the reproduction of the measured values of the quantity $r^{(D=0)}_{\tau}$, the semihadronic strangeless $\tau$ decay rate ratio (the leading-twist part, and with mass effects subtracted).151515Higher-twist contributions in $r_{\tau}$ decay ratio are known to be strongly suppressed. Its experimental value is approximately in the range $0.201\pm 0.002$. Its theoretical expression can be represented as a weighted countour integral of the (massless) Adler function161616This is the canonical Adler function with $N_{f}=3$, i.e., it is normalized in such a way that $d(Q^{2})_{\rm pt}=a(Q^{2})+{\cal O}(a^{2})$. $d(Q^{2})$

[TABLE]

We refer to 2danQCD ; 3l3danQCD ; 4l3danQCD for details. In the $n=2$ case (2$\delta$ ${\mathcal{A}}$QCD), we still have the freedom of choosing the renormalization scheme. We took it as the Lambert scheme ($c_{j}=c_{2}^{j-1}/c_{1}^{j-2}$, for $j\geq 3$), with $c_{2}=-4.9$ mathprg2 . It is possible to vary the value of $c_{2}$, but when it is different by several units from this value, either the pQCD-onset scale $M_{0}$ becomes appreciably higher than $\approx 1$ GeV, or the value ${\mathcal{A}}(0)$ becomes larger than one, cf. Table 2 of mathprg2 .

When we take $n=3$, i.e., three delta functions in the spectral function (38), there are two additional parameters to be fixed. These two parameters are fixed by the condition ${\mathcal{A}}(Q^{2})\sim Q^{2}$ when $Q^{2}\to 0$ and the local maximum of ${\mathcal{A}}(Q^{2})$ achieved at $Q^{2}\approx 0.135\ {\rm GeV}^{2}$. These conditions are motivated by the results of lattice calculations of the gluon and ghost dressing functions in the Landau gauge at low positive $Q^{2}$ LattcoupNf0 (cf. also LattcoupNf0b ).171717The calculation in LattcoupNf0 ; LattcoupNf0b were performed for $N_{f}=0$ case. They are similar to the results when $N_{f}=2,4$ LattcoupNf2 ; LattcoupNf4 although these results are not so precise. If the running is defined via the mentioned dressing functions, the mentioned conditions, with ${\mathcal{A}}(0)=0$, follow. Different definitions involving, in addition, a dynamical gluon mass DSEdecFreez ; PTBMF , would imply ${\mathcal{A}}(0)>0$, which holds also in 2$\delta$ ${\mathcal{A}}$QCD and in (F)APT, with the values ${\mathcal{A}}(0)=2.0713,1.3963$, respectively (for $N_{f}=3$). For a discussion on these issues, we refer to 4l3danQCD . In our couplings we use throughout $N_{f}=3$, which makes them applicable in the regions $|Q^{2}|<3\ {\rm GeV}^{2}$. The lattice calculations were performed in the MiniMOM renormalization scheme (MM) MiniMOM ; BoucaudMM ; CheRet ; KatMol . Our coupling ${\mathcal{A}}(Q^{2})$ was constructed in the $N_{f}=3$ MM scheme, but rescaled from the MM scale convention ($\Lambda_{\rm MM}$) to the usual scale convention ($\Lambda_{\overline{\rm MS}}$), the latter representing what we call Lambert MM (LMM) scheme. In 3l3danQCD we used the three-loop LMM, and in 4l3danQCD the four-loop LMM scheme. In the present work we use the latter (four-loop LMM) scheme, i.e., the coupling from 4l3danQCD .

In Table 8 we specify the parameters of the 2$\delta$ ${\mathcal{A}}$QCD 2danQCD ; mathprg2 and 3$\delta$ ${\mathcal{A}}$QCD coupling 4l3danQCD , used in the present work.

The values of the parameters of the 2$\delta$ coupling, with $c_{2}=-4.9$, are slightly different from the corresponding values in mathprg2 (Table 2 there, third line), because there we used $\alpha_{s}(M_{Z}^{2};\overline{\rm MS})=0.1184$ and $r_{\tau}^{(D=0)}=0.203$, while here we use $\alpha_{s}(M_{Z}^{2};\overline{\rm MS})=0.1185$ and $r_{\tau}^{(D=0)}=0.201$. We note that the maximal value of ${\mathcal{A}}(Q^{2})$ for positive $Q^{2}$ is achieved in the 2$\delta$ case at $Q^{2}=0$: $\pi{\mathcal{A}}^{(2\delta)}(0)=2.0713$; and in the 3$\delta$ case at $Q^{2}=Q_{0}^{2}=0.1348\ {\rm GeV}^{2}$: $\pi{\mathcal{A}}^{(3\delta)}(Q_{0}^{2})=0.9156$. The pQCD-onset scales are $M_{0}=\sqrt{s_{0}}\Lambda_{{\rm L.}}=1.298$ GeV and $2.953$ GeV, respectively.

Other QCD running couplings ${\mathcal{A}}(Q^{2})$ without Landau singularities have been constructed in the literature, some of them having at $Q^{2}=0$ zero value ${\mathcal{A}}(0)=0$ ArbZaits ; Boucaud ; mes2 , finite nonzero value (for reviews, cf. BrazJP ; Brodrev ), or infinite value Nest1 . The construction of such couplings is mostly based on the dispersive approaches similar to the ones described in this Appendix; however, such kind of dispersive approaches can also be applied to entire physical observables, cf. ShS98 ; MSS2 ; MagrGl ; mes2 ; DeRafael ; MagrTau ; Nest3a ; Nest3b ; NestBook . Yet another approach is the Light Front Holographic (LFH) QCD LFH ; LFHBSR , where the coupling ${\mathcal{A}}^{(\rm(LFH})(Q^{2})$ has the form $\exp(-Q^{2}/2\kappa^{2})$ at low positive $Q^{2}$, cf. Eq. (19), and can be extended to higher positive $Q^{2}$ by matching of ${\mathcal{A}}(Q^{2})$ and $d{\mathcal{A}}(Q^{2})/dQ^{2}$ at a matching scale $Q_{0}\sim 1\ {\rm GeV}^{2}$ to pQCD LFHmatch .

Appendix D Uncertainties of extracted values of $f^{p-n}_{2}$ and $\mu_{6}$: statistical and systematic

The experimental data in the fit procedure are $\Gamma_{1}^{p-n}(Q_{j}^{2};{\rm inel.})_{\rm exp.}$. The leading-twist part of the (theoretical) OPE expression (10) has no fit parameters in our approach, apart from the renormalization scale parameter $k\equiv\mu^{2}/Q^{2}$. The latter parameter will be kept fixed at its value obtained by the approach of minimization of $\chi^{2}$. We will use the following notation:

[TABLE]

In order to estimate the uncertainties of the extracted parameters $f^{p-n}_{2}(1)$ and $\mu_{6}$, we will consider for simplicity that not just $\mu_{6}$ but also $f_{2}^{p-n}$ does not run (i.e., that $\mu_{4}$ does not run) with $Q^{2}$. Further, we will make the following approximations: the statistical uncertainties at different points will be considered as completely uncorrelated; the systematic uncertainties at different points will be considered as completely correlated within one experiment, and completely uncorrelated between two different experiments. The systematic uncertainties $\delta\Gamma_{1}^{p-n}(Q_{j}^{2};{\rm inel.})_{\rm sys}$ are in general significantly larger than the statistical uncertainties $\delta\Gamma_{1}^{p-n}(Q_{j}^{2};{\rm inel.})_{\rm stat}$, cf. Fig. 8.

D.1 $\mu_{6}=0$ case

First we will consider the case when $\mu_{6}=0$ (and $f_{2}^{p-n}$ is not running with $Q^{2}$). We will estimate the uncertainty of $f_{2}^{p-n}$ due to the statistical uncertainties $\sigma_{j,{\rm stat}}$ at $Q_{j}^{2}$’s and then due to the systematic uncertainties $\sigma_{j,{\rm sys}}$. The renormalization scheme parameter $k$ will be kept fixed during these estimations.

In this approach, the considered $\chi^{2}$ is

[TABLE]

where we will denote

[TABLE]

i.e., the squared statistical uncertainty at the measured point $Q^{2}_{j}$. The minimization of $\chi^{2}$ then gives the extracted value of $\mu_{4}$

[TABLE]

where we will use throughout the notation for the unnormalized “average”

[TABLE]

Since the deviations of $y_{j}\equiv y(z_{j})$ due to statistical uncertainties are considered independent at different $z_{j}$, the rule of the sum of the squares of standard deviations is valid. It is then straightforward to deduce the square of the standard deviation for ${\hat{\mu}}_{4}$.

[TABLE]

Then it can be checked, by Taylor expansion of $\chi^{2}(\mu_{4})$ around the point ${\hat{\mu}}_{4}$ up to the terms $(\delta\mu_{4})^{2}$, that the following (approximate) relation holds:

[TABLE]

Any of the two relations (48)-(49) can be used to evaluate the uncertainty $(\delta{\hat{\mu}}_{4})_{stat}\equiv\sigma({\hat{\mu}}_{4})_{\rm stat}$. We used the relation (49), which gives us in practice somewhat higher values of the uncertainty.181818This difference, of a few percent ($\leq 5\%$), is presumably principally due to the effect of the running of $\mu_{4}(Q^{2})$ with $Q^{2}$, the effect not accounted for in the formula (48). Further, using the relation (6) at $Q^{2}=1\ {\rm GeV}^{2}$, this then gives us the uncertainty $(\delta{\hat{f}}_{2})_{\rm stat}$191919We denote, from here on, by ${\hat{f}}_{2}$ the value of $f_{2}^{p-n}(1{\rm GeV}^{2})$ extracted by the fit procedure.

[TABLE]

The systematic uncertainties were estimated in the following way. For simplicity we consider for this only two experiments, namely EG1b EG1b (experiment “1”) and the newer JLAB EG1-DVCS results DVCS (experiment”2”), in the interval $0<Q^{2}<3\ {\rm GeV}^{2}$. These two experiments represent most of the available data points, and each of them covers most of the mentioned $Q^{2}$-interval. In each experiment, the systematic uncertainties at different points will be considered as strongly correlated, in the sense that in each experiment we will estimate the systematic uncertainty of ${\hat{f}}_{2}$ as the one obtained by averaging the deviations

[TABLE]

where $E$ is the experiment ($E=1$ or $E=2$); “UP” refers to the value of $f_{2}^{p-n}(1{\rm GeV}^{2})$ extracted from the data of the experiment increased by $\sigma_{j,{\rm sys}}$ at the points $Q_{j}^{2}$ of the experiment, i.e., $\Gamma_{1}^{p-n}(Q_{j}^{2};{\rm inel.})+\sigma_{j,{\rm sys}}$; analogously, “DO” refers to the value extracted from $\Gamma_{1}^{p-n}(Q_{j}^{2};{\rm inel.})-\sigma_{j,{\rm sys}}$; and ${\hat{f}}_{2}^{(E)}$ is the value of $f_{2}^{p-n}(1{\rm GeV}^{2})$ extracted from the central points $\Gamma_{1}^{p-n}(Q_{j}^{2};{\rm inel.})$ of the experiment.

Having these estimates, the question is with what relative weights to combine the two systematic uncertainties $(\delta{\hat{f}}_{2})_{sys}^{(E)}$ for $E=1$ and $E=2$. This will be obtained by combining the unnormalized averages of the expressions appearing in ${\hat{\mu}}_{4}^{(E)}$ of the two experiments, and assuming that the (systematic) deviations for ${\hat{\mu}}_{4}^{(E)}$ of the two experiments are mutually independent. Namely, using the definition of the (unnormalized) averages, Eq. (47), for the two experiments, we have the identities

[TABLE]

Applying these identities to ${\overline{yz}}$ and $\overline{z^{2}}$ in Eq. (46), we obtain

[TABLE]

where

[TABLE]

The assumption (approximation) that the two systematic deviations of ${\hat{\mu}}_{4}$ of experiments 1 and 2 are independent, then leads us to the standard deviation $(\delta{\hat{\mu}}_{4})_{sys}$

[TABLE]

and by relation (6) analogously

[TABLE]

where the estimates $\sigma({\hat{f}}_{2})^{(E)}_{\rm sys}$ ($E=1,2$) are given in Eq. (51). Eq. (56) represents thus an estimate of the uncertainty of the extracted value of $f_{2}^{p-n}(1{\rm GeV}^{2})$ due to systematic uncertainties of the experimental data. In practice, it turns out that this uncertainty is dominated by the results of the experiment 2 (JLAB EG1-DVCS) DVCS , i.e., $(1-\alpha)\approx 1$. This is so because $\overline{z^{2}}^{(2)}\gg\overline{z^{2}}^{(1)}$, since the experiment 2 has significantly larger values of $w_{j}$, i.e., significantly smaller values of $\sigma^{2}_{j,{\rm stat}}$.

D.2 $\mu_{6}\not=0$ case

When the coefficient $\mu_{6}$ is included in the truncated OPE (10) as a fit parameter, the analysis is analogous to the previous Sec. D.1, except that now the algebra is more involved. The values of ${\hat{\mu}}_{4}$ and ${\hat{\mu}}_{6}$ are obtained by simultaneous minimization of

[TABLE]

with respect to $\mu_{4}$ and $\mu_{6}$. This gives

[TABLE]

where

[TABLE]

The corresponding squares of the standard deviations are

[TABLE]

When using Taylor expansion of $\chi^{2}(\mu_{4},\mu_{6})$ around the point $({\hat{\mu}}_{4},{\hat{\mu}}_{6})$ up to the terms quadratic in the deviations, it can be checked that the following (approximate) relations hold:

[TABLE]

We determined the values $\sigma({\hat{\mu}}_{4})_{\rm stat}$ and $\sigma({\hat{\mu}}_{6})_{\rm stat}$ from Eqs. (61).202020Eq. (61b) and Eqs. (60) give the same result for $\sigma({\hat{\mu}}_{6})_{\rm stat}$. For $\sigma({\hat{\mu}}_{4})_{\rm stat}$ [and $\sigma({\hat{f}}_{2})_{\rm stat}$, cf. Eq. (50)], the result of Eq. (61a) differs from that of Eqs. (60) by a few percent ($\leq 5\%$), principally due the the effect of the running of $f_{2}^{p-n}(Q^{2})$ with $Q^{2}$.

Also the systematic uncertainties were estimated analogously to the case of $\mu_{6}=0$, cf. Sec. D.1, only the algebra is now more involved. The basis is again the identity (52), but this time for the quantities $\overline{yz^{2}}$ and $\overline{yz}$ which appear in ${\hat{\mu}}_{4}$ and ${\hat{\mu}}_{6}$, Eqs. (58) and can thus be expressed by ${\hat{\mu}}_{4}$ and ${\hat{\mu}}_{6}$

[TABLE]

Using this, we obtain

[TABLE]

where

[TABLE]

As earlier, we denote by ’(1)’ and ’(2)’ the experiments 1 (EG1b, EG1b ) and 2 (EG1-DVCS, DVCS ), and the constants appearing in Eqs. (64) are

[TABLE]

where the $D$-terms are defined as

[TABLE]

where in the last expression on the right-hand side, the unnormalized averages are made over the experiments 1 and 2.212121This is somewhat different from the unnormalized averages $\overline{A}$ appearing in Eqs. (58)-(60) which are over all the experimental points of the fit, i.e., over more than two experiments.

The estimates of the systematic uncertainties for the quantities ${\widetilde{\mu}}_{4}^{(E)}$ and ${\widetilde{\mu}}_{6}^{(E)}$, for experiments $E=1,2$, are evaluated in complete analogy with Eq. (51) of the previous Sec. D.2

[TABLE]

where $N=4$ or $N=6$. For explanations of the notation ’UP’ and ’DO’ we refer to the previous Sec.D.2. In our approximation, we consider that the systematic uncertainties in the experiments 1 and 2 are mutually independent; this, in conjunction with the relations (63), then gives

[TABLE]

The systematic uncertainty of the extracted value $f_{2}^{p-n}(1{\rm GeV}^{2})$ is obtained then from Eq. (68a) via the relation (6)

[TABLE]

D.3 “Massive” OPE case

When the truncated OPE has, instead, a “massive” term Eq. (16), a similar estimation of the statistical and systematic uncertainties of the extracted parameters $f_{2}^{p-n}(1{\rm GeV}^{2})$ and $M^{2}$ can be made. This is so because, when expanding the “massive” higher-twist term in powers of $1/Q^{2}$, we obtain

[TABLE]

In the first approximation, the “massive” case can thus be considered as the case of truncated OPE (10). To estimate the statistical uncertainties of the extracted $f_{2}^{p-n}(1{\rm GeV}^{2})$ and $M^{2}$, we decided to apply the relations of the type (61)

[TABLE]

To estimate the systematic uncertainties in the “massive” case, the relations (68) for $\mu_{4}$ and $\mu_{6}$ were used. Namely, differentiation of the relation (70) gives an approximate relation

[TABLE]

where in the first two terms on the right-hand side we can use the expressions (68), and in the last term the correlator $\langle\delta{\hat{\mu}}_{4}\delta{\hat{\mu}}_{6}\rangle_{sys}$ can be estimated in a similar way as other correlators described above. Motivated by the relations (63), we chose the following estimate for this correlator:

[TABLE]

We note that in the “massive” case the systematic uncertainties are often quite large, so the estimate (72) should be regarded often as only a rough approximation.

Appendix E Charm mass effects in BSR

The contributions of the finite charm quark mass, i.e., the effects beyond $N_{f}=3$, in the polarized BSR at NLO were evaluated in BFF . When ignoring the $b$-quark contributions (considering $m_{b}\to\infty$), the mentioned effects at NLO can be expressed with the function $C_{\rm pBJ}^{\rm mass.,(2)}(\xi_{c})$, where $\xi_{c}\equiv Q^{2}/m_{c}^{2}$ which is assumed in BFF to be $\xi_{c}\gtrsim 1$, and $m_{c}\approx 1.59$ GeV is the pole mass. This function appears in the NLO coefficient (at $a^{2}$) when the perturbation expansion (4) is reexpressed in terms of the $N_{f}=4$ coupling $a_{N_{f}=4}$

[TABLE]

where $N_{f}=4$ and

[TABLE]

When $Q^{2}\gg m_{c}^{2}$ ($\xi\gg 1$), this function approaches unity quite slowly

[TABLE]

In this limit ($\xi_{c}\to\infty$) this then gives in the expansion (74) the massless expression for the NLO coefficient $d_{1}(N_{f})$ with $N_{f}=4$

[TABLE]

As noted, the convergence to the $N_{f}=4$ case (four massless quarks) is rather slow in BSR. Specifically, $3+C_{\rm pBJ}^{\rm mass.,(2)}(Q^{2}/m_{c}^{2})\approx 3.13$, $3.36$, $3.73$, $3.83$ for $Q^{2}=5$, $10$, $50$, $100\ {\rm GeV}^{2}$, respectively. This indicates that at $Q^{2}=5\ {\rm GeV}^{2}$ (the highest considered experimental $Q^{2}$ in BSR) we are still rather far away from the $N_{f}=4$ case. Therefore, it appears reasonable to use $N_{f}=3$ (i.e., decoupled charm quark) in the polarized BSR for the interval $0<Q^{2}<5\ {\rm GeV}^{2}$ considered in the present work.