text:
# 1 Introduction

## 1 Introduction

As the predictions of QCD become increasingly precise and more high quality data is available, the theoretical uncertainties associated with the analysis of the data are often found to be dominant and thus have come under increasing scrutiny. In this work we continue previous efforts  to optimize the analysis of the information contained in the experimental data, taking into account errors and correlations, while introducing the smallest possible theoretical bias. We consider the determination of the QCD vacuum condensates $`𝒪_6,𝒪_8`$ and higher dimensional condensates, which in principle can be extracted in a theoretically solid and experimentally clean way from the hadronic decays of the tau lepton. In practice, however, the situation is far from satisfactory, as revealed by the lack of stability of the value of the nonperturbative condensates obtained from this process by different procedures.

The main source of difficulties can be traced to the fact that conventional extractions of the QCD vacuum condensates from hadronic tau decays involve convolutions of the difference $`v_1(s)a_1(s)`$ of the isovector vector and axial vector spectral functions, which is a purely nonperturbative quantity, that does not converge to the perturbative result (these spectral functions are degenerate within perturbation theory) at energies $`sM_\tau ^2`$ and moreover has large uncertainties in the high $`s`$ region. To obtain a reliable value for the nonperturbative condensate, some convergence method must be applied, implying that the error on the condensates gets tangled with uncertainties within the theoretical assumptions and subject therefore to a variety of sources of theoretical bias. This is a consequence of the fact that the kinematics of the hadronic tau decays constrains the range of energies in which we can evaluate the spectral functions. The main difficulty is that any method to estimate non-perturbative condensates exploits the shape of the spectral functions near and beyond the boundary of the region where the data is available. Final results are then subject to systematics errors associated to the extrapolation of data as well as the way global theoretical constraints, e. g. Weinberg sum rules, are imposed.

In this paper, we approach the determination of nonperturbative condensates, in a way which tries to bypass these difficulties, by combining two techniques. First, a novel bona fide method to take into account experimental errors was proposed and implemented in the context of analysis of Deep Inelastic Scattering data to produce a probability measure in the space of deep inelastic structure functions by means of neural networks . Here, we adapt this method to the parametrization of spectral functions. The second technique we use refines the training of neural networks so as to implement the constraints that represent the QCD chiral sum rules in our neural network parametrization of the spectral function $`v_1(s)a_1(s)`$.

The representation of the probability density given in Ref.  takes the form of a set of neural networks, trained on an ensemble of Monte Carlo replicas of the experimental data, which reproduce their probability distribution. The parametrization is unbiased in the sense that neural networks do not rely on the choice of an specific functional form, and it interpolates between data points, imposing smoothness constraints in a controllable way. Information on experimental errors and correlations is incorporated in the Monte Carlo sample. Errors on physical quantities and correlations between them can then be determined without the need of linearized error propagation. Our final parametrization combines all the available experimental information, as well as constraints from different convolutions of the data, i.e. chiral sum rules, must verify. In this way statistical errors can be estimated and the loss of accuracy due to the extrapolation outside the kinematical region where the data is available is also analyzed.

Hence, we can obtain a determination of the QCD vacuum condensates which is unbiased with respect to the parametrization of the spectral function and the error and correlations propagation. We also try to keep under control all sources of uncertainty related to the method of analysis, and estimate their contribution to the total error. This gives us a determination of the nonperturbative condensates, and simultaneously illustrates the power of a method of analysis based on direct knowledge of a probability density in a space of functions.

This paper is organized as follows: in section 2 we review theoretical tools used in the analysis of non-perturbative effects in hadronic tau decays; in section 3 we present the experimental data that is used in our analysis and in section 4 we introduce the neural network parametrization of spectral functions. Section 5 contains the details and results of our extractions of the vacuum condensates: we explain our choice of training parameters, our error estimations and the consistency tests that we performed; finally section 6 summarizes our conclusions.

## 2 Spectral functions and hadronic tau decays

The tau particle is the only lepton massive enough to decay into hadrons. Already before its discovery, it was predicted to be important for the study of hadronic physics , a study that has been performed extensively at the LEP accelerator. Its semileptonic decays are therefore an ideal tool for studying the hadronic weak currents under clean conditions, both theoretically and experimentally, thanks to the high quality data from LEP. In this first section we will briefly introduce the theoretical foundations that form the basis of the QCD analysis of hadronic tau decays. There exists a huge literature on tau hadronic physics to which the interested reader is directed (see for instance Ref.  and references therein).

Spectral functions are the observables that give access to the inner structure of hadronic tau decays. As parity is maximally violated in $`\tau `$ decays, the spectral functions will have both vector and axial vector contributions. As spontaneous chiral symmetry breaking is a nonperturbative phenomena, this spectral functions are degenerate in perturbative QCD with massless light quarks, so any difference between vector and axial-vector spectral functions is necessarily generated by non-perturbative dynamics, that is, long distances resonance phenomena, being the most relevant the $`\rho (770)`$ and the $`a_1(1260)`$ in the vector and axial vector channels respectively. Therefore, the difference of these spectral functions is generated entirely from nonperturbative QCD dynamics, and provides a laboratory for the study of these perturbative contributions, which have resulted to be small and therefore difficult to measure in other processes where the perturbative contribution dominates. An accurate extraction of the QCD vacuum condensates is not only important by itself but also has many important phenomenological applications, for example in the evaluation of matrix elements in weak decays .

The ALEPH collaboration at LEP measured  these spectral functions from hadronic tau decays with great accuracy, providing an excellent source of precision analysis of nonperturbative effects. As it is well known , the basis of the comparison of theory with data is the fact that unitarity and analyticity connect the spectral functions of hadronic tau decays to the imaginary part of the hadronic vacuum polarization,

$$\mathrm{\Pi }_{ij,U}^{\mu \nu }(q)d^4xe^{iqx}0|T\left(U_{ij}^\mu (x)U_{ij}^\nu (0)^{}\right)|0,$$
(1)
of vector $`U_{ij}^\mu V_{ij}^\mu =\overline{q}_j\gamma ^\mu q_i`$ or axial vector $`U_{ij}^\mu A_{ij}^\mu =\overline{q}_j\gamma ^\mu \gamma _5q_i`$ color singlet quark currents in corresponding quantum states. After Lorentz decomposition is used to separate the correlation function into its $`J=1`$ and $`J=0`$ components,

$$\mathrm{\Pi }_{ij,U}^{\mu \nu }(q)=\left(g^{\mu \nu }q^2+q^\mu q^\nu \right)\mathrm{\Pi }_{ij,U}^{(1)}(q^2)+q^\mu q^\nu \mathrm{\Pi }_{ij,U}^{(0)}(q^2),$$
(2)
for non-strange quark currents one identifies

$$\mathrm{Im}\mathrm{\Pi }_{\overline{u}d,V/A}^{(1)}(s)=\frac{1}{2\pi }v_1/a_1(s).$$
(3)
This relation allows us to implement all the technology of QCD vacuum correlation functions to hadronic tau decays, and provides the basis of the comparison of theory with data.

The basic tool to study in a systematic way the power corrections introduced by nonperturbative dynamics is the operator product expansion. Since the approach of Ref. , the operator product expansion (OPE) has been used to perform calculations with QCD on the ambivalent energy regions where nonperturbative effects come into play but still perturbative QCD is relevant. In general, the OPE of a two point correlation function $`\mathrm{\Pi }^{(J)}(s)`$ takes the form

$$\mathrm{\Pi }^{(J)}(s)=\underset{D=0,2,4,\mathrm{}}{}\frac{1}{(s)^{D/2}}\underset{\mathrm{dim}𝒪=D}{}C^{(J)}(s,\mu )𝒪(\mu ),$$
(4)
where the arbitrary parameter $`\mu `$ separates the long distance nonperturbative effects absorbed into the vacuum expectation elements $`𝒪(\mu )`$, from the short distance effects which are included in the Wilson coefficient $`C^{(J)}(s,\mu )`$. The operator of dimension $`D=0`$ is the unit operator (perturbative series), and we are interested in the dimension $`D6`$ operators. What is relevant for us is that $`D=6`$ is the first non-vanishing non-perturbative contribution, in the limit of massless light quarks, to the $`v_1(s)a_1(s)`$ spectral function and, moreover, it has been shown to be the dominant one. This dominant contribution carries non-trivial four-quark dynamical effects of the form $`\overline{q}_i\mathrm{\Gamma }_1q_j\overline{q}_k\mathrm{\Gamma }_2q_l`$. Additional contributions from a mixed quark gluon condensate as well as a triple gluon condensate are assumed to be small. Therefore, this spectral function should provide a source for a clean extraction of the value of the nonperturbative contributions.

Finally, we review the important paper that QCD chiral sum rules play in the analysis of this process. Sum rules have always been an important tool for studies of non-perturbative aspects of QCD, and have been applied to a wide variety of processes, from Deep Inelastic Scattering to Heavy Quark systems ,. Now we will review one of the classical examples of low energy QCD sum rules, the chiral sum rules. The application of chiral symmetry together with the optical theorem leads to low energy sum rules involving the difference of vector and axial vector spectral functions,

$$\rho _{VA}(s)v_1(s)a_1(s).$$
(5)
These sum rules are dispersion relations between real and absorptive parts of a two point correlation function that transforms symmetrically under $`SU(2)_LSU(2)_R`$ in the case of non strange currents. Corresponding integrals are the Das-Mathur-Okubo sum rule

$$\frac{1}{4\pi }_0^{s_0\mathrm{}}𝑑s\frac{1}{s}\rho _{VA}(s)=\frac{f_\pi ^2<r_\pi ^2>}{3}F_A,$$
(6)
as well as the first and second Weinberg sum rules (WSR)

$$\frac{1}{4\pi ^2}_0^{s_0\mathrm{}}𝑑s\rho _{VA}(s)=f_\pi ^2,$$
(7)
$$_0^{s_0\mathrm{}}𝑑ss\rho _{VA}(s)=0,$$
(8)
where in eq. (7) the RHS term comes from the integration of the spin zero axial contribution, which for massless non-strange quark currents consists exclusively of the pion pole. Finally, there is the chiral sum rule associated with the electromagnetic splitting of the pion masses ,

$$\frac{1}{4\pi ^2}_0^{s_0\mathrm{}}𝑑ss\mathrm{ln}\frac{s}{\lambda ^2}\rho _{VA}(s)=\frac{4\pi f_\pi ^2}{3\alpha }(m_{\pi ^\pm }^2m_{\pi ^0}^2),$$
(9)
where $`f_\pi =(92.4\pm 0.3)`$ MeV  obtained from the decays $`\pi ^{}\mu ^{}\overline{\nu }_\mu `$ and $`\pi ^{}\mu ^{}\overline{\nu }_\mu \gamma `$, $`F_A=0.0058\pm 0.00008`$ is the pion axial vector form factor111Note that our definition of $`F_A`$ agrees with that of Ref.  but differs by a factor of $`1/2`$ from that given in Ref.  obtained from the radiative decays $`\pi ^{}l^{}\overline{\nu }_l\gamma `$ and $`r_\pi ^2=(0.439\pm 0.008)`$ fm2 is the pion charge radius squared. From now on these four chiral sum rules will be denoted by SR1, SR2, SR3 and SR4 respectively. It could be argued that as long that these chiral sum rules are taken in the chiral limit, the value of the pion decay constant should be the chiral limit value, $`f0.94f_\pi `$ . However, as long as the experimental data consists of real pions, we consider that it is more reasonable to use the real world value for the pion decay constant.

When switching quark masses on, only the first WSR remains valid while the second breaks down due to contributions from the difference of non-conserved vector and axial vector currents of order $`m_q^2/s`$ leading to a quadratic divergence of the integral. This is not numerically relevant in our analysis because we deal with finite energy sum rules, and in this case the contribution from non-zero quark masses is negligible.

The QCD vacuum condensates can be determined by virtue of the dispersion relation from another sum rule, that is, a convolution of the $`\rho _{VA}(s)`$ spectral function with an appropriate weight function. Let us define the operator product expansion of the chiral correlator in the following way

$$\mathrm{\Pi }(Q^2)|_{VA}=\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{Q^{2n+4}}C_{2n+4}(Q^2,\mu ^2)𝒪_{2n+4}(\mu ^2)\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{Q^{2n+4}}𝒪_{2n+4}.$$
(10)
The Wilson coefficients, including radiative corrections, are absorbed into the nonperturbative vacuum expectation values, to facilitate comparison with the current literature. The analytic structure of the $`\mathrm{\Pi }`$ is subject to the dispersion relation

$$\mathrm{\Pi }_{VA}(Q^2)=_0^{\mathrm{}}ds\frac{1}{s+Q^2}\frac{1}{\pi }\mathrm{Im}\mathrm{\Pi }_{VA}(s).$$
(11)
Condensates of arbitrary dimension are simply given by

$$𝒪_{2n+2}=(1)^n_0^{s_0}𝑑ss^n\frac{1}{2\pi ^2}\left(v_1(s)a_1(s)\right),n2,$$
(12)
which,if the asymptotic regime has reached, should be independent of the upper integration limit for large enough $`s_0`$. As can be seen from the experimental data, errors in the large $`s`$ region are very important, so large errors are expected in the evaluation of the condensates. The analysis of these sources of errors is one of our main goals in the present analysis, which will be obtained thanks to the natural capability of neural networks of smooth interpolating while implementing all the experimental information on errors and correlations.

### 2.1 Finite energy sum rules

As long as all previous integrals have to be cut at some finite energy $`s_0M_\tau ^2`$, since no experimental information on $`v_1(s)a_1(s)`$ is available above $`M_\tau ^2`$, we must perform a truncation that competes with all other sources of statistical and systematic errors, introducing a theoretical bias which is difficult to estimate. Many techniques have been developed to deal with this finite energy integrals, leading to the so-called Finite Energy Sum Rules (FESR). The paradigmatic example is the calculation of spectral moments , that is weighted integrals over spectral functions. Choosing appropriate weights allows to extract the maximum information possible from the experimental data while minimizing the contribution from the region with larger errors. This techniques allow a comparison of the same quantity evaluated on one side with experimental data and on the other side with theoretical input, basically the Operator Product Expansion with perturbative QCD corrections. The general expression that takes advantage of the analyticity properties of the chiral correlators is given by

$$_0^{s_0}𝑑sW(s)\mathrm{Im}\mathrm{\Pi }_{VA}^{(J)}(s)=\frac{1}{2i}_{|s|=s_0}𝑑sW(s)\mathrm{\Pi }_{VA}^{(J)}(s),$$
(13)
where $`W(s)`$ is an analytic function and $`s_0`$ is large enough for the OPE series to converge. The LHS of eq. (13) can be evaluated using the experimental input from spectral functions as determined in hadronic tau decays, while the RHS can be evaluated using the OPE representation of the chiral correlator. Finally, a fit is performed to extract the OPE parameters from the experimental data on spectral functions.

A common hypothesis in the majority of this kind of analysis is that the difference of the OPE representation for the chiral correlator from the full expression,

$$R[s_0,W]\frac{1}{2\pi i}_{|s|=s_0}𝑑s\left(\mathrm{\Pi }_{VA,OPE}(s)\mathrm{\Pi }_{VA}(s)\right)W(s),$$
(14)
can be neglected. This quantity is a measure of the OPE breakdown, also known as duality violation222 For a review of the current theoretical status of the quark-hadron duality violations, see Ref. . It is necessary to take into account that $`\mathrm{\Pi }_{VA,OPE}`$ fails at least in some region of the integration contour. This was shown in Ref. , where it was demonstrated that the OPE representation breaks down near the timelike real $`s`$ axis for insufficiently large $`s_0`$. The neglect of the duality violation component of the OPE is a key dynamical assumption and there exists several strategies to minimize its impact, as working with duality points  or using pinched Finite Energy Sum Rules , with polynomial weights that vanish at the upper integration limit. All these techniques yield different although compatible values for the $`𝒪_6`$, whereas non-compatible results are obtained for higher condensates.

Other types of finite energy sum rules have been used to extract the values of the condensates and other phenomenologically relevant related quantities. Borel sum rules and Gaussian sum rules  take advantage of certain combination of the condensates that theoretically optimize the accuracy of the extraction. Inverse moment sum rules  techniques make a connection between the phenomenological parameters of the QCD effective Lagrangian and the nonperturbative condensates. In section 5 we will compare our extraction of the condensates to those obtained with all these methods and argue why ours has a reasonable control of the different theoretical uncertainties.

Our approach is different with respect to previous determinations of nonperturbative condensates. First of all, we use the smooth interpolation of the neural network to extend the range of integration, so that our determination of the condensates corresponds to the asymptotic energy region $`s_0\mathrm{}`$. The second point is that the training method allows the incorporation of the chiral sum rules to the neural network parametrization of the spectral function, and therefore the specific weight used to determine the condensate turns out not to be relevant: different choices of weights differ by chiral sum rules that are already verified our parametrization. Therefore, the final results for the non-perturbative condensates that emerge from the neural network parametrization of the spectral function $`v_1(s)a_1(s)`$ are determined by

$$𝒪_{2n+2}=(1)^n_0^{\mathrm{}}𝑑ss^n\frac{1}{2\pi ^2}\left(v_1(s)a_1(s)\right)n2.$$
(15)
In the next sections will be argued why this choice is the most reasonable one, showing that all relevant constraints are verified.

## 3 Experimental data

Since the relevant spectral function for the determination of the condensates is the $`v_1(s)a_1(s)`$ spectral function, we need a simultaneous measurement of the vector and axial-vector spectral functions. Data from the ALEPH Collaboration ,  and from the OPAL collaboration  will be used, which provide a simultaneous determination of the vector and axial vector spectral functions in the same kinematic region and also provide the full set of correlated uncertainties for these measurements. Although the ALEPH data is of a higher quality due to the smaller errors, see fig. (1), the input from OPAL is complementary and will provide a cross check for our extractions of the nonperturbative condensates. There exists additional data on spectral functions coming from electron-positron annihilation, but their quality is lower than the data from hadronic tau decays and will be ignored here.

ALEPH experimental data consists on the invariant square-mass spectra for both the vector+axial vector and vector-axial vector components, that are related to the spectral functions by a kinematic factor and a branching ratio normalization

$$v_1(s)/a_1(s)\frac{M_\tau ^2}{6|V_{ud}|^2S_{EW}}\frac{B(\tau ^{}V^{}/A^{}\nu _\tau )}{B(\tau ^{}e^{}\overline{\nu }_e\nu _\tau )}\frac{dN_{V/A}}{N_{V/A}ds}\left[\left(1\frac{s}{M_\tau ^2}\right)\left(1\frac{2s}{M_\tau ^2}\right)\right]^1.$$
(16)
Altogether our parametrization is based on $`N_{dat}=`$ 61 experimental points for ALEPH, although the full experimental data consists in 70 points uniformly distributed between 0 and 3.5 GeV2, because only points with $`sM_\tau ^2=3.16`$ GeV2 are physically meaningful, and before this kinematic threshold is reached the invariant mass-squared spectrum vanishes due to phase-space suppression. For OPAL the data sample is a bit larger, $`N_{dat}=97`$, with the same restrictions as in the ALEPH case. Henceforth, $`\rho _{VA,i}^{(exp)}`$ will denote the $`i`$-th data point $`\rho _{VA}(s_i)v_1(s_i)a_1(s_i)`$. Figure (1) shows the experimental data used together with diagonal errors.

Note that errors are small in the low and middle $`s`$ regions and that they become larger as we approach the tau mass threshold. The last points are almost zero in the invariant mass spectrum, and are only enhanced in the spectral functions due to the large kinematic factor for $`s`$ near $`M_\tau ^2`$, so special care must be taken with the physical relevance of these points.

It is clear that the vanishing of the spectral function is not reached for $`sM_\tau ^2`$ and must be enforced artificially on the parametrization we are constructing, that is we must device a technique to impose the asymptotic constraint that at high $`s`$ this spectral function vanishes. The method we use takes advantage of the smooth, unbiased interpolation capability of the neural network: artificial points are added to the data set with adjusted errors in a region where $`s`$ is high enough that the $`\rho _{VA}(s)`$ spectral function should vanish. Once these artificial points are included, in a way to be discussed later, the neural network will smoothly interpolate between the real and artificial data points, also taking into account the constraints of the sum rules, as explained below.

The experimental data points are highly correlated, because the majority of the covariance matrix is composed of nonzero entries, so it is therefore crucial to take into account all their correlations, which are specially relevant in the high $`s`$ region. This is important because this region dominates the sum rule, eq. (12), that determines the vacuum condensates. As we shall discussed shortly, correlated errors are incorporated as a measure on the space of neural network parametrizations of the spectral functions using Monte Carlo statistical replicas of the experimental data.

## 4 Neural network parametrization

Ideally, a parametrization of spectral functions must incorporate all the information contained in the experimental measurements, i.e. their central values, their statistical and systematic errors and their correlations, furthermore, it must interpolate between them without introducing any bias. We will follow the method of Ref. , where an unbiased extraction of the probability measure in the space of structure functions of deep-inelastic scattering is performed, based on a coordinated use of Monte Carlo generation of data and neural network fits.

### 4.1 Probability measure in the space of spectral functions

The experimental data gives us a probability measure in an $`N_{dat}`$ dimensional space, assumed to be multigaussian. In order to extract from it a parametrization of the desired structure function we must turn this measure into a measure $`𝒫\left[\rho _{VA}\right]`$ in a space of functions. Once such a measure is constructed, the expectation value of any observable $`\left[\rho _{VA}(s)\right]`$ can be found by computing the weighted average

$$\left[\rho _{VA}(s)\right]=𝒟\rho _{VA}\left[\rho _{VA}(s)\right]𝒫\left[\rho _{VA}\right].$$
(17)
Errors and correlations can also be obtained from this measure, by considering higher moments of the same observable with respect to the probability distribution.

The determination of an infinite-dimensional measure from a finite set of data points is an ill-posed problem, unless one introduces further assumptions. In the approach of Ref. , neural networks are used as interpolating functions, so that the only assumption is the smoothness of the spectral function. Neural networks can fit any continuous function through a suitable training; smoother functions require a shorter training and less complex networks. Hence, an ideal degree of smoothness can be established on the basis of a purely statistical criterion without the need for further assumptions.

### 4.2 Fitting strategy

The construction of the probability measure is done in two steps: first, a set of Monte Carlo replicas of the original data is generated. This gives a representation of the probability density $`𝒫\left[\rho \right]`$ at points $`(s_i)`$ where data is available. Then a neural network is fitted to each replica. The ensemble of neural networks gives a representation of the probability density for all $`s`$: when interpolating between data the uncertainty will be kept under control by the smoothness constraint, but it will become increasingly more sizable when extrapolating away from the data region.

The $`k=1,\mathrm{},N_{rep}`$ replicas of the data are generated as

$$\rho _{VA,i}^{(art)(k)}=\rho _{VA,i}^{(exp)}+r_i^{(k)}\sigma _i,$$
(18)
where $`\rho _{VA,i}^{(exp)}=\rho _{VA}(s_i)`$ are the original data, $`\sigma _i`$ is the diagonal error, and $`r_i^{(k)}`$ are univariate gaussian random numbers whose correlation matrix equals that of the experimental data. The fact that the correlation matrix of the $`r_i^{(k)}`$ equals that of the experimental data is crucial to retain all the experimental information in our treatment. Then a set of $`N_{rep}=1000`$ replicas of this form is generated, and is verified that the central values, errors and correlations of the original experimental data are well reproduced by taking the relevant averages over a sample of this size. A explained above, the asymptotic constraint that $`\rho _{VA}(s\mathrm{})=0`$ has been implemented by adding a number of artificial data points with adjusted errors.

To verify that the central values, errors and correlations of the original experimental data are well reproduced, we can define statistical estimators that measure the deviations from the original correlations. A suitable one is the scatter correlation, which measure the deviations of the averages over the replica set from the original experimental values, and are defined as follows for the central value

$$r[\rho _{VA}^{(art)}]=\frac{\rho _{VA}^{(exp)}\rho _{VA}^{(art)}_{rep}_{dat}\rho _{VA}^{(exp)}_{dat}\rho _{VA}^{(art)}_{rep}_{dat}}{\sigma _s^{(art)}\sigma _s^{(exp)}},$$
(19)
where the scatter variances are defined by

$$\sigma _s^{(exp)}=\sqrt{\left(\rho _{VA}^{(exp)}\right)^2_{dat}\left(\rho _{VA}^{(exp)}_{rep}\right)^2},$$
(20)
$$\sigma _s^{(art)}=\sqrt{\left(\rho _{VA}^{(exp)}_{rep}\right)^2_{dat}\left(\rho _{VA}^{(art)}_{rep}_{dat}\right)^2},$$
(21)
and similarly for the diagonal errors and the correlations. In table 1 we show the scatter correlations for the central values, the errors and the correlations. We observe that we need $`N_{rep}=100`$ replicas to maintain the correlations of the original data, which is the main purpose of our analysis. We have checked that increasing the number of training replicas does not decrease the errors in the extraction of the condensate further, meaning that we have reached a faithful representation of errors.

Each set of generated data is fitted by an individual neural network. A neural network , is a function of a number of parameters, which fix the strength of the coupling between neurons and the threshold of activation of each neuron. The architecture of the network has been chosen to be 1-4-4-1, small enough to avoid overlearning and large enough to capture the non-linear structure of experimental data. The networks that we use are multilayer feed-forward neural networks constructed according to the following recursive relation

$$\xi _i^{(l)}=g\left(h_i^{(l)}\right),$$
(22)
$$h_i^{(l)}=\underset{j=1}{\overset{n_{l1}}{}}\omega _{ij}^{(l)}\xi _j^{(l1)}\theta _i^{(l)},$$
(23)
where $`\omega _{ij}^{(l)}`$ is the weight, the strength of the connection between two neurons, $`\theta _i^{(l)}`$ are the thresholds of each neuron, $`\xi _i^{(l)}`$ is the activation state of each neuron and $`g`$ is the activation function of the neurons.

We divide the training of the neural network in two epochs. In a first epoch, the training method is done by backpropagation, where the parameters of the network are fitted by minimizing the error function:

$$E_{\mathrm{err}}^{(k)}=\underset{i=1}{\overset{N_{dat}}{}}\frac{\left(\rho _{VA,i}^{(art)(k)}\rho _{VA,i}^{(net)(k)}\right)^2}{\sigma _i^{(exp)^2}},$$
(24)
where $`\rho _i^{(net)(k)}`$ is the prediction of the $`i`$th data point from the net trained on the $`k`$th replica of the data. A more detailed review of neural networks learning techniques is presented in the appendix A.

In a second training epoch, a different training technique called genetic algorithms training is used to implement the constraints from the sum rules. As explained below and in the appendix A, this technique allows us to implement in our training non-local constraints, as convolutions of the neural network output, with adjusted weights so that the chiral sum rules control the neural network interpolation in the data region where errors are greater. The error function eq. (24) is modified by adding a contribution proportional to the difference of the chiral sum rules, evaluated with the output of the trained networks, and their theoretical values, that is

$$E_{\mathrm{tot}}=E_{\mathrm{err}}+E_{\mathrm{sr}}=E_{\mathrm{err}}+\underset{i=\mathrm{sr}_1}{\overset{\mathrm{sr}_4}{}}w_{\mathrm{sr}_i}\left(_0^{s_0}dsf_i\left(\rho _{VA}(s)\right)A_i\right)^2,$$
(25)
where $`w_{\mathrm{sr}_i}`$ is the relative weight of each sum rule and $`A_i`$ is the theoretical value of the corresponding sum rule, eqns. (6-9). We note that this definition introduces a new set of, in principle, arbitrary parameters, that is the relative weights of the sum rules. As explained below, these are determined by stability criteria, demanding that the contribution from $`E_{\mathrm{sr}}`$ is similar to that of $`E_{\mathrm{err}}`$ and that the final result is not sensitive to the specific values of these parameters.

The basic idea of the genetic algorithms training, also known as natural selection training, works as follows. The training is divided in generations. For each generation the parameters of the network (weights and thresholds) are arranged to form a chain, called the ADN chain. This chain is replicated many times, creating a population of identical individuals. Later, random mutations are applied to each individual, where by mutation we mean a small change in one of the bits of his ADN chain. Then, the error function associated with each mutated individual is computed, which implies passing back the ADN bits to their original status of weights and thresholds and calculating the output of this new network. Only the best individuals are kept while discarding the rest, mimicking natural selection. This method provides a suitable technique to implement the effect of the chiral sum rules on our neural network training. Note that this technique leads to an important increase on the computing time, due to the fact that the chiral sum rules must be numerically evaluated many times each generation.

The main advantage of genetic algorithms is that they allow neural networks to learn from error functions that may be as complicate as making impossible the use of backpropagation training. Furthermore, genetic algorithms can be proven to search efficiently the parameter space of solutions, exploring exponentially many more times reasonable outputs as compare to manifestly wrong ones. Genetic algorithms can also handle the training of very large neural networks.

The parametrization obtained by means of the genetic algorithm training is represented in figure (2) where the output of the trained neural network, without and with the inclusion of the chiral sum rules, is compared with the experimental data points together with the corresponding statistical errors. It is clear that the effect of the chiral sum rules on the trained neural network output is forcing it to reach faster the asymptotic behavior of the spectral function $`\rho _{VA}(s)`$.

A common problem in the genetic algorithm learning techniques is getting stuck in a local minimum of the error function, far from the absolute minimum. In our training this difficulty has been bypassed by means of different simple modifications of the basic training procedure. First, within each generation large additional mutations are performed that allow the network configuration to escape from local minima. Secondly, as the training advances, the rate of the mutations decreases, allowing for way a better local learning. These modifications are instrumental to decrease the large duration of the training.

### 4.3 Results and validation

Once all the parameters of the training have been determined by stability criteria, an independent set of neural networks is trained on the spectral functions $`v_1(s)a_1(s)`$. The length of the training is fixed by studying the behavior of the error function $`E^{(0)}`$, as defined in eq. (24) for the neural net fitted to the central experimental values, and asking that $`E^{(0)}/N_{dat}`$ stabilizes to a value close to one, which can be considered a good training.

A number of checks is then performed in order to make sure that an unbiased representation of the probability density has been obtained. First, we have verified that the covariance of two data points computed from the Monte Carlo sample of nets is on average very close to the corresponding covariance matrix element of the data. Since correlations of the data are entirely due to systematic errors, this indicates that these errors are correctly reproduced. Statistical estimators are then constructed as in eq. (19) but now referred to the trained neural networks over the replicas, to explicitly verify that the training maintains all the experimental information. The scatter correlation is now defined as

$$r[\rho _{VA}^{(net)}]=\frac{\rho _{VA}^{(exp)}\rho _{VA}^{(net)}_{rep}_{dat}\rho _{VA}^{(exp)}_{dat}\rho _{VA}^{(net)}_{rep}_{dat}}{\sigma _s^{(net)}\sigma _s^{(exp)}},$$
(26)
and the corresponding values for the training of $`N_{rep}=100`$ replicas are presented in table 2. It is seen that the central values and the correlations are well reproduced, whereas this is not the case for the diagonal errors.

The average standard deviation for each data point computed from the Monte Carlo sample of nets is substantially smaller than the experimental error. This is due to the fact that the network is combining the information from different data points by capturing and underlying law, or that it is introducing a smoothing bias. This effect is enhanced by the inclusion of sum rules constraints. All networks have to fulfill these constraints which forces the fit to behave smoothly in a region where experimental data are very large. This should be understood as a success of the fitting procedure.

The final set of neural networks $`\rho _i^{(net)(k)}`$ provides a representation of the probability measure in the space of structure functions, which can be used to estimate any functional average, defined as in. eq (17) using

$$\left[\rho _{VA}(s)\right]=\frac{1}{N_{rep}}\underset{i=1}{\overset{N_{rep}}{}}\left[\rho _{VA}^{(net)(k)}(s)\right].$$
(27)
In particular, the average and standard deviation of the nonperturbative condensates computed using the Monte Carlo sample will provide a determination of the central values and errors of these condensates.

### 4.4 Details of the genetic algorithm training

As explained above, in the second part of the training the chiral sum rules eqs. (6-9) are incorporated to the error functional, eq. (24). These sum rules act as constraints on the neural network output, that is, the main contribution to the error function (which determines the learning of the network) still comes from the diagonal errors, and the sum rules are only relevant in the region where the errors are larger. The relative weights of the chiral sum rules will be chosen according to a stability analysis. The effect of including sum rules in the learning procedure is responsible for enforcing the desire vanishing oof the $`\rho _{VA}(s)`$ spectral function, which is badly needed for a reliable extraction of the nonperturbative vacuum condensates.

Obtaining maximum stability in our output is crucial for a proper parametrization of the spectral function. In the case of the relative weights of the chiral sum rules, we train the same network with different relative weights and search for the region where both contributions to the error function, the contribution from the errors and the contribution from the sum rules, are comparable. This is shown in figure (3). The observed behavior is not surprising: for large relative weights the contribution from the sum rule is small because the training forces its learning but, as a consequence, the contribution from the experimental data increases. This behavior can be observed explicitely if we plot the evolution of the sum rules, as calculated with the network output of the trained network, as a function of their relative weights. We observe in figure (4) that, as the relative weight increases, the network output better verifies the corresponding chiral sum rule.

Genetic algorithms thus allow to implement additional constraints in the training of neural networks in a smooth and efficient way. Its main drawback is the increase of the training time, because it is a random rather than a deterministic learning technique. In figure (5) we represent an example of the training of one replica. It can be observed a sharp transition when the sum rules constraints are introduced, but later the training forces the error function to stabilize to a situation similar to the initial training epoch. This sudden jump of the error function can be understood as follows: when the sum rules constraints are introduced, the training tends to verify them, causing that the net output does not follow the experimental values. Nevertheless, as generations go on, the net output begins to recover the original situation, while maintaining the verification of the sum rules constraints. When the number of generations is large, the error function approaches a value close to one, as is needed to keep systematic errors under control.

A key issue in this procedure is to guarantee stability of results with respect to the relative weights of the chiral sum rules. In our training normalized sum rules are used, that is, if $`A_j`$ is the theoretical value of the $`j`$th chiral sum rule, the corresponding contribution to the error function will be

$$w_{\mathrm{sr}_j}\left(_0^{s_0}dsf_j\left(\rho _{VA}(s)\right)A_j\right)^2/\left(A_j\right)^2,$$
(28)
therefore we expect the relative weights in the stability region of order 1. The only exception is the second Weinberg sum rule, eq. (8) whose relative weight has to be determined demanding stability of the network training, and that turns out to be around $`w_{\mathrm{sr}_3}=10^1`$.

Let us emphasize that the two Weinberg chiral sum rules are well verified by our neural network parametrization, and thus have been incorporated to the information contained on the experimental data. This fact will be crucial later because different extraction methods, differing in combinations of these chiral sum rules, can be shown to be equivalent in the asymptotic region $`s_0\mathrm{}`$. In figure (6) the two Weinberg sum rules, eqs. (7,8) evaluated with the neural network parametrization of the spectral function $`\rho _{VA}(s)`$ are represented. Both chiral sum rules are well verified in the asymptotic region, beyond the range of available experimental data. This also will ensure the stability of the evaluation of the condensates with respect to the specific value of $`s_0`$ chosen as long as it stays in the asymptotic region.

## 5 Determination of the nonperturbative condensates

Using the neural parametrization of spectral functions, we can compute for each trained replica any given sum rule. Because the neural parametrization retains all the experimental information (it even allows for a determination of errors and correlations), we can view values coming from the neural networks as direct experimental determinations of convolutions of the spectral function $`\rho _{VA}(s)`$. The value of the condensates $`𝒪_6,𝒪_8`$ and higher dimensional condensates is then extracted from the value of an appropriate sum rule, eq. (12). The method we will follow is the evaluation of the vacuum condensates as a function of the upper limit of integration for each replica and compute the mean and standard deviation. As has been explained before, it is crucial to represent the value of the different sum rules as a function of the upper limit of integration, to check both its convergence and its stability.

Our method works as follows. First, we train a neural network on each replica. On a first training epoch we do not use the sum rules, so that the training can arrive to the best possible minimum. This is important when training neural networks because when further constraints to the training are added, as in our case the chiral sum rules, the goodness of the fit will be better if we start from a deep local minimum. On a second training epoch, we add to the fitness the contribution from the sum rules, where the relative weights are chosen so that the sum rules never represent more than the contribution from the experimental errors to the total fitness. Then, sum rules act as a smooth constraint on the network training, being more relevant in the regions with larger errors and thus enforcing the asymptotic vanishing of the spectral function. This technique prevents the contribution of the chiral sum rules to become so strong that overcome the experimental data with the corresponding errors.

### 5.1 Central values

The first criterion to judge the reliability of a QCD sum rule is its independence, at large values of $`s`$ from the value of the upper integration limit, that its, its saturation. We then need to explore the values for the final condensates which are stable against the limit of integration of the sum rule. This stability criterium is completed with demanding independence of the results on the specific polynomial entering the sum rule. Further criteria are stability with respect to the precise artificial endpoints added to the data and with respect to the relative weights in the error function used to train the neural networks.

Stable results are obtained for the dimension six condensate $`𝒪_6`$ whereas higher condensates e. g. $`𝒪_8`$ are less stable. Fig. 7 shows the outcome for $`𝒪_6`$ and $`𝒪_8`$ including the propagation of statistical errors. It is clearly seen that convergence in the limit of integration $`s_0`$ is obtained due to the addition of sum rules and endpoints in the learning procedure. The central values for the condensates in the asymptotic limit come out to be: $`s\mathrm{}`$:

$`𝒪_6=4.210^3\mathrm{GeV}^6,`$
$`𝒪_8=12.710^3\mathrm{GeV}^8.`$ (29)
The value of the $`𝒪_6`$ is a cross check of the validity of our treatment: not only there are strong theoretical arguments that support the fact that $`𝒪_6`$ is negative , but also all previous determinations with different techniques yield negative results, being the majority of them compatible with ours within errors.

We note that our evaluation of the condensates is compatible with some of our previous evaluations and has a similar error. This is though misleading as the error quoted here is only statistical and a discussion on systematic errors is needed (and done below). We can also obtain values for the higher dimensional nonperturbative condensates:

$`𝒪_{10}=7.810^2\mathrm{GeV}^{10},`$
$`𝒪_{12}=2.610^1\mathrm{GeV}^{12}.`$ (30)
Although stability deteriorates as compared to the case of the lower dimensional condensates, these central values for the condensates are alternated in sign.

### 5.2 Discussion of errors

The discussion of the various sources of errors is crucial to our treatment. We enumerate and discuss then in turn:

1. Statistical error propagation from the experimental covariance matrix.
This is the best understood and treated error source in our analysis. As explained above, the neural network parametrization defines an unbiased probability measure in the space of spectral functions that provides a nonlinear error propagation. This source of error is kept under control by using the averages over Monte Carlo replica. The band for central values of the condensates allowed by this error propagation can be visualized in fig. (7). Numerically, the contribution to the experimental error (statistics and correlations) to the central values is

$`𝒪_6=\left(4.2\pm 1.1\right)10^3\mathrm{GeV}^6,`$
$`𝒪_8=\left(12.7\pm 6.4\right)10^3\mathrm{GeV}^8,`$
$`𝒪_{10}=\left(7.8\pm 2.4\right)10^2\mathrm{GeV}^{10},`$
$`𝒪_{12}=\left(2.6\pm 0.8\right)10^1\mathrm{GeV}^{12}.`$ (31)
Note that the sign obtained for each condensate remains unaltered within error.
2. Choice of the polynomial in the finite energy sum rule.
In principle, there are potentially important systematic uncertainities coming from the method of extraction of the condensates. These are much more difficult to assess, as noted when looking through the extense available literature. The extraction of $`𝒪_6`$ turns out to be clean and its errors are essentially of statistical nature. The uncertainty increases with the dimension of the condensate. Let us elaborate further these statements.

Consider the following convolutions

$$Z_{6a}_0^{s_0}𝑑s\frac{1}{2\pi }s^2\rho _{VA}(s)$$
(32)
$$Z_{6b}s_0^2_0^{s_0}𝑑s\frac{1}{2\pi }\left(1\frac{s}{s_0}\right)^2\rho _{VA}(s)f_\pi ^2s_0^2$$
(33)
The second equation is only equivalent to the first if, for some $`s_0`$, both Weinberg sum rules are satisfied. Although experimental data on tau decays do not exactly saturate these sum rules, the neural network parametrization trained to obey all the sum rules showed that Weinberg sum rules can indeed be well verified in the asymptotic region. We should then expect that $`Z_{6a}`$ and

$$Z_{6b}_0^{s_0}𝑑s\frac{1}{2\pi }\left(s^22ss_0\right)\rho _{VA}(s)$$
(34)
yields similar results for $`𝒪_6`$ within errors, as can be seen in fig. (8). The same applies to the dimension 8 condensate $`𝒪_8`$, where we now define the following finite energy sum rules:

$$Z_{8a}_0^{s_0}𝑑s\frac{1}{2\pi }s^3\rho _{VA}(s),$$
(35)
$$Z_{8b}_0^{s_0}𝑑s\frac{1}{2\pi }\left(s^3ss_0^2\right)\rho _{VA}(s),$$
(36)
$$Z_{8c}_0^{s_0}𝑑s\frac{1}{2\pi }\left(s^3+3s_0^2s\right)\rho _{VA}(s).$$
(37)
We conclude that the neural network parametrization of spectral functions properly trained to accomodate for all sum rules provides estimates for condensates which are independent of the choice of a specific finite energy sum rule.
3. Dependence on the endpoints.
Our fit implements the asymptotic constraint that $`\rho _{VA}(s\mathrm{})=0`$ by adding artificial endpoints. It is then necessary to verify the degree of sensitivity of our ouput to the precise location of these endpoints. As shown if Fig. (9), the sign of the dimension six and eight condensates remains unaltered when endpoints range between 3.5 and 4 GeV2. We also observe relatively large, although compatible within errors, fluctuations of the central values. This effect may be related to the presence of small wiggles in the spectral function $`\rho _{VA}(s)`$ for large $`s`$. The contribution of this source of uncertainty to $`𝒪_8`$ turns to dominate over statistical uncertainties and can be estimated to be

$$2.510^2𝒪_8510^3$$
(38)
while for $`𝒪_6`$ is comparable to the uncertainty due to the statistical errors of the experimental data, that is,

$$610^3𝒪_6210^3.$$
(39)
4. Chiral sum rules.
This turns out to be the main source of systematic errors for the dimension 6 condensate. Chiral sum rules can be forced to be fulfilled by adequate training to any desired degree of precision. This, though, introduces a large increase in the total error function, eq. (25), coming from the experimental error piece. It is then necessary to make an appropriate choice of relative weights between the error associated to experimental data and the error associated to the fulfillment of chiral sum rules.

We may advocate that the most appropiate relative weights for the normalized chiral sum rules are $`𝒪(1)`$. This is due to the fact that the total error function jumps above 1 when too large relative weights are considered, as seen in figs. (3, 4). We have thus performed a multi-dimensional stability analysis searching for the relative weights that produce a minimum sensitive final result, supplemented with the condition that in any case the contribution from the experimental errors to the error function can be greater than 1. The most suitable relative weights for the chiral sum rules turn out to be

$$w_{sr1}=1.0w_{sr2}=510^2w_{sr3}=0.3w_{sr4}=110^2.$$
(40)
This stability analysis shows that the sign of the central values of the condensates is not very sensitive to the relative weights for the chiral sum rules. The estimation of the error associated with this uncertainty leads to the following range of values for the lowest dimensional condensate

$`210^3𝒪_6610^3\mathrm{GeV}^6.`$ (41)
For $`𝒪_8`$ the statistical errors and the systematic error due to the position of the artificial endpoints turn out to dominate over this source of uncertainty. Similar estimates for the condensates of higher dimensions turn out not to be reliable, and therefore we present only the central values obtained in this analysis together with the statistical errors.

### 5.3 Analysis of the $`s_0=1.5`$ GeV2 duality point

Some values of previous extraction of the condensates  are based on the existence of a duality point around $`s_01.5`$ GeV2. Our neural network parametrization is such that the second Weinberg sum rules is indeed verified around this point. Consequently, the values of the condensates computed by truncating different finite energy sum rules at this point do agree, as can be seen in fig. (6). Nevertheless, as shown in fig. (8), the value of the condensates at this duality point is different than the asymptotic one. Both results share the same sign but not the same absolute magnitude for $`𝒪_6`$, and with the opposite sign for $`𝒪_8`$ (which is precisely the results that the authors in Ref.  obtain).

These extractions using the first duality point can be justified by some resonance models of the hadronic spectrum , inspired in the large$`N_C`$ limit of QCD with the additinal assumption of the validity of the Minimal Hadronic Approximation (MHA), in which the spectral functions are saturated by the pion pole, the first axial vector and the first rhoo vector resonances. Experimentally, it is observed that at different duality points, even as low as $`1.5`$ GeV2, there is a local quark-hadron duality, meaning that the OPE at the quark- gluon level and that evaluated with the entire resonance hadronic spectrum coincide. Whether this apparent duality point at $`s_0=1.5`$ GeV2 is an accident only due to the fullfillment of the second Weinberg sum rule, or it is really a consequence of the full QCD hadronic spectrum remains to be understood. What it is clear from our analysis is that the condensates evaluated at this first duality point are different to those obtained in the asymptotic regime $`s_0\mathrm{}`$, where the validity of the OPE is less questioned.

### 5.4 Spectral functions and the electroweak penguin $`Q_7`$

As a byproduct of our analysis additional sum rules of the spectral function $`\rho _{VA}`$ which are relevant to phenomenology can be estimated. As an example333In this section the work of Ref.  is followed, the reader is directed to this reference for definitions and notation, the sum rule

$$𝒜_{LR}(\mu ^2)_0^{s_0}𝑑ss^2\mathrm{ln}\left(\frac{s}{\mu ^2}\right)\frac{1}{2\pi ^2}\rho _{VA}(s),$$
(42)
will be considered, where $`\mu ^2`$ is a arbitrary factorization scale that cancels in the computation of physical observables. In this analysis the value $`\mu ^2=2`$ GeV2 will be used. Eq. (42) is relevant to the evaluation of $`\mathrm{Im}G_E`$, where $`G_E`$ is one of the couplings of the low energy chiral Lagrangian describing $`|\mathrm{\Delta }S|=1`$ transitions . The importance of a precise determination of this coupling relies on the fact that $`\mathrm{Im}G_E`$ is one of the most important contributions to $`ϵ^{}`$ in the Standard Model. In turn, its value is dominated by the electroweak penguin contributions $`Q_7`$ and $`Q_8`$, which explains why the data on spectral functions from hadronic tau decays is important in its determination.

Following the same steps that lead to the determination of the vacuum condensates, the same procedure for the sum rule eq. (42) is repeated. The result that is obtained in the asymptotic limit $`s_0\mathrm{}`$ is

$$𝒜_{LR}=\left(6.9\pm 1.6\right)10^3\mathrm{GeV}^6,$$
(43)
as can be seen in figure (10). It should be noted that the present determination is in good agreement with that obtained in the original work444Note however that a different normalization for the spectral correator is used. . The quoted error only refers to the propagation of experimental uncertainties.

### 5.5 Results and comparison with other determinations

Our determination of the nonperturbative condensates including the statistical error coming from the experimental errors and correlations was given in eq. (31). However, the error which dominates the determination of $`𝒪_6`$ comes from the relative weights of the chiral sum rules to be obeyed. We have performed a stability analysis on these relative weights that produces a final result:

$$𝒪_6=\left(4.0\pm 2.0\right)10^3\mathrm{GeV}^6.$$

As explained above, for the dimension 8 condensate, the systematic error associated with the endpoint position is comparable to the statistical uncertainty, that combine to yield a value

$$𝒪_8=\left(12\genfrac{}{}{0pt}{}{+7}{11}\right)10^3\mathrm{GeV}^8.$$
(44)
For higher dimensional condensates it is much difficult to estimate the different sources of systematic uncertainties. We, then, quote the central values we obtained and their statistical error:

$`𝒪_{10}=\left(7.8\pm 2.4\right)10^2\mathrm{GeV}^{10},`$
$`𝒪_{12}=\left(2.6\pm 0.8\right)10^1\mathrm{GeV}^{12}.`$ (45)
A similar analysis has been performed with the OPAL data, yielding equivalent results but with larger errors, due to the larger statistical uncertainties as compared to the ALEPH experimental data. The values of the QCD nonperturbative condensates have been previously extracted from the ALEPH and OPAL data, with different techniques and different results, as summarized in table 3.

Note that our results agree, at least on the sign, with that of Refs. ,,. This is also true for the higher dimensional condensates of Ref. , where the authors obtain:

$`𝒪_{10}=\left(4.8\pm 1.0\right)10^2\mathrm{GeV}^{10},`$
$`𝒪_{12}=\left(1.6\pm 0.26\right)10^1\mathrm{GeV}^{12}.`$ (46)
in agreement with eq. (31) , although the errors in our determination are somewhat larger.

There are some differences between these previous determinations and the present one, eq. (45). The first one is that we do not make any assumption on the values of the higher dimensional nonperturbative condensates. In many analysis the effect of $`𝒪_D`$ for $`D10`$ is simply neglected to get closed expressions for the condensates. In our analysis, though, we do not need to make this hypothesis. Neither we need to assume that the chiral sum rules are verified, because the chiral sum rules enter as constraints in the genetic algorithm training (see fig. (6)). This is relevant because previous analysis showed that the chiral sum rules are not verified for $`s_0=M_\tau ^2`$, except the DMO sum rule, implying that one must be extremely careful when dealing with them. A second main difference is the absence of theoretical bias introduced in other analysis with a choice of a given finite energy sum rules. Moreover, the smooth interpolating capability of the neural network lets the integration range to be taken up to arbitrarily high energies.

## 6 Conclusions

We have presented a determination of the nonperturbative vacuum condensates $`𝒪_6`$ and $`𝒪_8`$ from the spectral functions from hadronic tau decays aimed at minimizing the sources of theoretical bias which might be cause of concern in existing determinations of these condensates from spectral functions. This determination is based on a bias-free neural network parametrization of the $`v_1(s)a_1(s)`$ spectral function, inferred from the data, which retains all the information on experimental errors and correlations, and supplemented with the additional theoretical input of the chiral sum rules.

Our final results give negative central values for the dimension 6 and 8 condensates. These results take into account the propagation of statistical errors and their correlations. Morevover, the main source of systematic error in our procedure is identified as the choice of relative weights assigned to chiral sum rules in the fitness function used to train the neural networks. In the case of the dimension 6 condensate a stability analysis can be performed. Higher dimension condensates carry larger errors, although the sign of the condensates seem to remain unaltered.

The sign of the dimension eight condensate $`𝒪_8`$ deserves further comments. Our central value is negative within statistical errors but is sensitive to the position of artificial endpoints added to enforce the vanishing of the spectral function for large values of $`s`$. This produces a systematic bias as possible wiggles of the spectral function around the asymptotic zero value are suppreseed. Those wiggles may indeed produce a change of sign of $`𝒪_8`$. This is not the case in our approach as the smoothness of the neural network tends to avoid such wiggles, which might lead to a systematic error. Although our results seem to point in the same direction of other recent previous extractions of $`𝒪_8`$ we consider the issue of the sign of the vacuum condensate $`𝒪_8`$ to remain open.

Another result of this work is the implementation of a technique based on genetic algorithm neural network training, which extends the capabilities of neural network data analysis allowing to incorporate non-local constraints like convolutions in the training. This technique extends previous efforts  oriented to the improvement of high energy physics data analysis, specially for the strong interaction sector.

Acknowledgments

It is a pleasure to thank S. Peris and E. de Rafael for suggesting the application of Ref.  to this problem and for fruitful discussions, and J. Prades and A. Pich for useful comments and suggestions. We would also like to thank A. Hoecker and W. Mader for their help with the experimental data of the ALEPH and OPAL collaborations respectively. This work was supported by the projects MCYT FPA2001-3598, GC2001SGR-00065 and by the Spanish grant AP2002-2415.

## Appendix A Neural network techniques

In this section a brief review of the neural network training algorithms that have been used in the this analysis is presented. Two different learning algorithms have been used for the neural network training: in a first epoch the only contribution to the error function comes from the statistical experimental errors, and the used learning algorithm is known as backpropagation. Then, in a second training epoch, the contribution from the chiral sum rules is added to the error function. As long as convolutions over the neural network output are non local constraints, the previous technique is no longer useful and another learning algorithm is needed, which is called genetic algorithms training. Now each of this techniques will be introduced, directing the interested reader to the standard reviews and textbooks  on neural networks and applications.

Learning by backpropagation allows to train multilayer neural networks and has proved to be an excellent tool in classification, interpolation and prediction tasks. It is a standard technique that has been recently applied to data analysis in high energy physics, see Ref. . The starting point is a set of input-output patters,

$$(𝐱^\mu ,𝐳)R^n\times R^m,\mu =1,\mathrm{},p,$$
(47)
that network must learn. In our case, each input-output pattern consists of a single data point, the input being the energy $`s`$ and the output the spectral function $`\rho _{VA}(s)`$.

The basis of the network learning is the error function, also known as fitness functional. The error function is defined as the difference between the actual and the desired output of the net, measured over the training set, and weighted with the experimental errors. It is given by

$$E=\underset{i=1}{\overset{N_{dat}}{}}\frac{\left(\rho _i^{(exp)}\rho _i^{(net)}\right)^2}{\sigma _i^{(exp)^2}}.$$
(48)
Applying the gradient descent minimization procedure, that is, looking for the direction of steepest descent of the error function, the appropriate changes in the network parameters such that the error function decreases can be determined. The error is introduced in the units of the last layer by

$$\mathrm{\Delta }_i^{(L),\mu }=g^{}(h_i^{(L)\mu })\left[o_i(\stackrel{}{x}^\mu )z_i^\mu \right],$$
(49)
and then backpropagated to the rest of the network

$$\mathrm{\Delta }_j^{(l1),\mu }=g^{}(h_j^{(l1)\mu })\underset{i=1}{\overset{n_l}{}}\mathrm{\Delta }_i^{(l),\mu }\omega _{ij}^{(l)},$$
(50)
and the last step consistes on the update of the weights and thresholds of the net

$$\delta \omega _{ij}^{(l)}=\eta \underset{\mu =1}{\overset{p}{}}\mathrm{\Delta }_i^{(l)\mu }\xi _j^{(l1)\mu }+\alpha \delta \omega _{ij}^{(l)}(\mathrm{last}),$$
(51)
$$\delta \theta _i^{(l)}=\eta \underset{\mu =1}{\overset{p}{}}\mathrm{\Delta }_i^{(l)\mu }+\alpha \delta \theta _i^{(l)}(\mathrm{last}),$$
(52)
where $`\eta `$ is the learning rate parameter which controls the velocity of the training and the term with $`\alpha `$ is called a momentum term, which improves the algorithm so that the training does not get stuck in a local minimum. The main advantage of this method is that it is deterministic and it has been used repeatedly in different situations, always successfully.

As stated above, the second of the neural networks techniques that are used in our analysis is called genetic algorithm learning, also known as natural selection learning. As long as the chiral sum rules are convolutions of the network output, that is are non local functions on the error function, meaning that they depend not on a single network output but on the global output, the usual backpropagation training techniques are not useful in this second epoch of our training. Now the error functional has the form of eq. (48) but with additional convolutions of the neural network output

$$E=\underset{i=1}{\overset{N_{dat}}{}}\frac{\left(\rho _i^{(exp)}\rho _i^{(net)}\right)^2}{\sigma _i^{(exp)^2}}+\underset{i=1}{\overset{N_{conv}}{}}w_i\left(_0^{s_0}𝑑sf_i\left(\rho ^{(net)}(s)\right)A_i\right),$$
(53)
where $`A_i`$ is the theoretical value of the $`i`$th sum rule and $`w_i`$ is the relative weight of this sum rule. In this case genetic algorithms are used, a training method inspired in the evolutionary theories in biology. In this method, the network parameters are transformed into bits in an ADN chain. These chains are replicated, and then some bits are mutated with a certain probability, and only those chains with the smallest contribution to the error functional survive. By analogy, this method is also known as natural selection learning. A simple scheme of the recursive process can be seen as follows: the starting piont is the set of parameters that define the neural network,

$$\omega _{11}^{(1)},\omega _{11}^{(2)},\mathrm{},\theta _1^{(1)},\theta _1^{(2)},\mathrm{}$$

$$$$

Creation of ADN chain

$$$$

$$ADN_1=(\omega _{11}^{(1)},\omega _{11}^{(2)},\mathrm{},\theta _1^{(1)},\theta _1^{(2)},\mathrm{})$$

$$$$

Replication of ADN chains

$$$$

$$ADN_1=(\omega _{11}^{(1)},\omega _{11}^{(2)},\mathrm{},\theta _1^{(1)},\theta _1^{(2)},\mathrm{}),ADN_2=(\omega _{11}^{(1)},\omega _{11}^{(2)},\mathrm{},\theta _1^{(1)},\theta _1^{(2)},\mathrm{})\mathrm{}$$

$$$$

Random mutation of bits in the ADN chains

$$$$

$$ADN_1=(\omega _{11}^{(1)},\omega _{11}^{(2)},\mathrm{},\theta _1^{(1)},\theta _1^{(2)},\mathrm{}),ADN_2=(\omega _{11}^{(1)}+\delta _2\omega _{11}^{(1)},\mathrm{},\theta _1^{(1)}+\delta _2\theta _1^{(1)},\mathrm{})\mathrm{}$$

$$$$

Selection of best ADN chain

$$$$

$$ADN_j=(\omega _{11}^{(1)}+\delta _j\omega _{11}^{(1)},\mathrm{},\theta _1^{(1)}+\delta _j\theta _1^{(1)},\mathrm{})$$

$$$$

New weights and thresholds and decrease of fitness

$$$$

$$\omega _{11}^{(1)}+\delta _j\omega _{11}^{(1)},\mathrm{},\theta _1^{(1)}+\delta _j\theta _1^{(1)},\mathrm{}$$

This is a very simple genetic algorithm, more modifications could be added to improve its efficiency like crossing between individuals (characterized by an ADN chain), but our analysis showed that this was not necessary. Its main drawbacks are that it is random rather than deterministic as is the backpropagation algorithm, and that requires to carefully adjust many parameters (rate of mutations, size of the population). These parameters have been adjusted following two requirements: efficiency of the learning and stability of the result. Learning by genetic algorithms allows therefore to impose the theoretical constraints from the chiral sum rules in a natural way.

Finally, we would like to comment on a new learning algorithm that implements the main advantages of the two methods, that is, it is deterministic and therefore the simulation time is smaller, but at the same time it supports non local contributions to the error function. During the realization of this work, a novel technique was developed that allowed to use the backpropagation learning algorithms in the case of eq. (53), when the error functional contains convolutions. This allowed to check that the results obtained with the genetic algorithm approach were correct. In brief, this technique consists in noticing that an integral can be determined up a any desired accuracy by a finite sum of local contributions, when in this context local means that only depends on one network output. In fact, this is what any numerical integration method does, so it is clear that training algorithms for backpropagation learning can be implemented. The result is that for a discretization of the integral of the form

$$𝑑xf(x)=\underset{j=1}{\overset{n}{}}c_jf_j,$$
(54)
where the coefficients depend on the method, applying the usual backpropagation condition (variation of weights and thresholds in the direction of steepest descent of the error function) to the convolution term one finds that corresponding equations are the backpropagation equations but with eq. (49) replaced by

$$\mathrm{\Delta }_1^{(L)k}=2\left(\underset{j=1}{\overset{n_s}{}}c_jf_bA\right)c_a\frac{df}{do_k}zg^{}(h_{1,k}^{(L)}),$$
(55)
where for simplicity we have only considered one sum rule and $`z`$ is a normalization factor present due to the fact that the inputs and outputs of the neural network are normalized, so that the activation function of the neurons are always within the sensibility range of the activation function. In eq. (55) $`o_k`$ means the output of the network when the input $`x_k`$ is introduced, $`f(o(x))`$ is the convolution that we want the network to learn and A is its theoretical value. In this equation each term should be understood as a new pattern for the neural network training. From here the usual backpropagation techniques apply as usual. This novel technique was implemented in the present analysis but it did not improve neither the quality nor the speed of the training, so the genetic algorithms technique was mantained for the training with convolutions. This technique, that is called backpropagation for convolutions, has many applications, and will be the subject of future work.


format: text

text:
# XTRA: The fast X-ray timing detector on XEUS

## 1 Introduction

The X-rays generated in the inner accretion flows around black holes and neutron stars carry information about regions of strongly curved spacetime. This is a regime in which there are important predictions of general relativity still to be tested, such as the existence of an innermost stable circular orbit. X-ray spectroscopy and fast timing studies can both be used to diagnose the orbital motion of the accreting matter in the immediate vicinity of the collapsed star, where the effects of strong gravity become important.

With the discovery of millisecond aperiodic X-ray time variability from accreting black hole and neutron star X-ray binaries (quasi-periodic oscillations: QPOs), and brightness burst oscillations in neutron stars, RXTE (Bradt et al. (1993)) has clearly demonstrated that fast X-ray timing has the potential to constrain the mass and radius of neutron stars and measure the motion of matter in strong gravity fields (see van der Klis (2000) and e.g. Lamb these proceedings). As pointed out by van der Klis, it is now time to turn these diagnostics into true tests of general relativity. For this, a follow-up to the RXTE Proportional Counter Array (PCA) should have a collecting area of at least ten times larger (more than 6 m2) and an energy resolution at least five times better (i.e. $`200`$ eV at the Iron line).

There are two possible implementations for such an instrument. One is collimated, in which the collecting area is the same as the detector area (as in the PCA), and one in which the two are decoupled through the use of focussing optics. In the first case, beside proportional counters (e.g. Zhang, these proceedings), large area thick Silicon PIN diodes could be considered (e.g. Barret (2001); Kaaret (2001) and Kaaret, these proceedings), although it is unclear whether they will provide an improvement in energy resolution and whether the background could be kept as low as in PCA-like proportional counters. On the other hand, with the use of focussing optics, the immediate advantage is that the detector can be made very small with low background and good energy resolution, allowing timing and spectroscopy to be carried out simultaneously.

This paper describes the XEUS Timing for Relativistic Astrophysics (XTRA) instrument considered for the focal plane instrumentation of the XEUS mission Barret (2003a). Combining the huge collecting area (30 m2) and broad band pass (up to $`30`$ keV) of the XEUS mirror optics with Silicon drift chambers, XTRA will fulfill the two main requirements for a follow-up to the PCA. The science case for such an instrument was presented by van der Klis in the context of XEUS van der Klis (2003), and will be emphasized again in many papers in these proceedings (see the contributions by Kaaret, van der Klis, Lamb, Miller, Psaltis, Strohmayer, …). It will not be repeated here. Instead in this paper, first I will briefly summarize the current XEUS mission profile, then I will present some sensitivity estimates for XTRA and show the results of some simulations of QPOs and burst oscillations, using the expected count rates for XTRA. Finally, I will describe the detector implementation (see also Barret (2003b); Staubert (2003)).

## 2 XEUS in brief

XEUS: The X-ray Evolving Universe Spectroscopy mission represents a potential follow-on mission to the ESA XMM$``$Newton cornerstone observatory currently in orbit. The XEUS mission was considered as part of the ESA Horizon 2000 plus program within the context of the International Space Station (ISS). XEUS is the next logical step forward in X-ray astrophysics after the current set of great observatories (XMM$``$Newton and Chandra), have completed their operational lives. The scientific objectives of XEUS are however so demanding that the mission will clearly represent a major technological challenge compared to past astrophysics missions. Its development and ultimate success rely on the capability to achieve a key breakthrough in the size of an optic capable of entering orbit.

The primary aim of XEUS is the astrophysics of the most distant discrete objects in the Universe (see Hasinger (2000) for details). The specific scientific issues, which XEUS will address can be summarized as follows: 1) To measure the spectra of objects with a redshift z$`>4`$ at flux levels below 10-18 ergs cm-2 s-1, which is just about 1000 times fainter than XMM$``$Newton 2) To determine from the X-ray spectral lines the redshift and thus age of these very faint objects that may not have easily identified optical counterparts 3) To establish the cosmological evolution of matter in the early Universe through the very clear means of the study of heavy element abundances as a function of redshift, i.e. the role of element evolution as the Universe aged through galaxy formation in the associated early stellar processes.

To meet these scientific objectives, XEUS consists of two free flying satellites: a mirror spacecraft (MSC) and a detector spacecraft (DSC) separated by 50 meters and aligned in orbit by an active orbital control and alignment system Bavdaz (2003); Parmar (2003). In the baseline mission scenario, the required large aperture mirrors (see Figure 1) cannot be deployed in a single launch, making XEUS a two step mission, in which the growing of the mirrors involves the facilities in place on the ISS. However, given the uncertainties related to the future of the ISS, alternative mission scenarios are currently being considered at ESA.

In the baseline, XEUS$``$1 comprising the mated mirror and detector satellites will be launched by an Ariane 5 type of launcher in a low earth orbit, with an inclination similar to the ISS. After having completed a few years of observations, MSC1 will go docking on the ISS, where the external mirror segments are waiting. Using the robotic arms on the ISS, the mirror segments will be added to MSC1 to reach the final configuration of MSC2. A second detector spacecraft will then be launched to operate at the focus of MSC2: XEUS$``$2 will then be born. The effective areas of XEUS$``$1 and XEUS$``$2 Aschenbach (2001) are shown in Figure 1, and compared with the effective area of the RXTE/PCA (5 units combined). The proposed high energy extension in which the inner mirror shells are coated with supermirrors is also shown. This yields $`20000`$ cm2 at $`9`$ keV and still $`1700`$ cm2 at 30 keV.

Although XEUS was designed to explore the most distant regions of the Universe, we now show how such an ambitious mission also has great potential for studying the brightest sources in the sky.

## 3 XTRA sensitivity for timing studies

For the sensitivity computations, we have assumed that the timing detector is made of 300 microns of Silicon lying above 2 mm of CdZnTe (see below). This produces a more or less flat detector response up to 80 keV. Table 1 gives the count rates expected from some sources, accounting for typical source spectra and interstellar absorption (for galactic sources, this suppresses most of the photons below $`1`$ keV). The Crab produces about 14000 counts/s in the PCA, which is just about 17 times less than for XEUS$``$1 and 57 times less than for XEUS$``$2 (larger conversion factors apply for spectra softer than the Crab spectrum; 20 is a good average).

### 3.1 QPO detection sensitivity

Let us now compute the sensitivity for QPO detection. The signal to noise ratio $`n_\sigma `$ at which a QPO is detected in a photon counting experiment is approximately:

$$n_\sigma =\frac{1}{2}\frac{S^2}{B+S}r_S^2\left(\frac{T}{\mathrm{\Delta }\nu }\right)^{1/2}$$

where $`S`$ and $`B`$ are the source and background count rates, respectively, $`r_S`$ is the (RMS) root mean squared amplitude of the variability expressed as a fraction of $`S`$, $`T`$ the integration time and $`\mathrm{\Delta }\nu `$ the bandwidth of the variability. The coherence time of the signal is related to the width of the QPO as $`\tau =1/\pi \mathrm{\Delta }\nu `$.

From the above formulae, assuming B $`0`$ appropriate for XTRA, one can estimate the RMS amplitude corresponding to a $`5\sigma `$ QPO detection as a function of the source count rate (Figure 2). One can also do the same computations for the PCA (assuming a background of 100 counts/s) and for an EXTRA-like instrument of the type proposed to ESA in response to the call for F2/F3 mission proposals Barret (2001). EXTRA was designed as a collimated instrument with a detector made of large area Si PIN diodes, covering a total effective area of 6.7 m2 (see also Kaaret, these proceedings for a similar design). Estimated background rate for such an instrument in a low earth orbit gave about 2500 count/s. Figure 2 shows that i) XTRA provides better than one order of magnitude sensitivity improvement in terms of RMS for QPO detections over the PCA ii) thanks to the negligible background, XTRA will enable the timing of very faint X-ray sources (down to less than 0.1 mCrab; e.g. binaries in external galaxies), whereas the background of an EXTRA like mission is of the order of 5 mCrab.

Given the scaling of the above formula, and the increase of count rates between the PCA and XTRA, a QPO detected at $`5\sigma `$ with the PCA will be detected at $`100\sigma `$ with XEUS$``$1. Similarly, XTRA on XEUS$``$1 will detect signals at the same level of significance as the PCA but for an observing time $`400`$ times shorter, and for the strongest signals over tens of milliseconds. XTRA will thus enable us to detect kilo-Hz QPOs at much lower amplitudes, especially in the domain in which the PCA loses the signal towards the highest frequencies van der Klis (2000). This is a region where we might be able to see the saturation of the frequency at the innermost stable circular orbit (which gives a measure of the neutron star mass), if the orbital interpretation for the signal is correct. One might also detect QPOs at frequency above 1330 Hz which is the maximum observed today van Straaten (2000). As stressed by Miller in these proceedings, detection of QPOs at frequencies as high as 1800 Hz would allow us to eliminate all standard nucleonic or hybrid quark matter equations of state, leaving only strange stars.

As stressed above, QPOs will be detected on very short timescales, and depending on the nature of the signal, possibly within their coherence times. This would by itself allow us to discriminate between the various classes of models van der Klis (2003), yielding in some cases constraints on the compact object mass, angular momentum and orbital radius at which the QPOs are produced (see e.g. Abramowicz & Kluzniak, Miller these proceedings). For illustrative purposes, we have first carried out a simulation of a QPO, which is a pulsar like signal of constant amplitude (6.5% RMS for a source count rate of 3000 count/s in the PCA) for which the phase is changed randomly every 200 cycles. The frequency of the signal shifts at a constant rate of 0.25 Hertz per second. This produces a QPO of width $`4`$ Hz corresponding to a quality factor Q$`=\nu /\mathrm{\Delta }\nu 200`$. Q values as large as 200 have already been reported for QPOs measured in power density spectra integrated over 100 seconds or more (e.g. Berger (1996) for 4U1608$``$52). Here we have computed power density spectra over 1 second time intervals, and produced a dynamical power density spectrum for both the PCA and XTRA, assuming a count rate of 45000 counts/s for the latter. The result is shown in Figure 3. It clearly shows that, although the frequency evolution of the QPO can be tracked in the PCA image, only XTRA can measure its frequency and width on short time scales.

We have carried out a second set of illustrative simulations to demonstrate that under somewhat favorable conditions, one can localize the signal directly in the time domain. This time the QPO is made of a succession of wavetrains of finite lifetime (QPO at 830 Hz keeping the same frequency for 200 cycles, the RMS of the QPO is again fixed at 6.5% for a source at 3000 count/s in the PCA and there is one wavetrain per second). To explore the signal in the time domain, we slide a window of 0.4 second width over the data with a time step of 0.1 second. Then we search for a periodic signal using a standard Rayleigh test. The periodograms (2.5 Hz of resolution) so obtained are stacked together to form a time-frequency image. One expects that when the window is centered on the QPO wavetrains, the periodogram will reach a maximum at the QPO frequency, and this will appear as a strong excess in the image. A comparison between a time-frequency image obtained for PCA and XTRA count rates is shown in Figure 4. Whereas it is impossible to locate the signal in the PCA image, it stands out very clearly in the XTRA image. Once located in the time domain, one can look at the time profile of the oscillations (in the present simulation folding the data produces a nice sinusoid) and study what causes the signal to lose its coherence. One can also do waveform fitting, correlated timing and spectral analysis and so on. With sophisticated time-frequency analysis tools (singular spectrum analysis, maximum enthropy methods, …), even under less favorable conditions (i.e. shorter coherence times, rapid frequency drifts, superposition of wavetrains, …), with XTRA it will be possible to investigate the true nature of the signal by looking for the first time at its properties in the time domain.

In general, when plotting the QPO RMS as a function of energy, it is found that it increases, to reach a saturation around 10$``$20 keV (e.g. Figure 5 from Berger (1996)). One limitation of the use of focussing optics is their energy band pass which drops off quickly after 10 keV (see Figure 1). In the case of XEUS however, the use of advanced super mirrors extends the energy response up to 30 keV. In fact, the effective area of the XEUS mirrors is comparable or even larger than the PCA (5 units) up to 30 keV. In the current detector implementation for XTRA, the use of CdZnTe below the Silicon would match the high energy response of the mirrors. One can compute the limiting RMS (using equation above), as a function of energy for a given source energy spectrum. This is shown in Figure 5 and compared with the observed data points for the neutron star X-ray binary 4U1608$``$52 Berger (1996). This figure shows that thanks to the high energy extension of the mirrors, XTRA will be more sensitive than the PCA over the whole energy range QPOs were detected.

### 3.2 Coherent signal detection sensitivity

For a coherent signal ($`T`$$`<`$$`1/\mathrm{\Delta }\nu `$), the familiar exponential detection regime applies, with false-alarm probability $``$$`\mathrm{exp}[S^2r_S^2T/2(B+S)]`$. One can compute the RMS for the detection of a coherent signal at a given false alarm probability (Figure 6). The figure shows that XTRA would be able to detect pulsations at the 0.01% RMS level in Sco X$``$1. Why so far no classical low-mass X-ray binaries, such as Sco X$``$1, have shown pulsations in their persistent emission remains very puzzling. Several competing explanations have been put forward, as for example that they may be rotating at sub-millisecond periods or may contain compact neutron stars van der Klis (2003). Thanks to its improved sensitivity for periodic signals XTRA will enable us to test these hypothesis, which have their own constraints on the equation of state of dense matter.

In addition, it has been recently suggested that the neutron star spin frequencies inferred from burst oscillations have an upper limit of $``$750 Hz Chakrabarty et al. (2003). If this is because of angular momentum removed by gravitational radiation, then searches for periodic gravitational waves emitted by these systems require highly accurate measurements of their spin periods. XTRA could naturally achieve this task.

Concerning burst oscillations, XTRA with its optimum band pass for X-ray bursts, will allow the oscillations to be detected within one cycle (see Figure 7). These oscillations are probably caused by rotational modulation of a hot spot on the stellar surface. The emission from the hot spot is affected by Doppler boosting, relativistic aberration and gravitational light bending (whose magnitude increases with the neutron star compactness, see Figure 7). By fitting the waveform, it will be possible to investigate the spacetime around the neutron star, and simultaneously constrain its mass and radius, and hence determine the equation of state of its high density core (see Strohmayer, these proceedings, but see also Poutanen for an example of what can be done by studying the energy dependence of pulse profiles of accreting millisecond pulsars).

## 4 Detector implementation

To meet the scientific objectives of XTRA, its detector must be capable of handling up to 3 Mcts/s (XEUS$``$1) and 10 Mcts/s (XEUS$``$2) (equivalent to a 10 Crab source, see Table 1) with a timing resolution of $`10\mu `$s and a deadtime as low as possible, with the requirement that the latter must be measured accurately. In addition, the detector energy range should match closely the high energy response of the mirrors.

### 4.1 Silicon Drift Detector

Among the fast X-ray detectors currently available, Silicon Drift Detectors (SDDs) are the most promising Strüder (2000); Lechner (2001). The SDD is a completely depleted volume of Silicon in which an arrangement of increasingly negative biased rings drive the electrons generated by the impact of ionising radiation towards a small readout node in the center of the device. The time needed for the electrons to drift is much less than 1 $`\mu `$s. The main advantage of SDDs over conventional PIN diodes is the small physical size and consequently the small capacitance of the anode, which translates to a capability to handle high count rates simultaneously with good energy resolution. To take full advantage of the small capacitance, the first transistor of the amplifying electronics is integrated on the detector chip (see Figure 8). The stray capacitance of the interconnection between the detector and amplifier is thus minimized, and furthermore the system becomes practically insensitive to mechanical vibrations and electronic pickup.

A spectrum of a 5 mm2 SDD, with a customized amplifying electronics, is shown in Figure 9 to illustrate that energy resolution as good as 130 eV can be achieved, even at room temperature. Energy resolution of better than $`200`$ eV can be maintained for count rates up to 105 cts/s (e.g. Lechner (2001)).

With such an energy resolution (comparable to what CCDs are currently providing), spectroscopy and timing can be used as independent tools to probe the strong gravity regions around compact objects. A particular emphasis should be put on the relativistically broadened Fe K$`\alpha `$ lines, recently detected by XMM$``$Newton and Chandra (e.g. Miller (2003)). The line profile which is affected by several physical processes (Doppler shifts, beaming, gravitational redshifts) gives complementary information about the same region, as explored with millisecond timing. The lines are relatively broad (1 keV or more) so that the 200 eV energy resolution of XTRA is sufficient. This again offers independent opportunities to constrain the compact object mass and spin, as well as the radius at which the line is produced. Similarly, as pointed out by Strohmayer (these proceedings) lines in X-ray bursts Cottam et al. (2002) might be detectable, even with a 200 eV spectral resolution. Measurement of a gravitational redshift yields a direct estimate of the compactness of the neutron star, and yet another constraint on the equation of state of dense matter (see Strohmayer, these proceedings).

Finally with a rather conservative low energy threshold $`1`$ keV, XTRA will be very sensitive to very soft sources (e.g. quiescent neutron star transients, supersoft sources, millisecond puslars, isolated neutron stars …), in an energy domain never explored before (2.5 keV is the current energy threshold of the RXTE/PCA).

### 4.2 The XTRA detector

For timing studies, deadtime is always a critical issue. Deadtime will include contributions from the signal rise time, the charge sensitive amplifier, and the shaping amplifier. The first two of these can be very short, and the limiting contribution is that of the amplifier, where a trade-off between speed and energy resolution is necessary. Shaping time constants as short as 50 nanoseconds (ns) have been found to be usable Strüder (2000). The challenging requirement for XTRA is to achieve an accurately measured deadtime of $`500`$ ns per event. Using currently available devices and pipelining techniques, the analog-digital conversion stage is not a limiting factor at these speeds.

A 500 ns deadtime per event corresponds to a 5% deadtime for a source producing $`10^5`$ cps/s. To handle $`10^6`$ cps/s with a reasonable deadtime, one must therefore distribute the focal beam over $`10`$ detectors or more. The best and easiest solution could thus be a detector made of an ensemble of about $`10`$ separate SDDs on a single wafer. Such SDD arrays already exist, as shown in Figure 10. This detector should therefore be operated out of focus. For XEUS$``$1, the out of focus distance is of the order of 10 cm. This could be accomplished either by a mechanical construction, or by changing the distance between the mirror and detector satellites. Although this will require a careful study, both solutions appear to be feasible within the current XEUS mission design. The requirements in terms of real estate on the detector spacecraft are not constraining, in particular because no complicated cooling systems will be necessary.

As mentioned above, a high energy extension (above 10 keV) is proposed for the mirrors Aschenbach (2001). SDDs are currently produced with thicknesses of 300 microns, which is adequate to cover the energy range below 10 keV. Although there are on-going efforts to make thicker devices, the best match of the high energy response of the mirrors will require the SDD array to be associated with a higher density detector located underneath. Among the potential high energy semi-conductor detectors, CdZnTe Budtz-Jørgensen (2001) stands today as a very promising solution. Such a detector would both ensure the overlap in energy range with the SDD array, and provide a flat energy response up to $`80`$ keV and 10 microsecond timing resolution Budtz-Jørgensen (2001). Count rates should not exceed a few thousands counts/s in the CdZnTe detector.

The goal for XTRA is to send to the ground the time and energy information of every photon. For most sources, data compression will make this possible (within a 2 Mbits/s data rate) without compromising either time or energy resolution. For the very brightest sources, this can still be done with a restricted number of energy channels. However, within the timeframe of XEUS, on-board memory and telemetry rate should not be considered as real issues.

The detector will be exposed to high radiation doses during operations and one must therefore consider its radiation hardness. The main limitation in the maximum acceptable dose arises from the JFET connected to the collecting anode on the back of the device Leutenegger (2000). High energy photons absorbed in the transistor region increase the amount of oxyde charge and interface traps, thus reducing the charge carrier lifetimes, and thus contributing to increase the leakage current. Laboratory measurements indicate however that a 300 micron thick SDD survives a radiation dose of $`10^{13}`$ incoming high energy photons (E$`>12`$ keV) Leutenegger (2000). This is equivalent to a continuous exposition of 3 years at 105 photons/s. Concerning particles, the XMM$``$Newton EPIC pn cameras, which have similar detector technology is performing extremely well in space Strüder (2001). So, the device selected can be clearly considered as radiation hard.

To get ready when XEUS gets approved in the ESA science program, we have proposed to the French Space Agency an R&T program, aimed at characterizing the performance of the detector for XTRA and starting to work on deadtime issues and on-board data processing. The proposal was approved in 2002 (before CNES experienced internal problems). Hopefully funding should arrive in 2004. This hardware program will be carried out in close collaboration with MPE/MPI, University of Tübingen and the ROENTEC company which has a great expertise on high count rate electronics for Silicon drift detectors.

## 5 Conclusions

The XEUS Timing for Relativistic Astrophysics instrument, by combining the huge aperture of the XEUS mirrors with a silicon drift detector array in the focal plane meets the two most important science requirements for a follow-up to the RXTE/PCA. It will have a collecting area more than ten times larger than the one of the PCA, and a much improved energy resolution. Through fast X-ray timing and spectroscopy, XTRA will thus provide the first opportunity to test general relativity in the strong gravity field regime and constrain with unprecedented accuracy the equation of state of dense matter.

I am extremely grateful to the following colleagues who gave their support to the proposal for a fast X-ray timing detector in XEUS: J.L. Atteia, T. Belloni, H. Bradt, L. Burderi, S. Campana, A. Castro-Tirado, D. Chakrabarty, P. Charles, S. Collin, S. Corbel, C. Done, G. Dubus, M. Gierlinski, J. Grindlay, A. Fabian, R. Fender, E. Gourgoulhon, J.M. Hameury, C. Hellier, E. Kendziorra, W. Kluzniak, E. Kuulkers, S. Larsson, J.P. Lasota, T. Maccarone, D. de Martino, K. Menou, C. Miller, F. Mirabel, M. Nowak, J.F. Olive, S. Paltani, R. Remillard, J. Rodriguez, R. Rothschild, T. di Salvo, R. Sunyaev, M. Tagger, M. Tavani, L. Titarchuk, G. Vedrenne, N. White, R. Wijnands, J. Wilms, A. Zdziarski, W. Zhang. I am also grateful to G. Hasinger, A. Parmar and M. Bavdaz for helpful information about the XEUS mission in general. I thank J.F. Olive and J.P. Chabbert for helping me out to produce some of the simulations presented in this paper, and M. Ehanno and O. Godet for providing me with the spectrum of the SDD. Many thanks also to J. Wilms and C. Miller for comments on the manuscript.