Towards tests of quark-hadron duality with functional analysis and spectral function data

TL;DR

This paper develops and tests a functional analysis method to detect quark-hadron duality violations in spectral function data, demonstrating its effectiveness with simulated datasets inspired by experimental measurements.

Contribution

It introduces a new analytic expression for the minimal functional distance and applies it to spectral data, enabling statistically significant detection of duality violations.

Findings

Functional analysis tests can detect duality violations in spectral data.

Excluding high-error regions improves detection sensitivity.

The method is validated with realistic toy model pseudodata.

Abstract

The presence of terms that violate quark-hadron duality in the expansion of QCD Green's functions is a generally accepted fact. Recently, a new approach was proposed for the study of duality violations (DVs), which exploits the existence of a rigorous lower bound on the functional distance, measured in a certain norm, between a "true" correlator and its approximant calculated theoretically along a contour in the complex energy plane. In the present paper we pursue the investigation of functional-analysis based tests towards their application to real spectral function data. We derive a closed analytic expression for the minimal functional distance based on the general weighted $L^2$ norm and discuss its relation with the distance measured in $L^\infty$ norm. Using fake data sets obtained from a realistic toy model in which we allow for covariances inspired from the publicly available ALEPH spectral functions, we obtain by Monte Carlo simulations the statistical distribution of the strength parameter that measures the magnitude of the DV term added to the usual operator product expansion (OPE). The results show that, if the region with large errors near the end-point of the spectrum in $\tau$ decays is excluded, the functional-analysis based tests using either $L^2$ or $L^\infty$ norms are able to detect, in a statistically significant way, the presence of DVs in realistic spectral function pseudodata.

Tables

Table 1. Table 1 : Values of μ 𝜇 \mu from the minimization of δ 2 subscript 𝛿 2 \delta_{2} as a function of the number N 𝑁 N of coefficients included in the truncated sum ( 21 ) for s 0 = 2.76 ​ GeV 2 subscript 𝑠 0 2.76 superscript GeV 2 s_{0}=2.76\,\text{GeV}^{2} . The numerical integrals needed for the c n subscript 𝑐 𝑛 c_{n} are computed as in Eq. ( 48 ). The central μ 𝜇 \mu value is obtained from the median of the distribution.

$N$	$\mu$
5	$1.1\pm 1.4$
15	$1.00\pm 0.21$
100	$1.00\pm 0.22$
150	$1.00\pm 0.23$
300	$1.00\pm 0.23$
500	$1.00\pm 0.23$

Table 2. Table 2 : Values of μ 𝜇 \mu from the minimization of δ 2 subscript 𝛿 2 \delta_{2} and δ 0 subscript 𝛿 0 \delta_{0} for different values of s 0 subscript 𝑠 0 s_{0} . N = 300 𝑁 300 N=300 terms are included in the sum ( 21 ) and in the Hankel matrix used in ( 18 ). The numerical integrals needed for the c n subscript 𝑐 𝑛 c_{n} are computed as in Eq. ( 48 ). The central μ 𝜇 \mu value is obtained from the median of the distribution.

$s_0$	$\mu$ from ${\delta }_2$	$\mu$ from ${\delta }_0$
$2.76$ GeV2	$1.00\pm 0.23$	$1.00\pm 0.26$
$2.84$ GeV2	$1.00\pm 0.29$	$1.01\pm 0.30$
$3.00$ GeV2	$1.01\pm 0.56$	$1.01\pm 0.70$
$m_{\tau }^2$	$1.0\pm 1.6$	$1.0\pm 3.0$

Table 3. Table 3 : Optimal values of μ 𝜇 \mu from the minimization of the exact analytic expression of δ 2 , w subscript 𝛿 2 𝑤 \delta_{2,w} , for two values of s 0 subscript 𝑠 0 s_{0} . Three weights w 𝑤 w are used: w = 1 𝑤 1 w=1 corresponding to the standard L 2 superscript 𝐿 2 L^{2} norm, the weight ( 39 ) with a = 0.96 𝑎 0.96 a=0.96 and the pinched weight ( 50 ). The results are obtained with 5,000 toy data sets, the representative μ 𝜇 \mu value being the median of the distributions.

$s_0$	$w=1$	$w=\frac{\sqrt{1-a^2}}{1-as/s_0}$	$w={\left(1-\frac{s}{s_0}\right)}^2$
2.76  GeV2	$1.00\pm 0.23$	$1.00\pm 0.26$	$1.00\pm 0.70$
$m_{\tau }^2$	$1.0\pm 1.6$	$1.0\pm 2.6$	$1.0\pm 2.7$

Equations

$${\Pi }_{\mathrm{QCD}}(s)={\Pi }_{\mathrm{OPE}}(s)+{\Pi }_{\mathrm{DV}}(s),$$

$$\delta =\underset{\theta \in (0,2\pi )}{\sup }|\Pi (s_0{e}^{i\theta })-{\Pi }_{\mathrm{QCD}}(s_0{e}^{i\theta })|,$$

$$\text{mbox}Im\Pi (s+iϵ)=\sigma (s)$$

$${\delta }_0=\underset{\{\Pi \}}{\min }\underset{\theta \in (0,2\pi )}{\sup }|\Pi (s_0{e}^{i\theta })-{\Pi }_{\mathrm{QCD}}(s_0{e}^{i\theta })|,$$

$${\delta }_{\mathrm{exact}}\ge {\delta }_0.$$

$$z\equiv \tilde{z}(s)=\frac{s}{s_0},$$

$$\Vert F{\Vert }_{L^{\infty }}\equiv \underset{\theta \in (0,2\pi )}{\sup }|F(e^{i\theta })|.$$

$$\Vert F{\Vert }_{L^2}\equiv {\left[\frac{1}{2\pi }{\int }_0^{2\pi }|F(e^{i\theta }){|}^{2}d\theta \right]}^{1/2},$$

$$\Vert F{\Vert }_{L_w^2}\equiv {\left[\frac{1}{2\pi }{\int }_0^{2\pi }|w(e^{i\theta })F(e^{i\theta }){|}^{2}d\theta \right]}^{1/2},$$

$${\delta }_0=\underset{g\in H^{\infty }}{\min }\Vert g-h{\Vert }_{L^{\infty }}.$$

$$h(\zeta )=-\frac{\zeta }{\pi }{\int }_{x_{\pi }}^{1+\eta }dx\frac{\sigma (s_{0}x)}{x(x-\zeta -iϵ)}+{\Pi }_{\mathrm{QCD}}(s_0\zeta ),$$

$$\underset{g\in H^{\infty }}{\min }\Vert g-h{\Vert }_{L^{\infty }}=\underset{w,f\in S^2}{\sup }\left|\frac{1}{2\pi }\underset{|\zeta |=1}{\oint }w(\zeta )f(\zeta )h(\zeta )d\zeta \right|,$$

$$\Vert w{\Vert }_{L^2}\le 1,\Vert f{\Vert }_{L^2}\le 1,$$

$$w(z)=\sum _{n=0}^{\infty }{w}_n{z}^n,f(z)=\sum _{m=0}^{\infty }{f}_m{z}^m,$$

$$\sum _{n=0}^{\infty }{w}_n^2\le 1,\sum _{m=0}^{\infty }{f}_m^2\le 1.$$

$${\mathcal{H}}_{nm}=c_{n+m-1},n,m\ge 1,$$

$$c_n=\frac{1}{\pi }{\int }_0^{1}dx{x}^{n-1}\sigma (s_{0}x)+\frac{1}{2\pi }{\int }_0^{2\pi }d\theta e^{in\theta }{\Pi }_{\mathrm{QCD}}(s_0{e}^{i\theta }),n\ge 1.$$

$${\delta }_0=\Vert \mathcal{H}{\Vert }_{L^2}=\Vert \mathcal{H}\Vert ,$$

$${\delta }_2=\underset{\{\Pi \}}{\min }{\left[\frac{1}{2\pi }{\int }_0^{2\pi }d\theta |\Pi (s_0{e}^{i\theta })-{\Pi }_{\mathrm{QCD}}(s_0{e}^{i\theta }){|}^2\right]}^{1/2},$$

$${\delta }_2=\underset{g\in H^2}{\min }\Vert g-h{\Vert }_{L^2},$$

$${\delta }_2={\left[\sum _{n=1}^{\infty }{c}_n^2\right]}^{1/2},$$

$${\delta }_{2,w}=\underset{g\in H^2}{\min }\Vert w(g-h){\Vert }_{L^2},$$

$${\delta }_{2,w}={\left[\sum _{n=1}^{\infty }{d}_n^2\right]}^{1/2},$$

$$d_n=\frac{1}{\pi }\underset{0}{\overset{1}{\int }}{x}^{n-1}w(x)\sigma (s_{0}x)dx+\frac{1}{2\pi }\underset{0}{\overset{2\pi }{\int }}{e}^{in\theta }w(e^{i\theta }){\Pi }_{\mathrm{QCD}}(s_0{e}^{i\theta })d\theta ,$$

$$d_n=d_{n,\sigma }+d_{n,\mathrm{QCD}},$$

$${\delta }_{2,w}^2=\sum _{n=1}^{\infty }{d}_{n,\sigma }^2+2\sum _{n=1}^{\infty }{d}_{n,\sigma }{d}_{n,\mathrm{QCD}}+\sum _{n=1}^{\infty }{d}_{n,\mathrm{QCD}}^2.$$

$$\sum _{n=1}^{\infty }{d}_{n,\sigma }^2=\frac{1}{{\pi }^2}{\int }_0^{1}dx{\int }_0^{1}dy\frac{w(x)w(y)\sigma (s_{0}x)\sigma (s_{0}y)}{1-xy}.$$

$$d_{n,\mathrm{QCD}}=\frac{1}{4\pi }{\int }_0^{2\pi }\left[e^{in\theta }w(e^{i\theta }){\Pi }_{\mathrm{QCD}}(s_0{e}^{i\theta })+e^{-in\theta }{w}^{\ast }(e^{i\theta }){\Pi }_{\mathrm{QCD}}^{\ast }(s_0{e}^{i\theta })\right]d\theta ,$$

$$2\sum _{n=1}^{\infty }{d}_{n,\sigma }{d}_{n,\mathrm{QCD}}=\frac{2}{{\pi }^2}{\int }_0^{1}dxw(x)\sigma (s_{0}x)$$

$$\times {\int }_0^{\pi }d\theta \frac{\mathrm{Re}\left[w(e^{i\theta }){\Pi }_{\mathrm{QCD}}(s_0{e}^{i\theta })\right](\cos \theta -x)-\mathrm{Im}\left[w(e^{i\theta }){\Pi }_{\mathrm{QCD}}(s_0{e}^{i\theta })\right]\sin \theta }{1-2x\cos \theta +x^2}.$$

$$\sum _{n=1}^{\infty }{d}_{n,\mathrm{QCD}}^2$$

Full text

1 Introduction

The presence of additional terms in the QCD Green’s functions, beyond those generated by OPE (understood as perturbation theory and power corrections), is a generally accepted fact, with support both from theory and phenomenology. According to the standard terminology **[1, 2, 3]**, these terms are said to violate quark-hadron duality. We recall that, in its conventional sense, quark-hadron duality assumes that the description in terms of the OPE, valid on the Euclidian axis and at complex energies, can be analytically continued to match with the description in terms of hadrons, which live on the Minkowski axis.

DVs are supposed to arise from contributions of internal lines with soft momenta in the Feynman diagrams, which are not included in the OPE. Their calculation from first principles is, at least at present, impossible. Quantitative understanding must be based on realistic models, whose main features have been tested against experimental data. Two types of specific mechanisms have been suggested, one considering quarks in an instanton field **[1, 2, 3]**, the other based on narrow-resonance saturation in the large-$N_{c}$* limit [3, 4, 5, 6, 7]. Both mechanisms are materialized in exponentially suppressed terms on the spacelike axis, which exhibit oscillations when analytically continued to the timelike axis. More formal arguments in favour of the existence of DVs are provided by ideas of resurgence and the associated trans-series [8]. The assumption that the OPE, expected to be a divergent expansion with increasing large-order coefficients, is actually an asymptotic series, leads also naturally to the presence of additional exponentially suppressed contributions [9]. However, beyond these somewhat general arguments, no detailed dynamical calculation of the additional contributions present in the theoretical expression of the Green’s functions is available.*

The phenomenological extraction of DVs is far from trivial because one must detect terms exponentially suppressed as function of energy, while an infinity of terms logarithmically and power suppressed, i.e. larger in principle, are neglected in the standard truncated expansions of the Green functions. Since these expansions are actually divergent, the magnitude of the neglected terms can be quite substantial. Moreover, as mentioned above, the confrontation between theory and experiments implies an analytic continuation in the complex energy plane, with its known instabilities and pitfalls. Analyticity is usually exploited by means of a Cauchy integral relation along a contour in the complex plane for the QCD polarization amplitude of interest multiplied by a suitable weight. This allows one to build sum-rules that relate the integrated theoretical predictions on the contour to weighted integrals over the spectral function data on the positive Minkowski axis. The weight is chosen such as to enhance or to suppress the contribution of various terms in the theoretical expression of the amplitude. The impact of DVs for practical calculations is therefore sensitive to the weights that are employed and vary depending on the quantity of interest. When extracting QCD parameters, for example, different weight functions have been advocated. In some works, DVs are explicitly taken into account by means of realistic parametrizations **[10, 11, 13, 12, 14, 15]**, which allows for a quantitative control of their contribution, while in others DVs are ignored on the basis of their suppression by the weight functions employed  **[16, 17, 18]**. The reliability of the different approaches is still being questioned **[18, 19]** and, therefore, a better understanding of DVs would certainly contribute to the precision with which QCD parameters can be extracted. This is particularly true for the determination of the coupling $\alpha_{s}$ from the $\tau$-hadronic spectral functions.

In the recent paper **[20]**, a method based on functional analysis was proposed in order to test the presence of DVs in QCD. The method starts from the obvious remark that the “true” polarization amplitude and its approximate theoretical expression are entirely different functions, with different analytic properties, which cannot coincide in the complex energy plane. Moreover, defining a functional distance, measured in a certain norm, between these two functions along a contour in the complex plane, a rigorous nonzero lower bound on this distance can be shown to exist. In particular, for the functional distances defined in $L^{\infty}$ and in $L^{2}$ norms, the lower bound can be calculated by an explicit algorithm involving the QCD approximant in the complex plane and an infinity of Fourier coefficients obtained from the spectral function (“moments”) measured experimentally on a part of the timelike axis.

As argued in Ref. **[20]**, the minimal distance between the true function and its approximant can be used as a tool for detecting the presence of DVs. In particular, from the variation of the minimal distance with respect to a parameter $\mu$ that measures the strength of the duality violating contribution, one can infer the optimal value of this parameter. Formulated in this way, the problem becomes analogous to the search for new physics beyond the standard model (SM) in experiments at very high energies, where one tests for the presence of new physics through a “strength parameter” $\mu$ of the signal, while treating SM as a background. In our case, the “new physics” is represented by DV terms, while OPE is the background representing the “known physics”.

The application of these ideas to a toy model proposed in Ref. **[6]** indicated that the new approach might be useful for detecting the presence of DVs in QCD. The asymptotic expansion of the exact model contains, besides a purely perturbative term and an expansion identified with higher-dimension terms in the OPE, an additional term that can be interpreted as a DV contribution. The minimal functional distance defined in Ref. **[20]**, calculated with the spectral function of the model and a truncated OPE to mimic the physical cases, displayed a sharp minimum at the true value of the strength $\mu$ of the DV term. In particular, the functional distance measured in $L^{\infty}$ norm proved to be more sensitive to the variation of the parameter $\mu$ than the distance measured in $L^{2}$ norm. However, the effect of the experimental uncertainties inherent in the spectral function used as input was only barely touched in Ref. **[20]**. A detailed investigation of this aspect is crucial for assessing the usefulness of the method to detect DVs from real data. In the present paper we address precisely this problem.

We consider the same toy model proposed in Ref. **[6]**, assuming now that the spectral function is measured only in a finite number of bins with uncertainties and correlations similar to those reported in real experiments on hadronic $\tau$ decays. It turns out that a statistical interpretation of the minimal distances defined by functional analysis is difficult to assess theoretically. Therefore, we perform an empirical study based on pseudodata, where fake data on the spectral function are generated in a number of bins, with a multivariate Gaussian distributions with covariances inferred from the ALEPH covariance matrix for the vector ($V$) channel **[16]**. The statistical distribution of the parameter $\mu$ that measures the magnitude of the DV term added to the usual OPE is then derived by Monte Carlo simulations, allowing the extraction of a mean and a standard deviation. The aim is to establish if the method is able to detect, in a statistically significant way, the presence of DVs from error-affected experimental measurements. We also compare the procedures based on $L^{\infty}$ and $L^{2}$ norms and establish which is the most eficient tool when the uncertainties in the spectral function are taken into account. In the process, we give closed analytical expressions for the functional distances in a generalized weighted $L^{2}$ norm that interpolates almost exactly between $L^{2}$ and $L^{\infty}$.

The plan of the paper is as follows. We start, in Sec. 2, with a brief review of the approach proposed in Ref. **[20]**, defining the minimal functional distances in $L^{\infty}$ and $L^{2}$ norms and presenting the algorithms for their calculation. Sec. 3 contains two new mathematical developments important for applications: in subsection  3.2 we prove that the minimal distance based on the general weighted $L^{2}$ norm can be written down in a closed analytic form, and in subsection 3.2 we derive a suitable approximation of the minimal distance based on $L^{\infty}$ norm by a class of weighted $L^{2}$ norms. In Sec. 4 we briefly review the toy model and describe the data generation with ALEPH-based covariances. Sec. 5 contains our main results and Sec. 6 is devoted to our conclusions.

2 Theoretical framework

We begin with a short presentation of the work performed in Ref. **[20]**. The main idea is to quantify the difference, along a contour in the complex $s$ plane, between the QCD description $\Pi_{\rm QCD}(s)$ of a correlator of light quark currents and its true value $\Pi(s)$. By QCD description one understands the perturbative part, the contribution from OPE condensates and possible duality violations:

[TABLE]

where $\Pi_{\rm OPE}$ encompasses both the purely perturbative contribution (or dimension $D=0$ contribution) and the power corrections.

For definiteness the contour was taken as the circle $|s|=s_{0}$ shown in Fig. 1, where $s_{0}$ is sufficiently large such that the QCD approximant is valid. Measuring the distance by the $L^{\infty}$ norm **[21]**, we consider the quantity

[TABLE]

where $\Pi(s)$ is the “true”, physical function, known to be an analytic function in the $s$-plane cut along the real axis for $s\geq 4m_{\pi}^{2}$, which satisfies the Schwarz reflection condition $\Pi(s^{*})=\Pi^{*}(s)$. In addition, its discontinuity

[TABLE]

is known experimentally on a limited energy range, $4m_{\pi}^{2}\leq s\leq s_{0}$.

The exact value $\delta_{\rm exact}$ of $\delta$ cannot be computed in QCD for lack of the true $\Pi(s)$: the properties stated above do not specify uniquely the function, but define a whole class of admissible functions to which the physical one must belong. If we define

[TABLE]

where the minimization is with respect to all functions $\Pi$ in this admissible class, then it follows that

[TABLE]

As shown in Ref. **[20]**, the quantity $\delta_{0}$ can be calculated by applying a duality theorem in functional optimization **[21]**. For completeness, we present below the main steps of the proof. We make first the simple change of variable

[TABLE]

which maps the domain shown in Fig. 1 onto the unit disk $|z|\leq 1$ cut along a segment of the real axis. Various classes of analytic functions have been defined in the canonical domain $|z|\leq 1$. In particular, adopting the functional distance (4) we are lead naturally to the class $H^{\infty}$ of functions $F(z)$ analytic inside the disk and bounded on the circle $|z|=1$, with the $L^{\infty}$ norm defined as the supremum of the modulus along the boundary $|z|=1$:

[TABLE]

We consider also the class $H^{2}$ of functions analytic inside the disk and of finite $L^{2}$ norm on the frontier $|z|=1$, where

[TABLE]

and the more general class $H^{2}_{w}$ of analytic functions of finite $L^{2}_{w}$ norm, where $L^{2}_{w}$ is the more general norm defined as

[TABLE]

in terms of a weight $w(\zeta)$ given on the boundary $\zeta=e^{i\theta}$ of the unit circle.

As shown in Ref. **[20]**, the problem (4) can be written in the equivalent form

[TABLE]

Here the minimization is performed with respect to all the functions $g(z)$ analytic in the disk $|z|<1$ and bounded on the frontier, and $h$ is a known complex function defined on the boundary $\zeta=\exp(i\theta)$ of the unit circle by

[TABLE]

where $\eta>0$ is an arbitrary parameter introduced for technical reasons, which does not appear in the final result (see **[20]** for more explanations).

The solution of the problem (10) has been obtained in Ref. **[20]** by means of a duality theorem in functional optimization **[21]**. This theorem reads

[TABLE]

where the functions $w(z)$ and $f(z)$ belong to the unit sphere $S^{2}$ of $H^{2}$, i.e. are analytic in $|z|<1$ and satisfy the conditions

[TABLE]

where the norm is defined in Eq. (8).

We recall that all the functions considered here are real analytic, i.e. satisfy the reflection property written as $f(z^{*})=f^{*}(z)$. Therefore, if one writes the Taylor expansions

[TABLE]

coefficients will be real and, due to (13), will satisfy the conditions

[TABLE]

The supremum in the right-hand side of Eq. (12) can be calculated by means of a rather simple numerical algorithm, as shown in Ref. **[20]**. Namely, let ${\cal H}$ be the Hankel matrix defined as

[TABLE]

in terms of the real coefficients

[TABLE]

The coefficients $c_{n}$ defined in Eq. (17) are actually the negative-frequency Fourier coefficients which measure the “non-analytic” part in $|z|<1$ of the complex function $h(\zeta)$ defined in Eq. (11). One may recognize in them the moments used in traditional finite-energy sum rules based on a Cauchy integral relation for $\Pi$ multiplied with a power of $s$ along the contour of Fig. 1. Then, $\delta_{0}$ is obtained as the spectral norm

[TABLE]

i.e.  the square root of the greatest eigenvalue of the positive-semidefinite matrix ${\cal H}^{\dagger}{\cal H}$.

In the numerical calculations, the matrix ${\cal H}^{\dagger}{\cal H}$ is truncated at a finite order $m=n=N$, using the fact that for large $N$ the successive approximants tend toward the exact result (for a formal proof of convergence see Appendix E of Ref. **[22]** and for numerical tests see Ref. **[23]**). By the duality theorem, the initial functional minimization problem (4) is thus reduced to a rather simple numerical computation.

One can define also the minimal functional distance $\delta_{2}$ based on the $L^{2}$ norm

[TABLE]

which can be written in the $z$ variable as

[TABLE]

for the same function $h$ defined in Eq. (11). The solution of this minimization problem has the simple form **[20]**

[TABLE]

in terms of the same coefficients $c_{n}$ defined in Eq. (17).

More generally, we consider the functional distance based on the more general norm $L^{2}_{w}$ defined in Eq. (9), when instead of (20) we must solve the problem

[TABLE]

where $w$ is a suitable weight. It can be shown, without loss of generality, that $w$ can be taken as the boundary value of an outer function **[21]**, i.e. a function $w(z)$ analytic and without zeros in $|z|<1$. It is easy to show then that the solution of the problem (22) is

[TABLE]

where the real numbers $d_{n}$ are the weighted moments

[TABLE]

depending on the function $w(z)$. The quantity $\delta_{2}$ defined in the standard $L^{2}$ norm is obtained from these relations for $w(z)\equiv 1$.

In practice, as in the calculation of $\delta_{0}$ by means of Eq. (18), the infinite sums in Eqs. (21) and (23) are truncated after a finite number of terms and the convergence towards the values $\delta_{2}$ and $\delta_{2,w}$ is tested numerically. Actually, as we will show in the next section, the infinite summation in the general case (23) can be performed exactly and the minimal distance $\delta_{2,w}$ can be written in a closed analytic form.

3 New mathematical developments

3.1 Closed analytic form of $\delta_{2,w}$

A compact analytic form for the quantity $\delta_{2,w}$ can be obtained easily by performing the summation upon $n$ in the expression (23). For convenience we write the real coefficients $d_{n}$ defined in Eq. (24) as

[TABLE]

where the significance of the terms is obvious. Then we obtain from (23):

[TABLE]

Using the expression of $d_{n,\sigma}$ from (24), the first sum in Eq. (26) can be written immediately as

[TABLE]

The second sum in Eq. (26) is written in a convenient form by using for the coefficients $d_{n,{\rm QCD}}$ the expression

[TABLE]

which is explicitly real. Using further the reality property of the functions $w$ and $\Pi_{\rm QCD}$, i.e. the relation $w(z^{*})\Pi_{\rm QCD}(s^{*})=w^{*}(z)\Pi^{*}_{\rm QCD}(s)$, the integration interval can be reduced to $(0,\,\pi)$. Thus we obtain after a straightforward calculation

[TABLE]

We note that the end singularities at $x=y=1$ in the integrand of (27) and at $x=1$, $\theta=0$ in the integrand of (29) are logarithmically integrable.

The last sum in Eq. (26) can be written as

[TABLE]

where $P$ denotes the principal part and

[TABLE]

is the phase of the complex function $w(z)\Pi_{\rm QCD}(s)$ on the circle $|s|=s_{0}$, i.e. on $|z|=1$.

For the numerical evaluation it is more convenient to write the second term in the r.h.s. of Eq. (3.1) in the equivalent form:

[TABLE]

By collecting the terms in Eqs. (27), (29), (3.1) and (3.1) we obtain the final expression of the squared of $\delta_{2,w}$

[TABLE]

All the integration intervals can be further reduced to $(0,\,\pi)$ by taking into account, as explained above, the reality property of the functions, which implies that the real (imaginary) parts are even (odd) functions of $\theta$.

3.2 Approximation of $L^{\infty}$ norm by a suitable class of $L^{2}_{w}$ norms

We show now that it is possible to approximate the minimal distance $\delta_{0}$ measured by the $L^{\infty}$ norm by a class of minimal distances $\delta_{2,w}$ defined by the weighted $L^{2}_{w}$ norms. We follow an argument put forward for the first time in Refs. **[24, 25]**, which is based on the duality theorem Eq. (12) applied for solving the original minimization problem (10).

We note that the r.h.s. of Eq. (12) requires the calculation of the supremum upon two sets of functions, $w(z)$ and $f(z)$, which are analytic in the unit disk $|z|<1$ and of $L^{2}$ norm bounded by 1. The idea is to calculate first the supremum upon one class of functions, say $f$, keeping the other one fixed. We note that the r.h.s. of Eq. (12) can be written as

[TABLE]

where $f_{n}$ are the Taylor coefficients defined in Eq. (14) and $d_{n}$ are negative-frequency Fourier coefficients of the product $w(\zeta)h(\zeta)$, given by the weighted moments (24). Then Eq. (12) becomes

[TABLE]

The supremum upon the coefficients $f_{n}$ subject to the second condition (15) can be evaluated immediately by Cauchy-Schwarz inequality, leading to

[TABLE]

where the dependence of the coefficients $d_{n}$ on the weight $w$ is given in Eq. (24). Finally, by using (23), we write this relation as

[TABLE]

We emphasize that this is an exact relation, which states that the minimal distance $\delta_{0}$ in the $L^{\infty}$ norm is the largest value from the class of distances $\delta_{2,w}$ in the weighted $L^{2}_{w}$ norm, for all the weight $w$ subject to the first condition (13).

Of course, the problem is not yet solved, we still have to calculate the supremum in (37). The procedure makes sense if one can find a suitable, simple parametrization of the functions $w$, such that the maximization upon this limited class approximates well the exact $\delta_{0}$. It turns out that such a choice exists **[24, 25]**: the main observation is that one can obtain approximately the maximum modulus of a function on a certain interval by computing the normalized integral of its modulus squared in a variable that dilates the region where the modulus of the function reaches its maximum. Therefore, one can approximate the $L^{\infty}$ norm (7) of an arbitrary function by an $L^{2}$ norm (8) defined on the integration range distorted by a suitable change of variable. In order to obtain it, we consider the conformal mapping of the unit disc $|z|\leq 1$ onto itself, achieved by the so-called Blaschke transformation **[21]**

[TABLE]

where $a$ is an arbitrary parameter with $|a|<1$. Since we consider real analytic functions, one can restrict $a$ to real values. The transformation (38) maps in particular the unit circle $|z|=1$ onto itself. This change of variable in the $L^{2}$ norm (8) introduces the Jacobian $|dz^{\prime}/dz|$, which corresponds to a weight function $w(z)$ in the weighted $L^{2}_{w}$ norm (9), of the form

[TABLE]

where $a$ is a real parameter in the range $(-1,1)$. It is easy to check that this function satisfies the first condition (13).

By the above remark, the functional supremum in Eq. (37) was reduced to a maximization with respect to a single real parameter $a$. The minimal distance $\delta_{0}$ can thus be calculated approximately by a relatively simple algorithm: first one calculates the minimal distance $\delta_{2,w}$ given in (23), with the particular choice (39) of the weight. Then the parameter $a$ is varied in the range $(-1,1)$ and the largest value of $\delta_{2,w}$ is retained. This problem is numerically quite simple, especially since, as shown in the previous subsection, $\delta_{2,w}$ for an arbitrary weight $w$ can be written in an analytic compact form.

Some hints on the optimal value of the parameter $a$ are obtained from the specific properties of the input. Thus, we note that for values of $a$ close to 1, the function $w(z)$ is large near $z=1$, i.e. near $s=s_{0}$, both on the circle and on the real axis. Therefore, in this case the weighted norm (22) is dominated by the region of the circle near the timelike axis. One can expect that such values of $a$ would be useful in order to detect DVs that are large only near the timelike axis. We shall test these expectations in the numerical studies reported in Sec. 5.

4 Toy data generation

The main goal of the formal developments presented in the previous sections is to provide tools for the validation of DV models using information from the spectral function data. It is therefore instrumental to test the procedure with toy data sets generated from a realistic model for which the DVs are known exactly. Part of the work described here is an extension of analytical results of Ref. **[20]** to a more realistic situation, where the spectral function comes in the form of a binned data set subjected to statistical fluctuations. With the application to ALEPH data in mind, we shall consider data sets that are obtained from a realistic covariance matrix. In this section we discuss the central model used for the exercises performed in this work as well as how we construct our covariance matrix. With these two ingredients, we have full control over the problem and can perform simulations in order to understand how the procedure can be applied to real data.

We start with a brief review of the model that we employ for this exercise.111An extended discussion of the model in the present context can be found in Ref. [20]. The model was introduced in Ref. **[6]**, based on previous ideas from Refs. **[3, 2, 4]**. To be concrete, here we focus on the vector spectral function. The description is based on a “Regge tower” of resonances and upon including the $\rho$ meson pole into the tower, the correlator $\Pi(s)$ assumes the following exact form

[TABLE]

where we defined

[TABLE]

and $\psi(v)=\Gamma^{\prime}(v)/\Gamma(v)$ is the Euler digamma function. We employ the following set of parameters:

[TABLE]

which provides, for our purposes, a realistic description of the experimentally observed vector spectral function of the QCD correlator.

The asymptotic expansion of the digamma function can be used in order to obtain an OPE-type description of the correlator $\Pi_{\rm model}(s)$. Truncating the expansion at an order $N_{\rm OPE}$ it reads

[TABLE]

The first term corresponds to the “purely perturbative” part and the other terms are power corrections, akin to the condensate contributions of QCD. The explicit expression of the coefficients that appear in $\Pi_{\rm OPE}$ are

[TABLE]

with $B_{k}(x)$ representing Bernoulli polynomials.

The asymptotic expansion of Eq. (43) is not accurate near the timelike axis, as in the case for the OPE in QCD. For large enough $|s|$ and ${\rm Re}(s)>0$ the description can be improved taking into account the DVs. In practice, the DV term can be obtained from the reflection property of the digamma function **[6, 20]**. The following modified approximant is thus obtained

[TABLE]

valid for large enough $|s|$ and for ${\rm Re}(s)>0$. The DV contribution is given in the first quadrant (${\rm Im}(s)>0$* and ${\rm Re}(s)>0$) by*

[TABLE]

and can be defined in the lower half-plane using Schwarz reflection as $\Pi_{\rm DV}(s^{*})=\Pi^{*}_{\rm DV}(s)$. For ${\rm Re}(s)\leq 0$ this correction is assumed to vanish.

Comparing the modulus of the exact function, Eq. (40), along the upper semi-circle $s=s_{0}e^{i\theta}$, $\theta\in(0,\pi)$, with its approximants, Eqs. (43) and (45), one learns that the truncated OPE-type expansion of Eq. (43) provides an accurate description except close to the timelike axis ($\theta=0$), as expected in QCD. The addition of the DV term fixes this deficiency and the approximated description of Eq. (45) becomes excellent also in the vicinity of the timelike axis. (We refer to Ref. **[20]** for a visual account of this comparison.)

For the numerical exercises described in this work we use the model of Eq. (40) as our central description. Hence, the OPE for the model and the DV contribution are exactly known and are given by Eqs. (43) and (46), respectively. The values of the vector spectral function for toy data generation are obtained, thus, from

[TABLE]

In order to mimic the experimental situation, the interval $[0:m_{\tau}^{2}]$ is split in $N_{\rm b}$ bins and the central value of each bin is obtained from a statistical distribution that fluctuates around the values of Eq. (47) calculated at the center of each bin. We turn now to the issue of the covariance matrix that governs these fluctuations.

Our toy data generation is performed having in mind the application to the ALEPH spectral functions **[16]**. It is therefore desirable that the covariances used reflect those of ALEPH data sets. One could simply adopt the ALEPH covariances as such, since they are publicly available **[27]**, and generate toy data sets following a statistical distribution given by this matrix, together with the central values of Eq. (47). The price to pay is that one would have to use the ALEPH binning of the interval $[0:m_{\tau}^{2}]$. In the most recent version of the data sets, due to an improved unfolding procedure, an adaptive binning was used which results in bins with different widths, notably with larger bins towards the edge of the spectrum. Here we prefer to adopt a fixed bin width, for simplicity, and we choose $N_{\rm b}$ such as to have more bins than ALEPH towards the end-point of the spectrum. This allows us to have a finer description at higher energies. The accompanying realistic covariances are obtained from a numerical interpolation of the ALEPH covariance matrix for the vector channel. In this way, we preserve a fixed binning together with a covariance matrix that has all the main properties of ALEPH’s, namely, strong correlations between neighbouring bins, larger uncertainties towards the end-point of the spectrum and, of course, uncertainties that are of the same order of those of ALEPH’s data.

In the present work we adopt $N_{\rm b}=80$ (which is in line with what is used in the experimental analyses **[28, 29, 16]**), the central values of Eq. (47), and the covariance matrix obtained from a numerical interpolation of the ALEPH covariances obtained from **[27]**, as described above. Toy data sets can then be generated from a multivariate Gaussian distribution. An example of a data set generated in this way is displayed in Fig. 2. We show also the central values of Eq. (47) for comparison. In this figure, the strong correlations are clearly visible, mainly towards the end-point of the spectrum, where the uncertainties are also larger.

When using data sets for the calculation of the functional distances $\delta$ discussed in Secs. 2 and 3, weighted integrals over the spectral function such as those entering Eq. (17) must be discretized. We are going to adopt integration by rectangles, as is usual when dealing with this type of integrals of the spectral functions **[16, 15]**. However, weight functions $w_{n}$ with high powers of the energy variable appearing, for example, in Eq. (17), vary strongly within a bin. It is therefore necessary to average over the weight function inside a bin to improve the numerical result.222In the case of ALEPH data this prescription is sometimes used due to the large bin widths of the right-most bins [15]. The numerical counterpart of a typical integral reads then

[TABLE]

*where $\bar{s}_{i}$ is the value of $s$ at the center of the *$i$th bin, $\Delta_{b}$ is the fixed bin width, and $[s_{0}]$ represents the index of the last included bin — here we always work with $s_{0}$ values that correspond to the right edge of a bin. The same was applied for the calculation of the relevant integrals which appear in the analytic form of $\delta_{2,w}$ derived in Sec. 3.1. We have tested that this algorithm provides enough accuracy for the explorations performed in this work.

5 Results

We apply now the functional-analysis based tools to test in practice the description of DVs. To illustrate the potential of the method, a useful approach is to introduce in the approximate description (45) of the correlator a strength parameter $\mu$ that allows one to tune the contribution of the DVs. Formally, we do this by using for $\Pi_{\rm QCD}$ in the formalism presented in Secs. 2 and 3.1, instead of Eq. (1), the more general expression

[TABLE]

where the true value of the strength parameter is $\mu=1$. As in Ref. **[20]**, to simulate the situation of the light-quark correlators in QCD, we take $\Pi_{\rm OPE}$ as the asymptotic expansion (43) of the exact model truncated after $N_{\rm OPE}=5$ terms. For $\Pi_{\rm DV}$ we take the prediction (46) of the model.

In Ref. **[20]**, it was shown by means of analytical computation that $\delta_{0}$ has a sharp minimum at the correct value $\mu=1$, when one employs the description of DVs that follows from the model used for $\Pi(s)$. The alternative quantity $\delta_{2}$ displayed also a minimum at the correct value of $\mu$, but this minimum was found to be shallower **[20]**. In this section we investigate the impact of the use of spectral function data with realistic covariances to the above findings. It will be interesting to make use of the weighted $L^{2}_{w}$ norm, since it permits a continuous and almost exact interpolation between the $L^{2}$ and the $L^{\infty}$ norms, as well as a study of other weighted norms, such as “pinched” norms. The analytical results obtained for $\delta_{2,w}$ shall also be instrumental in this analysis.

5.1 Comparison between $L^{\infty}$ and $L^{2}_{w}$ norms

We check first on the toy model the approximation of the minimal distance $\delta_{0}$ based on $L^{\infty}$ norm by the distances $\delta_{2,w}$ based on the norms $L^{2}_{w}$, using the particular class of weights given in (39). In this discussion we use, as in Ref. **[20]**, the exact spectral function of the model, with no errors, and the OPE expansion truncated at $N_{\rm OPE}=5$.

From Fig. 3, which shows the modulus of $\Pi_{\rm DV}(s)$ as a function of $\theta$ on the first quadrant of the circle $s=s_{0}\exp(i\theta)$, one can see that the DV part of the model is strongly peaked towards the Minkowskian axis. Therefore, following the discussion at the end of Sec. 3.2, a weight strongly peaked towards $s=s_{0}$ (i.e. $\theta=0$) is expected to give the best approximation of $L^{\infty}$ norm by weighted $L^{2}_{w}$ norms for this model. As shown below, the expectation is confirmed.

In Fig. 4 we show the variation with $\mu$ of several functional distances, calculated with the algorithms based on Fourier coefficients truncated at $N=100$. We set in this exercise the radius of the circle in Fig. 1 to $s_{0}=2.76\,\,\text{GeV}^{2}$, but the results are similar for other choices, including $s_{0}=m_{\tau}^{2}$.

As usual, $\delta_{0}$ denotes the minimal distance measured in $L^{\infty}$ norm, calculated from the norm (18) of the Hankel matrix, Eq. (16). For calculating $\delta_{2,w}$, we used the truncated sum (23), with the expression (24) of $d_{n}$ and a weight $w$ of the form (39) with $a=0.96$. This weight drastically dilates the region near $\theta=0$ on the circle, increasing its contribution to the $L^{2}_{w}$ norm. One can see that, for this choice of the weight, the distance $\delta_{2,w}$ practically coincides with $\delta_{0}$. Both curves are steeper than the standard distance $\delta_{2}$, which corresponds to the weight $w=1$, as already remarked for $\delta_{0}$ in Ref. **[20]**. The figure shows also that the minimal distance $\delta_{2,{\rm pd}}$, calculated with a pinched weight of the form333One refers as “pinched” to weight functions that have a zero for $s=s_{0}$.**

[TABLE]

is much less sensitive to the variation of $\mu$, which is not surprising since this type of weight suppresses the region where the DV term is nonzero.

The optimal value of the parameter $a$ was found empirically, by computing $\delta_{2,w}$ for several values of $a$ close to 1, and keeping the value leading to the best approximation of $\delta_{0}$. Of course, the best value achieving the supremum in (36) depends also on the other ingredients of the input. Thus, for a different number $N$ of Fourier coefficients taken into account in the calculation of the norms a slightly different value of $a$ might yield the best approximation. Also, a slightly different optimal value of $a$ is expected if the input spectral function is slightly changed. This remark will be useful for understanding the results of the simulations performed below, which take into account the uncertainties on the spectral function.

5.2 Stability and comparison with the summed results

We start the simulations by investigating the computation of the distance $\delta_{2}$ measured with $L^{2}$ norm from the truncated version of Eq. (21) — which can be viewed as a special case of Eq. (23) with the appropriate choice of the weight $w$. For a given data set, one can compute the value of $\delta_{2}$, or more generically, of $\delta_{2,w}$, after a truncation of the infinite sums at the $N$-th term. The value of the norm can then be minimized numerically with respect to the strength parameter $\mu$. Due to the statistical fluctuations, each toy data set that is generated yields a different value of $\mu$. We repeat this procedure for 5,000 different data sets, in a reproducible way, in order to obtain the statistical distribution of the parameter $\mu$. The final value of $\mu$ can be read off from the distributions. We quote central values given by the medians and uncertainties defined by $68$% confidence levels, but the distributions are, to a very good approximation, Gaussian, as illustrated with the histogram shown in Fig. 5.

In Table 1, we show the dependence of the best values of $\mu$ on the number $N$ of included terms in the truncated version of Eq. (21). In this table we choose $s_{0}=2.76$ GeV2*, which avoids some of the bins with larger uncertainties (see Fig. 2). One can conclude from this table that the convergence of the results seems to be satisfactory; with a few hundred terms in the sum the results are already stable. Furthermore, the exercise indicates that with a realistic data set the error in $\mu$ is such that we are able to detect the presence of DVs, i.e. $\mu=1$, from their absence, $\mu=0$, in a statistically meaningful way.*

We now turn to the dependence of the predictions on the choice of $s_{0}$. The fact that the last few bins suffer from a much larger uncertainty, combined with the decrease of the DV contribution at higher $s$, has the consequence that the choice of larger $s_{0}$ produces less precise determinations of $\mu$. In Table 2 we compare values of $\mu$ obtained from $\delta_{0}$ and $\delta_{2}$ from the truncated versions of Eqs. (18) and (21), respectively, for different values of $s_{0}$ (we use $N=300$ terms in the sums). Two main conclusions can be drawn from this table. First, as expected, the uncertainties are larger when $s_{0}$ is chosen to be closer to the edge of the spectrum. For $\delta_{2}$, the determination loses statistical significancy rather fast when the last few bins are included, and at $s_{0}=m_{\tau}^{2}$ one can no longer distinguish in a meaningful way the central value $\mu=1.0$ from the absence of DVs. Second, the use of $\delta_{0}$ leads to broader $\mu$-distributions and hence to larger uncertainties. This effect is small for lower $s_{0}$ values and the results are essentially undistinguishable from those obtained using $\delta_{2}$. At $s_{0}=m_{\tau}^{2}$, however, the uncertainty is the double of the $\delta_{2}$ counterpart.444In order to obtain the results of Table 2 for $s_{0}=m_{\tau}^{2}$ it becomes important to allow for negative central values in the toy data spectral functions. These are rare, but do occur for the last few bins where the uncertainties, following the recent ALEPH reanalysis, are rather large. * We conclude that the deeper minimum of $\delta_{0}$ with respect to $\mu$ observed in the analytical calculations of Ref. **[20]** (and seen also above in Fig. 4) does not translate into a narrower $\mu$-distribution when the errors on the spectral function are taken into account.*

A final validation of the results obtained in Tabs. 1 and 2 can be obtained using the closed analytical form of $\delta_{2,w}$ derived in Sec. 3.1. The use of the weighted norm $L^{2}_{w}$ is particularly convenient as it allows for an almost exact interpolation between the $L^{2}$ and $L^{\infty}$ norms and, at the same time, is amenable to a fully analytical treatment of the minimization problem.

Using the decomposition (49) of $\Pi_{\rm QCD}$, we can write Eq. (3.1) as a quadratic polynomial of $\mu$ of the form

[TABLE]

where

[TABLE]

and the calculable coefficients $b_{i}$ can be read off from (3.1). In particular, it is easy to see that only $b_{0}$ depends on the spectral function, the coefficients $b_{1}$ and $b_{2}$ involving only the values of the theoretical expressions $\Pi_{\rm OPE}(s)$ and $\Pi_{\rm DV}(s)$ on the circle $|s|=s_{0}$.

The optimal value of $\mu$, which achieves the minimum of (51), is obtained in a straightforward manner as

[TABLE]

This formulation turns out to be very convenient for numerical simulations, which can be done directly on the spectral function, avoiding the calculation of many experimental moments and the issue of truncating infinite sums. In practice, the coefficients $b_{1}$ and $b_{2}$ and the function $F_{\rm DV}$ are calculated only once, being fixed during data generation. This amounts to a considerable reduction of the computational time required by the statistical simulations.

In Tab. 3 we display the ranges of $\mu$ obtained from simulations using the analytic result (53), for $s_{0}=2.76\,\text{GeV}^{2}$ and $s_{0}=m_{\tau}^{2}$ and three choices for the weight: $w=1$, which corresponds to the standard norm $L^{2}$ with the minimal distance $\delta_{2}$, the expression (39) with $a=0.96$, expected to approximate well the $L^{\infty}$ norm, and the pinched weight (50). We remark the perfect agreement between the results quoted in Table 3 for the weight $w=1$ and the values obtained from $\delta_{2}$ in Table 2 for the same values of $s_{0}$. This validates the convergence of the results based on truncated sums of Fourier coefficients. We also remark a very good agreement between the results from $\delta_{0}$ and those from its approximated version in the second column of Tab. 3. (We discuss these results further in the next section.)

As seen from Table  3, in all cases the central value of the parameter $\mu$ coincides with the true central value. This result was expected, having in view the precise theoretical input used along the circle in our study. The confirmation of this expectation is a test of the numerical algorithms used in the calculations. In particular, the integral in (53) had to be computed using the improved algorithm described in (48), i.e. the product $w(x)F_{\rm DV}(x)$ was integrated exactly over each bin. On the other hand, the uncertainties quoted in the various entries of Table  3 are quite different. The explanation of these results and their relevance for the application of the method to real data will be discussed in the next subsection.

5.3 Discussion

The two weights, used in the simulations with the analytic form of $\delta_{2,w}$ beside $w=1$, are quite extreme: the expression (39) with $a=0.96$ strongly enhances the contribution of the region near the point $s=s_{0}$ on the circle shown in Fig. 1. Since, as seen from Fig. 3, the magnitude of $\Pi_{\rm DV}$ is strongly peaked near the timelike axis, the corresponding minimal distance $\delta_{\rm 2,w}$ will be very sensitive to the variation of the strength parameter $\mu$, as seen from Fig. 4. The same figure shows also that for this weight the norm $L^{2}_{w}$ approximates well the $L^{\infty}$ norm. On the contrary, the weight (50) suppresses the region near $s=s_{0}$, which explains the low sensitivity of the corresponding distance $\delta_{\rm 2,pd}$ to the variation of $\mu$, visible in Fig. 4.

The above remarks refer to a fixed spectral function $\sigma(s)$. When this quantity is varied within errors during the simulations, the two extreme weights respond in a different way. The weight (39), which enhances also the region near $s=s_{0}$ on the real axis, will be more sensitive to the variation of the input data, since the errors are larger towards the end of the spectrum. For the lower value $s_{0}=2.76\,\text{GeV}^{2}$, when the errors are still moderate, the effect of the variation of the input data turns out to be comparable to the opposite effect produced by the larger sensitivity to the variation of $\mu$. As a consequence, the spread of the $\mu$-distribution for $w$ of the form (39) with $a=0.96$ is comparable to that obtained with $w=1$, as seen from the second and third columns of the first row of Table 3. On the other hand, the pinched weight $w_{\rm pd}$ ensures a low sensitivity of $\delta_{\rm 2,pd}$ to the variations of the spectral function produced by the errors. However, the low sensitivity of the same quantity to the variation of the strength parameter $\mu$ leads to an overall large spread of the statistical distribution, which explains the larger error quoted in the last column of the first row of Table 3.

For $s_{0}=m_{\tau}^{2}$, the large errors of the input data in the last bins lead to the large uncertainties on $\mu$ quoted in the second line of Tab. 3. In this case, the detection of DVs in a significant way from the pseudodata is not possible. For the weight (39) with $a=0.96$, the great sensitivity with respect to the input data near the upper end of the spectrum exceeds the opposite large sensitivity to the variation of $\mu$. The resulting $\mu$ has a larger uncertainty than that obtained with the standard $L^{2}$ norm. In the case of the pinched weight $w_{\rm pd}$, the suppressing effect on the large errors of the last bins compensates the spread produced by the low sensitivity to $\mu$ variation. The overall effect is that for $s_{0}=m_{\tau}^{2}$ the spreads on $\mu$ obtained with the two extreme weights are comparable.

A last remark concerns the relation between the weighted $L^{2}_{w}$ norms and the $L^{\infty}$ norm. One can see that for $s_{0}=2.76\,\text{GeV}^{2}$ the $\mu$ distribution quoted in Table 3 for the weight of the form (39) with $a=0.96$ coincides with that obtained from the distance $\delta_{0}$ measured in $L^{\infty}$ norm, given in Table 2. However, for $s_{0}=m_{\tau}^{2}$, the standard deviation on $\mu$ quoted in the third column of Table 3 is somewhat smaller than that obtained with $\delta_{0}$ in Table 2.

To understand this small difference, we recall that the simulations reported in Table 3 were performed with a fixed value, $a=0.96$, in the expression (39). But, as discussed in Sec. 3.2, the optimal choice of $a$ achieving the supremum in (36) depends also on the input spectral function. For $s_{0}=2.76\,\text{GeV}^{2}$, when the bins with large errors are excluded, this dependence affects in an almost unobservable way the simulations. However, for $s_{0}=m_{\tau}^{2}$ the variation of the input can be considerable due to the large errors in the last bins. In this case the weight (39) with $a=0.96$ is not always the optimal weight leading to the precise approximation of $\delta_{0}$ according to Eq. (36). Therefore, the result presented in Table 3 only illustrates the use of the norm $L^{2}_{w}$ for a rather extreme weight, inspired from the $L^{\infty}$ norm but not reproducing exactly its results. When the errors are large, the simulations using the $L^{\infty}$ norm must resort to the exact algorithm (18) with the Hankel matrix (16).

6 Summary and conclusions

In the present paper we continued the investigation of the functional-analyses tools proposed in Ref. **[20]** for detecting DVs from measurements of the spectral functions of the QCD correlators. The aim was to evaluate the potential of the method when the spectral function comes in the form of binned data with realistic covariances. We performed the analysis still in the context of the toy model for the correlator $\Pi(s)$ considered in **[20]**, in which we allowed for uncertainties described by the covariances obtained from the publicly available ALEPH spectral functions. In this way we had full control over the problem and the outcome of the method could be checked against the expected results.

The paper contains also some theoretical developments of the approach proposed in Ref. **[20]**. In addition to the functional distances based on $L^{2}$ and $L^{\infty}$ norms, already discussed in Ref. **[20]**, we introduced a general class of weighted norms $L^{2}_{w}$, which are instrumental for several reasons. First, as shown in Sec. 3.1, we were able to obtain a closed analytic expression for the minimal functional distance $\delta_{2,w}$ measured in this norm, thereby avoiding truncated sums of Fourier coefficients. Second, these norms provide an interpolation between two extreme cases: the pinched weights familiar from phenomenological works, and the opposite class of weights which, as discussed in Sec. 3.2, provide a good approximation of the functional distance measured in $L^{\infty}$ norm.

To investigate the potential of the method for the detection of DVs we introduced, in the spirit of Ref. **[20]**, a strength parameter $\mu$ that quantifies the DV contribution to $\Pi(s)$ according to (49), the true value of this parameter being $\mu=1$. As in **[20]**, we define the optimal $\mu$ as the value that achieves the minimum of the lower bounds on the functional distances $\delta_{0}$, $\delta_{2}$ or $\delta_{2,w}$, measured in the norms $L^{2}$, $L^{\infty}$ or $L^{2}_{w}$, respectively, between the true correlator and its approximant along the circle $|s|=s_{0}$ in the complex energy plane.

For want of a theoretical statistical interpretation of the minimal distances defined by functional analysis, we performed an empirical study where fake data on the spectral function have been generated in a number of bins. To mimic the experimental situation, we adopted a multivariate Gaussian distribution with covariances inferred from the ALEPH covariance matrix for the vector channel **[16]**. By simulations with 5,000 different data sets, we obtained the statistical distributions of the optimal parameter $\mu$, which were, to a very good approximation, Gaussian.

The main results of these investigations are displayed in Tabs. 1, 2 and 3, where we quote the central values given by the medians and the uncertainties defined by $68$% confidence levels from the corresponding distributions. One can see that the results based on the truncated computation of the norms converge relatively fast, which make their practical use feasible. This could be confirmed using the analytical determination of $\mu$ given in Eq. (53), that avoids the necessity of truncating the sums. We investigated in this framework three types of weights: $w=1$, which corresponds to the standard $L^{2}$ norm, the expression (39) with $a=0.96$, expected to approximate well the $L^{\infty}$ norm, and the pinched weight (50).

We note that in all cases the true value $\mu=1$ of the strength parameter is obtained with high accuracy. Since the theoretical input we use is quite precise, this result represents a good test of the numerical algorithms adopted. In particular, as discussed in Sec. 5, the refined integration rule (48) for calculating either the moments (17) or the quantity (52) must be used for reaching this level of accuracy. On the other hand, the standard deviations, crucial for the extraction of DVs in a significant way, differ for various tests. For the lower value of $s_{0}$, the tests based on the norms $L^{2}$ and $L^{\infty}$ produce comparable uncertainties, with a successful and statistically significant (by three standard deviations) detection of DVs. The test based on pinched weight (50) is however unable to detect DVs in a significant way even at low $s_{0}$. For $s_{0}=m_{\tau}^{2}$, due to the large uncertainties towards the edge of the spectrum, a statistically significant determination of $\mu$ is not possible. All the tests have very large uncertainties, although one may note that the performance of the $L^{2}$ norm is superior to those of the $L^{\infty}$ norm and the weighted $L^{2}_{w}$ norms with the other two weights considered in Table 3.

One might ask what is the relation of the present approach to the standard $\chi^{2}$-type fits used up to now for the phenomenological determination of DVs. The coefficients $c_{n}$ defined in Eq. (17) are actually the moments used in traditional finite-energy sum rules based on a Cauchy integral relation for the correlator $\Pi(s)$ multiplied with a power of $s$ along the contour of Fig. 1. Replacing the approximant $\Pi_{\rm QCD}$ by the exact $\Pi$ and using the exact spectral function $\sigma$, we would have $c_{n}=0$, by analyticity. In practice, the coefficients $c_{n}$ are not zero due to the imperfections of $\Pi_{\rm QCD}(s)$ and to the statistical fluctuations of experimental values of $\sigma(s)$.

In the standard analyses, starting from this remark a few moments are selected and combined for defining a certain “fit quality”, usually a $\chi^{2}$, with an assumed statistic distribution. This allows, by standard techniques of $\chi^{2}$ minimization, the extraction of the parameters of the DV models and their covariances together with the values of other parameters of the OPE, in particular the strong coupling constant and condensates, for example. The limitation of this approach is that only a small number of low-order moments, with known errors and possible correlations, can be included in the fit, due to the fact that in QCD only a small number of power corrections are available. The inclusion of high-order moments must be avoided, as it would introduce unknown high-order condensates **[19]**.

In the present approach, the fact that only a finite number of power corrections are known in the OPE is not essential, because no assumption about the vanishing of specific moments is made. The method exploits the obvious remark that the exact correlator and its approximate representation provided by the QCD calculations available at present are different functions, with different analytic properties. Therefore, the functional distance between them, measured in a certain norm, must exceed a rigorous lower bound. As seen from the algorithms of calculating this lower bound, all the moments contribute to it, irrespective of the number of power corrections available in the OPE. The individual moments are actually not relevant, since the minimal functional distance measured in the general $L^{2}_{w}$ norm is proved to have the analytic expression (3.1) directly in terms of the spectral function, which obviates the need for the computation of the moments that appear in Eq. (17).

Of course, there is a price to pay, and this is the fact, already mentioned above, that the functional distances lack a definite statistical interpretation. Therefore, in order to extract optimal parameters with definite confidence intervals, one must resort to simulations based on fake-data generation. Related to that point is the fact that the data covariances enter the procedure solely in the Monte Carlo error propagation and do not affect the evaluation of the functional distances (as they would in a $\chi^{2}$ analysis). The extraction of several free parameters in this framework is feasible in principle, but may be complicated in practice. Therefore, in this work we applied the functional approach for the extraction of only one parameter, the strength $\mu$ of the DV term, assuming that all the other parameters are known. The results show that, restricting $s_{0}$ to values slightly below $m_{\tau}^{2}$, in order to avoid the large errors near the end-point of the spectrum in $\tau$ decays, the functional approach is able to detect, in a statistically significant way, the presence of DVs in realistic spectral function pseudodata.

The present analysis paves the road for the next step, of testing DV models with real data. The first task will be to detect DVs in the real data employing the formalism described here. This can be carried out along the lines of our study with pseudodata, with the help of a strength parameter and using values of additional QCD parameters extracted from other processes. The method can also be used to compare different models and, although technically more challenging, also allows for the extraction of parameters entering the models. Having in view the importance of detecting and describing DVs in a reliable way, this framework, although somewhat limited, is of interest as a complementary approach to other phenomenological studies.

Acknowledgements

The authors would like to express special thanks to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and support during the commencement of this work. We would like to thank Maarten Golterman and Santi Peris for the careful reading of the manuscript. DB thanks the kind hospitality of IFAE and the Universitat Autònoma de Barcelona where this work was finalized. DB’s work was supported by the São Paulo Research Foundation (FAPESP) grant 2015/20689-9, by CNPq grant 305431/2015-3, and by the Alexander von Humboldt Foundation. IC’s work was supported by ANCS, Contract PN 16 42 01 01/2016.