text:
# Generalization properties of finite size polynomial Support Vector Machines

## I Introduction

In the last decade, the typical properties of neural networks that learn classification tasks from a set of examples have been analyzed using the approach of Statistical Mechanics. In the general setting, the value of a binary output neuron represents whether the input vector, describing a particular pattern, belongs or not to the class to be recognized. Manuscript character recognition and medical diagnosis are examples of such classification problems. The process of inferring the rule underlying the input-output mapping given a set of examples is called learning. The aim is to predict correctly the class of novel data, i.e. to generalize.

In the simplest neural network, the perceptron, the inputs are directly connected to a single output neuron. The output state is given by the sign of the weighted sum of the inputs. Then, learning amounts to determine the weights of the connexions in order to obtain the correct outputs to the training examples. Considering the weights as the components of a vector, the network classifies the input vectors according to whether their projections onto the weight vector are positive or negative. Thus, patterns of different classes are separated by the hyperplane orthogonal to the weight vector. Beyond these linear separations, two different learning schemes have been suggested. Either the input vectors are mapped by linear hidden units to so called internal representations that must be linearly separable by the output neuron, or a more powerful output unit is defined, able to perform more complicated functions than just the weighted sum of its inputs.

The first solution is implemented using feedforward layered neural networks. The classification of the internal representations, performed by the output neuron, corresponds in general to a complicated separation surface in input space. However, the relation between the number of hidden units of a network and the class of rules it can infer is still an open problem. In practice, the number of hidden neurons is either guessed or determined through constructive heuristics.

A solution that uses a more complex output unit, the Support Vector Machine (SVM) , has been recently proposed. The input patterns are transformed into high dimensional feature vectors whose components may include the original input together with specific functions of its coordinates selected a priori, with the aim that the learning set be linearly separable in feature space. In that case the learning problem is reduced to that of training a simple perceptron. For example, if the feature space includes all the pairwise products of the input vector, the SVM may implement any classification rule corresponding to a quadratic separating surface in input space. Higher order polynomial SVMs and other types of SVMs may be defined by introducing the corresponding features. A big advantage is that learning a linearly separable rule is a convex optimization problem. The difficulties of having many local minima, that hinder the process of training multilayered neural networks, are thus circumvented. Once the adequate feature space is defined, the SVM selects the particular hyperplane called Maximal Margin (or Maximal Stability) Hyperplane (MMH), which lies at the largest distance to its closest patterns in the training set. These patterns are called Support Vectors (SV). The MMH solution has interesting properties . In particular, the fraction of learning patterns that belong to the SVs provides an upper bound  to the generalization error, that is, to the probability of incorrectly classifying a new input. It has been shown  that the perceptron weights are a linear combination of the SVs, an interesting property in high dimensional feature spaces, as their number is bounded.

A perceptron can learn with very high probability any set of examples, regardless of the underlying classification rule, provided that their number does not exceed twice its input space dimension . However, this simple rote learning does not capture the rule underlying the classification. As it may arise that the feature space dimension of the SVM is comparable to, or even larger than, the number of available training patterns, we would expect that SVMs have a poor generalization performance. Surprisingly, this seems not to be the case in the applications .

Two theoretical papers  have recently addressed this interesting question. They determined the typical properties of a family of polynomial SVMs in the limit of large dimensional spaces, reaching completely different results in spite of the seemingly innocuous differences between the models. Both papers consider polynomial SVMs in which the input vectors $`𝐱\text{IR}^N`$ are mapped onto quadratic features $`𝚽`$. More precisely, the normalized mapping $`𝚽_n(𝐱)=(𝐱,x_1𝐱/\sqrt{N},x_2𝐱/\sqrt{N},\mathrm{},x_N𝐱/\sqrt{N})`$ has been considered in . The non-normalized mapping $`𝚽_{nn}(𝐱)=(𝐱,x_1𝐱,x_2𝐱,\mathrm{},x_k𝐱)`$ has been studied in  as a function of $`k`$, the number of quadratic features. For $`k=N`$ the dimension of both feature spaces is the same, corresponding to a linear subspace of dimension $`N`$, and a quadratic subspace of dimension $`N^2`$. The mappings only differ in the distributions of the quadratic components in feature space. Due to the normalization, those of $`𝚽_n`$ are squeezed by a normalizing factor $`a=1/\sqrt{N}`$ with respect to those of $`𝚽_{nn}`$. In the case of learning a linearly separable rule with the non-normalized mapping $`𝚽_{nn}`$, the generalization error at any given learning set size increases dramatically with the number $`k`$ of quadratic features included . On the contrary, in the case of mapping $`𝚽_n`$, the generalization error exhibits an interesting stepwise decrease, also found within the Gibbs learning paradigm in a quadratic feature space . If the number of training patterns scales with $`N`$, the dimension of the linear subspace, it decreases up to an asymptotic lower bound. If the number of examples scales proportionally to $`N^2`$, it vanishes asymptotically. In particular, if the rule to be inferred is linearly separable in the input space, learning in the feature space with the mapping $`𝚽_n`$ is harmless, as the decrease of the generalization error with the number of training patterns presents a slight slow-down with respect to that of a simple perceptron learning in input space.

As this stepwise learning is exclusively related to the fact that the normalizing factor of the quadratic features vanishes in the thermodynamic limit $`N\mathrm{}`$, in the present paper we determine the influence of the normalizing factor on the typical generalization performance of finite size SVMs. To this end, we introduce two parameters, $`\sigma `$ and $`\mathrm{\Delta }`$, caracterizing the mapping of the $`N`$-dimensional input patterns onto the feature space. The variance $`\sigma `$ reflects the width of the high-order features distribution and is related to the normalizing factor $`a`$. The inflation factor $`\mathrm{\Delta }`$ accounts for the proportion of quadratic features with respect to the input space dimension $`N`$. Actual quadratic SVMs are caracterized by different values of $`\mathrm{\Delta }`$ and $`\sigma `$, depending on $`N`$ and $`a`$. Keeping $`\sigma `$ and $`\mathrm{\Delta }`$ fixed in the thermodynamic limit allows us to determine the typical properties of actual SVMs, which have finite compressing factors and inflation ratios.

In fact, the behaviour of the SVMs is the same as that of a simple perceptron learning a training set with patterns drawn from a highly anisotropic probability distribution, such that a macroscopic fraction of components have a different variance from the others. Not surprisingly, we find that the asymptotic behaviour corresponding to both the small and large training set size limits, is the same as the one of the perceptron’s MMH. Only the prefactors depend on the mapping used by the SVM.

As expected, the stepwise learning obtained with the normalized mapping in the thermodynamic limit becomes a crossover. Upon increasing the number of training patterns, the generalization error first present an abrupt decrease, that corresponds to learning the weight components in the linear subspace, followed by a slower decrease corresponding to the learning of the quadratic components. The steepness of the crossover not only depends on $`\mathrm{\Delta }`$ and $`\sigma `$, but also on the task to be learned. The agreement between our analytic results and numerical simulations is excellent.

The paper is organized as follows: in section II we introduce the model and the main steps of the Statistical Mechanics calculation. Numerical simulation results are compared to the corresponding theoretical predictions in section III. The two regimes of the generalization error and the asymptotic behaviours are discussed in section IV. The conclusion is left to section V.

## II The model

We consider the problem of learning a binary classification task from examples with a SVM in polynomial feature spaces. The learning set contains $`M`$ patterns $`(𝐱^\mu ,\tau ^\mu )`$ ($`\mu =1,\mathrm{},M`$) where $`𝐱^\mu `$ is an input vector in the $`N`$-dimensional input space, and $`\tau ^\mu \{1,1\}`$ is its class. We assume that the components $`x_i^\mu `$ ($`i=1,\mathrm{},N`$) are independent identically distributed (i.i.d.) random variables drawn from gaussian distributions having zero-mean and unit variance:

$$P(𝐱)=\underset{i=1}{\overset{N}{}}\frac{1}{\sqrt{2\pi }}\mathrm{exp}\left(\frac{x_i^2}{2}\right).$$
(1)
In the following we concentrate on quadratic feature spaces, although our conclusions are more general, and may be applied to higher order polynomial SVMs, as discussed in section IV. The mappings $`𝚽_{nn}(𝐱)=(𝐱,x_1𝐱,x_2𝐱,\mathrm{},x_N𝐱)`$ and $`𝚽_n(𝐱)=(𝐱,x_1𝐱/\sqrt{N},x_2𝐱/\sqrt{N},\mathrm{},x_N𝐱/\sqrt{N})`$ are particular instances of mappings of the form $`𝚽(𝐱)=(\varphi _1,\varphi _2,\mathrm{},\varphi _N,\varphi _{11},\varphi _{12},\mathrm{},\varphi _{NN})`$ where $`\varphi _i=x_i`$, and $`\varphi _{ij}=ax_ix_j`$, where $`a`$ is the normalizing factor of the quadratic components: $`a=1`$ for mapping $`𝚽_{nn}`$ and $`a=1/\sqrt{N}`$ for $`𝚽_n`$. The patterns probability distribution in feature-space is:

$$P\left(𝚽\right)=\underset{i=1}{\overset{N}{}}\frac{dx_i}{\sqrt{2\pi }}\mathrm{exp}\left(\frac{x_i^2}{2}\right)\delta (\varphi _ix_i)\underset{j=1}{\overset{N}{}}\delta \left(\varphi _{ij}ax_ix_j\right).$$
(2)
Clearly, the components of $`𝚽`$ are not independent random variables. For example, a number $`O(N^3)`$ of triplets of the form $`\varphi _{ij}\varphi _{jk}\varphi _{ki}`$ have positive correlations. These contribute to the third order moments, which should vanish if the features were gaussian. Moreover, the fourth order connected correlations  do not vanish in the thermodynamic limit. Nevertheless, in the following we will neglect these and higher order connected moments. This approximation, used in  and implicit in , is equivalent to assuming that all the components in feature space are independent gaussian variables. Then, the only difference between the mappings $`𝚽_n`$ and $`𝚽_{nn}`$ lies in the variance of the quadratic components distribution. The results obtained using this simplification are in excellent agreement with the numerical tests described in the next section.

Since, due to the symmetry of the transformation, only $`N(N+1)/2`$ among the $`N^2`$ quadratic features are different, hereafter we restrict the feature space and only consider the non redundant components, that we denote $`𝝃=(𝝃_u,𝝃_\sigma )`$. Its first $`N`$ components $`𝝃_u=(\xi _1,\mathrm{},\xi _N)`$ hereafter called u-components, represent the input pattern of unit variance, lying in the linear subspace. The remaining components $`𝝃_\sigma =(\xi _{N+1},\mathrm{},\xi _{\stackrel{~}{N}})`$ stand for the non redundant quadratic features, of variance $`\sigma `$, hereafter called $`\sigma `$-components. $`\stackrel{~}{N}`$ is the dimension $`\stackrel{~}{N}=N(1+\mathrm{\Delta })`$ of the restricted feature space, where the inflation ratio $`\mathrm{\Delta }`$ is the relative number of non-redundant quadratic features per input space dimension. The quadratic mapping has $`\mathrm{\Delta }=(N+1)/2`$.

According to the preceding discussion, we assume that learning $`N`$-dimensional patterns selected with the isotropic distribution (1) with a quadratic SVM is equivalent to learning the MMH with a simple perceptron in an $`\stackrel{~}{N}`$-dimensional space where the patterns are drawn using the following anisotropic distribution,

$$P\left(𝝃\right)=\underset{i=1}{\overset{N}{}}\frac{1}{\sqrt{2\pi }}\mathrm{exp}\left(\frac{\xi _i^2}{2}\right)\underset{j=N+1}{\overset{\stackrel{~}{N}}{}}\frac{1}{\sqrt{2\pi \sigma ^2}}\mathrm{exp}\left(\frac{\xi _j^2}{2\sigma ^2}\right).$$
(3)
The second moment of the u-features is $`𝝃_u^2=N`$ and that of the $`\sigma `$-features is $`𝝃_\sigma ^2=N\mathrm{\Delta }\sigma ^2`$. If $`\sigma ^2\mathrm{\Delta }=1`$, we get $`𝝃_\sigma ^2=𝝃_u^2`$, which is the relation satisfied by the normalized mapping considered in . The non-normalized mapping corresponds to $`\sigma ^2\mathrm{\Delta }=N`$. In the following, instead of selecting either of these possibilities a priori, we consider $`\mathrm{\Delta }`$ and $`\sigma `$ as independent parameters, that are kept constant when taking the thermodynamic limit.

Since the rules to be inferred are assumed to be linear separations in feature space, we represent them by the weights $`𝐰^{}=(w_1^{},w_2^{},\mathrm{},w_{\stackrel{~}{N}}^{})`$ of a teacher perceptron, so that the class of the patterns is $`\tau =\mathrm{sign}(𝝃𝐰^{})`$. Without any loss of generality we consider normalized teachers: $`𝐰^{}𝐰^{}=\stackrel{~}{N}`$. The training set in feature space is then $`_M=\{(𝝃^\mu ,\tau ^\mu )\}_{\mu =1,\mathrm{},M}`$.

In the following we study the typical properties of polynomial SVMs learning realizable classification tasks, using the tools of Statistical Mechanics. If $`𝐰=(w_1,\mathrm{},w_{\stackrel{~}{N}})`$ is the student perceptron weight vector, $`\gamma ^\mu =\tau ^\mu 𝝃^\mu 𝐰/\sqrt{𝐰𝐰}`$ is the stability of pattern $`\mu `$ in feature space. The pertinent cost function is :

$$E(𝐰,\kappa ;_M)=\underset{\mu =1}{\overset{M}{}}\mathrm{\Theta }(\kappa \gamma ^\mu ).$$
(4)
$`\kappa `$, the smallest allowed distance between the hyperplane and the training patterns, is called the margin. The MMH corresponds to the weights with vanishing cost (4) that maximize $`\kappa `$.

The typical properties of cost (4) in the case of isotropic pattern distributions have been exhaustively studied . The case of a single anisotropy axis has also been investigated . Here we study the case of the anisotropic distribution (3), where a macroscopic fraction of components have different variance from the others, which is pertinent for understanding the properties of the SVM.

Considering the cost (4) as an energy, the partition function at temperature $`1/\beta `$ writes

$$Z(\kappa ,\beta ;_M)=\mathrm{exp}[\beta E(𝐰,\kappa ;_M)]p(𝐰)𝑑𝐰.$$
(5)
Without any loss of generality, we assume that the a priori distribution of the student weights is uniform over the hypersphere of radius $`\stackrel{~}{N}^{1/2}`$, i.e. $`p(𝐰)=\delta (𝐰𝐰\stackrel{~}{N})`$, meaning that the student weights are normalized in feature space. In the limit $`\beta \mathrm{}`$, the corresponding free energy $`f(\kappa ,\beta ;_M)=(1/\beta N)\mathrm{ln}Z(\kappa ,\beta ;_M)`$ is dominated by the weights that minimize the cost (4).

The typical properties of the MMH are obtained by looking for the largest value of $`\kappa `$ for which the quenched average of the free energy over the patterns distribution, in the zero temperature limit $`\beta \mathrm{}`$, vanishes. This average is calculated by the replica method, using the identity

$$f(\kappa ,\beta )=\frac{1}{N\beta }\overline{\mathrm{ln}Z(\kappa ,\beta ;_M)}=\frac{1}{N\beta }\underset{n0}{lim}\frac{\mathrm{ln}\overline{Z^n(\kappa ,\beta ;_M)}}{n},$$
(6)
where the overline represents the average over $`_M`$, composed of patterns selected according to (3).

We obtain the typical properties of the MMH corresponding to given values of $`\mathrm{\Delta }`$ and $`\sigma `$ by taking the thermodynamic limit $`N\mathrm{}`$, $`M\mathrm{}`$, with $`\alpha M/N`$, $`\mathrm{\Delta }`$ and $`\sigma `$ constant. Notice that the relation between the number of training examples and the feature space dimension, $`\stackrel{~}{\alpha }M/\stackrel{~}{N}=\alpha /(1+\mathrm{\Delta })`$, is finite. Thus, not only are we able to study the dependence of the learning properties as a function of the training set size as usual, but also of the inflation factor that characterizes the SVM, as well as of the variance of the quadratic components. As we only consider realizable rules, i.e. classification tasks that are linearly separable in feature space, the energy (4) is a convex function of the weights $`𝐰`$, and replica symmetry holds.

For any $`\kappa <\kappa _{max}`$, there are a macroscopic number of weights that minimize the cost function (4). In particular, in the case of $`\kappa =0`$, the cost is the number of training errors, and is minimized by any weight vector that classifies correctly the training set. The typical properties of such solution, called Gibbs learning, may be expressed in terms of several order parameters . Among them, $`q_u^{ab}=_{i=1}^N\overline{w_i^aw_i^b}/\stackrel{~}{N}`$, $`q_\sigma ^{ab}=_{i=N+1}^{\stackrel{~}{N}}\overline{w_i^aw_i^b}/\stackrel{~}{N}`$ and $`Q^a=_{i=N+1}^{\stackrel{~}{N}}\overline{w_i^aw_i^a}/\stackrel{~}{N}`$, where $`ab`$ are replica indices and $`\mathrm{}`$ stands for the usual thermodynamic average (with Boltzmann factor corresponding to the partition function (5)). $`q_u^{ab}`$ and $`q_\sigma ^{ab}`$ represent the overlaps between different solutions in the u\- and the $`\sigma `$\- subspaces respectively. $`\stackrel{~}{N}Q^a`$ is the typical norm of the $`\sigma `$-components of replica $`a`$. Because of replica symmetry we have $`Q^a`$=$`Q^b=Q`$, $`q_\sigma ^{ab}=q_\sigma `$ and $`q_u^{ab}=q_u`$ for all $`a`$, $`b`$. Upon increasing $`\kappa `$, the volume of the error-free solutions in weight space shrinks, and vanishes when $`\kappa `$ is maximized. Correspondingly, $`q_u1Q`$ and $`q_\sigma Q`$, with $`xlim_{\kappa \kappa _{max}}(1q_u/(1Q))/(1q_\sigma /Q)`$ finite. In the limit of $`\kappa \kappa _{max}`$, the properties of the MMH may be expressed in terms of $`x`$, $`\kappa _{max}`$ and the following order parameters,

$`Q`$ $`=`$ $`{\displaystyle \frac{1}{\stackrel{~}{N}}}{\displaystyle \underset{i=N+1}{\overset{\stackrel{~}{N}}{}}}\overline{w_i^2},`$ (7)
$`R_u`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{(1Q)(1Q^{})}}}{\displaystyle \frac{1}{\stackrel{~}{N}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\overline{w_iw_i^{}},`$ (8)
$`R_\sigma `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{QQ^{}}}}{\displaystyle \frac{1}{\stackrel{~}{N}}}{\displaystyle \underset{i=N+1}{\overset{\stackrel{~}{N}}{}}}\overline{w_iw_i^{}},`$ (9)
where $`Q^{}=_{i=N+1}^{\stackrel{~}{N}}(w_i^{})^2/\stackrel{~}{N}`$ is the teacher’s squared weight vector in the $`\sigma `$-subspace. $`Q`$ is the corresponding typical value for the student. $`R_u`$ and $`R_\sigma `$ are proportional to the overlaps between the student and the teacher weights in the u\- and the $`\sigma `$\- subspaces respectively. The factors in the denominators arise because the weights are not normalized in each subspace.

The saddle point equations corresponding to the extremum of the free energy for the MMH are

$`2{\displaystyle \frac{\alpha }{\mathrm{\Delta }}}\mathrm{\Delta }_\sigma I_1`$ $`=`$ $`(1R_\sigma ^2){\displaystyle \frac{(x+\mathrm{\Delta }_\sigma )^2}{1+\mathrm{\Delta }_\sigma }},`$ (10)
$`2{\displaystyle \frac{\alpha }{\mathrm{\Delta }}}I_2`$ $`=`$ $`\sqrt{{\displaystyle \frac{1+\mathrm{\Delta }_\sigma ^{}}{\mathrm{\Delta }_\sigma ^{}}}}R_\sigma {\displaystyle \frac{x+\mathrm{\Delta }_\sigma }{\sqrt{\mathrm{\Delta }_\sigma (1+\mathrm{\Delta }_\sigma )}}},`$ (11)
$`2{\displaystyle \frac{\alpha }{\mathrm{\Delta }}}Q(1\sigma ^2)I_3`$ $`=`$ $`\left(1x{\displaystyle \frac{1R_\sigma ^2}{1R_u^2}}\right){\displaystyle \frac{x+\mathrm{\Delta }_\sigma }{1+\mathrm{\Delta }_\sigma }},`$ (12)
$`{\displaystyle \frac{R_\sigma ^2}{1R_\sigma ^2}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }_\sigma ^{}}{\mathrm{\Delta }}}{\displaystyle \frac{R_u^2}{1R_u^2}},`$ (13)
$`{\displaystyle \frac{R_u}{R_\sigma }}`$ $`=`$ $`x{\displaystyle \frac{\mathrm{\Delta }}{\sqrt{\mathrm{\Delta }_\sigma \mathrm{\Delta }_\sigma ^{}}}}.`$ (14)
where $`\mathrm{\Delta }_\sigma \sigma ^2Q/(1Q)`$ and $`\mathrm{\Delta }_\sigma ^{}\sigma ^2Q^{}/(1Q^{})`$. The integrals in the left hand side of equations (10-12) are

$`I_1`$ $`=`$ $`{\displaystyle _{\stackrel{~}{\kappa }}^{\mathrm{}}}Dt(t+\stackrel{~}{\kappa })^2H\left({\displaystyle \frac{tR}{\sqrt{1R^2}}}\right),`$ (15)
$`I_2`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}\left[{\displaystyle \frac{\sqrt{1R^2}\mathrm{exp}(\stackrel{~}{\kappa }^2/(2(1R^2)))}{\sqrt{2\pi }}}+\stackrel{~}{\kappa }H\left({\displaystyle \frac{\stackrel{~}{\kappa }}{\sqrt{1R^2}}}\right)\right],`$ (16)
$`I_3`$ $`=`$ $`{\displaystyle _{\stackrel{~}{\kappa }}^{\mathrm{}}}Dt\stackrel{~}{\kappa }(t+\stackrel{~}{\kappa })H\left({\displaystyle \frac{tR}{\sqrt{1R^2}}}\right),`$ (17)
with $`Dtdt\mathrm{exp}(t^2/2)/\sqrt{2\pi }`$, $`H(x)=_x^{\mathrm{}}Dt`$, and

$`\stackrel{~}{\kappa }`$ $`=`$ $`{\displaystyle \frac{\kappa _{max}}{\sqrt{(1Q)(1+\mathrm{\Delta }_\sigma )}}},`$ (18)
$`R`$ $`=`$ $`{\displaystyle \frac{R_u+\sqrt{\mathrm{\Delta }_\sigma \mathrm{\Delta }_\sigma ^{}}R_\sigma }{\sqrt{(1+\mathrm{\Delta }_\sigma )(1+\mathrm{\Delta }_\sigma ^{})}}}.`$ (19)
The value of $`R`$ determines the generalization error through $`ϵ_g=(1/\pi )\mathrm{arccos}(R)`$.

After solving the above equations for $`Q`$, $`R_u`$, $`R_\sigma `$, $`x`$ and $`\stackrel{~}{\kappa }`$, it is straightforward to determine $`\rho _{SV}`$, the fraction of training patterns that belong to the subset of SV :

$$\rho _{SV}=2_{\mathrm{}}^{\stackrel{~}{\kappa }}H\left(tR/\sqrt{1R^2}\right)Dt.$$
(20)
In summary of this section, instead of considering a particular scaling of the fraction of high order features components and their normalization with $`N`$, we analyzed the more general case where these quantities are kept as free parameters. We determined the saddle point equations that define the typical properties of the corresponding SVM. This approach allows us to consider several learning scenarios, and more interestingly, to study the crossover between the different generalization regimes.

## III Results

We describe first the experimental data, obtained with quadratic SVMs, using both mappings, $`𝚽_{nn}`$ and $`𝚽_n`$, which have normalizing factors $`a=1`$ and $`a=1/\sqrt{N}`$ respectively, where $`N`$ is the input space dimension. The $`M=\alpha N`$ random input examples of each training set were selected with probability (1) and labelled by teachers of normalized weights $`𝒘^{}(𝒘_l^{},𝒘_q^{})`$ drawn at random. $`𝒘_l^{}`$ are the $`N`$ components in the linear subspace and $`𝒘_q^{}`$ are the $`N^2`$ components in the quadratic subspace. Notice that, because of the symmetry of the mappings, teachers having the same value of the symmetrized weights in the quadratic subspace, $`(w_{q,ij}^{}+w_{q,ji}^{})/2`$, are all equivalent. The teachers are characterized by the proportion of (squared) weight components in the quadratic subspace, $`Q^{}=𝒘_q^{}𝒘_q^{}/𝒘^{}𝒘^{}`$. In particular, $`Q^{}=0`$ and $`Q^{}=1`$ correspond to a purely linear and a purely quadratic teacher respectively.

The experimental student weights $`𝒘(𝒘_l,𝒘_q)`$ were obtained by solving numerically the dual problem , using the Quadratic Optimizer for Pattern Recognition program , that we adapted to the case without threshold treated in this paper. We determined $`Q`$, and the overlaps $`R_l`$ and $`R_q`$ in the linear and the quadratic subspaces, respectively. For each value of $`M`$, averages were performed over a large enough number of different teachers and training sets to get the precision shown in the figures.

Experiments were carried out for $`N=50`$. The corresponding feature space dimension is $`N(N+1)=2550`$. The restricted feature space considered in our model is composed of the $`N`$ (linear) input components, which define the u-subspace of the feature space, and the $`N\mathrm{\Delta }`$ non redundant quadratic components of the $`\sigma `$-subspace. For the sake of comparison with the theoretical results determined in the thermodynamic limit, we caracterize the actual SVM by its (finite size) inflation factor $`\mathrm{\Delta }=(N+1)/2`$, and the variance $`\sigma ^2`$ of the components in the $`\sigma `$-subspace, related to the normalizing factor $`a`$ of the new features through $`\sigma ^2=Na^2/\mathrm{\Delta }`$. In our case, since $`N=50`$, $`\mathrm{\Delta }=25.5`$ and $`\sigma ^2=1.960784a^2`$, that is $`\sigma ^2=1.960784`$ for the non-normalized mapping and $`\sigma ^2=0.039216`$, for the normalized one.

The values of $`Q`$, the fraction of squared student weights in the $`\sigma `$-subspace, and the teacher-student overlaps $`R_u`$ and $`R_\sigma `$, normalized within the corresponding sub-space, are represented on figures 1 to 4 as a function of $`\alpha M/N`$, using full and open symbols for the mappings $`𝚽_n`$ and $`𝚽_{nn}`$ respectively. Notice that the abscissas correspond to the fraction of training patterns per input space dimension. Error bars are smaller than the symbols’ size. The lines are not fits, but the theoretical curves corresponding to the same classes of teachers as the experimental results. The excellent agreement with the experimental data is striking. Thus, the high order correlations of the features, neglected in the theoretical models, are indeed negligible.

Fig. 1 corresponds to a purely linear teacher ($`Q^{}=0`$), i.e. to a quadratic SVM learning a rule linearly separable in input space. As in this case $`R_\sigma =0`$, only $`R_u`$ and $`Q`$ are represented. In the case of a purely quadratic rule, $`Q^{}=1`$, represented on fig. 2, $`R_u=0`$. Notice that the corresponding overlaps, $`R_u`$ and $`R_\sigma `$, do not have a similar behaviour, as the latter increases much slower than the former, irrespective of the mapping. This happens because, as the number of quadratic components scales like $`N\mathrm{\Delta }`$, a number of examples of the order of $`N\mathrm{\Delta }`$ are needed to learn them. Indeed, $`R_u`$ reaches a value close to $`1`$ with $`\alpha O(1)`$ while $`R_\sigma `$ needs $`\alpha O(\mathrm{\Delta })`$ to reach similar values.

Fig. 3 shows the results corresponding to the isotropic teacher, having $`Q^{}=Q_{iso}^{}\mathrm{\Delta }/(1+\mathrm{\Delta })`$. For $`\mathrm{\Delta }=25.5`$ we have $`Q_{iso}^{}=0.962`$ A particular case of such a teacher has all its weight components of equal absolute value, i.e. $`(w_i^{})^2=1/\stackrel{~}{N}`$, and was studied in  and . Finally, the results corresponding to a general rule, with $`Q^{}=0.5`$, are shown in fig. 4. Notice that at fixed $`\alpha `$, $`R_u`$ decreases and $`R_\sigma `$ increases with $`Q^{}`$ at a rate that depends on the mapping. These quantities determine the student’s generalization error through the combination (19). The fact that they increase as a function of $`\alpha `$ with different speed is a signature of hierarchical learning.

The generalization error $`ϵ_g`$ corresponding to the different rules is plotted against $`\alpha `$ on fig. 5, for both mappings. At any fixed $`\alpha `$, the performance obtained with the normalized mapping is better the smaller the value of $`Q^{}`$. The non-normalized mapping shows the opposite trend: its performance for a purely linear teacher is extremely bad, but it improves for increasing values of $`Q^{}`$ and slightly overrides that of the normalized mapping in the case of a purely quadratic teacher. These results reflect the competition on learning the anisotropically distributed features. In the case of the normalized mapping, the $`\sigma `$-components are compressed ($`\sigma ^2=0.039`$) with respect to the u-components, which have unit variance. This is advantageous whenever the linear components carry the most significant information, which is the case for $`Q^{}1`$. When $`Q^{}=1`$, the linear components only introduce noise that hinders the learning process. As the number of linear components is much smaller than the number of quadratic ones, their pernicious effect should be more conspicuous the smaller the value of $`\mathrm{\Delta }`$. Conversely, the non-normalized mapping has $`\sigma ^2=1.96`$, meaning that the compressed components are those of the u-subspace. Therefore, this mapping is better when most of the information is contained in the $`\sigma `$-subspace, which is the case for teachers with large $`Q^{}`$ and, in particular, with $`Q^{}=1`$.

Finally, for the sake of completeness, the fraction of support vectors $`\rho _{SV}M_{SV}/M`$, where $`M_{SV}`$ is the number of training patterns with maximal stability, is represented on figure 6. This fraction is an upper bound to the generalization error. Notice that these curves present qualitatively the same trends as $`ϵ_g`$. Interestingly, $`\rho _{SV}`$ is smaller for the normalized mapping than for the non-normalized one for most of the rules. Since the student’s weights can be expressed as a linear combination of SVs , this result is of practical interest.

## IV Discussion

In order to understand the results obtained in the previous section, we first analyze the relative behaviour of $`R_u`$ and $`R_\sigma `$, which can be deduced from equation (13). If $`\mathrm{\Delta }_\sigma ^{}\mathrm{\Delta }`$, which is the case for sufficiently small $`Q^{}`$, we get that $`R_\sigma R_u`$. This means that the quadratic components are more difficult to learn than the linear ones. On the other hand, if the teacher lies mainly in the quadratic subspace, $`\mathrm{\Delta }_\sigma ^{}\mathrm{\Delta }`$, and then $`R_\sigma >R_u`$. The crossover between these different behaviours occurs at $`\mathrm{\Delta }_\sigma ^{}=\mathrm{\Delta }`$, for which equation (13) gives $`R_\sigma =R_u`$. For $`N=50`$, which is the case in our simulations, this arises for $`Q_n^{}=0.998`$ or $`Q_{nn}^{}=0.929`$, depending on whether we use the normalized or the non-normalized mapping. In the particular case of the isotropic teacher and the non-normalized mapping, $`Q^{}>Q_{nn}^{}`$, so that $`R_\sigma >R_u`$, as shown on figure 3. These considerations alone are not sufficient to understand the behaviour of the generalization error, which depends on the weighted sum of $`R_\sigma `$ and $`R_u`$ (see equation (19)).

The behaviour at small $`\alpha `$ is useful to understand the onset of hierarchical learning. A close inspection of equations (10-13) shows that in the limit $`\alpha 0`$, $`x=\sigma ^2`$ and $`Q\mathrm{\Delta }\sigma ^2/(\mathrm{\Delta }\sigma ^2+1)`$ to leading order in $`\alpha `$. This results may be understood with the following simple argument: if there is only one training pattern, clearly it is a SV and the student’s weight vector is proportional to it. As a typical example has $`N`$ components of unit length in the u-subspace and $`N\mathrm{\Delta }`$ components of length $`\sigma `$ in the $`\sigma `$-subspace, we have $`Q=N\mathrm{\Delta }\sigma ^2/(N\mathrm{\Delta }\sigma ^2+N)`$. With the normalized mapping, $`lim_{\alpha 0}Q=1/2`$. In the case of the non normalized one $`lim_{\alpha 0}Q=(2\mathrm{\Delta }1)/2\mathrm{\Delta }`$, which depends on the inflation factor of the SVM. In this limit, we obtain:

$`\kappa _{max}`$ $``$ $`{\displaystyle \frac{1+\sigma ^2\mathrm{\Delta }}{\sqrt{1+\sigma ^4\mathrm{\Delta }}}}{\displaystyle \frac{1}{\sqrt{\alpha }}},`$ (21)
$`R_u`$ $``$ $`\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{1}{\sqrt{1+\mathrm{\Delta }_\sigma ^{}}}}\sqrt{\alpha },`$ (22)
$`R_\sigma `$ $``$ $`\sqrt{{\displaystyle \frac{2}{\pi }}}\sqrt{{\displaystyle \frac{\mathrm{\Delta }_\sigma ^{}}{1+\mathrm{\Delta }_\sigma ^{}}}}\sqrt{{\displaystyle \frac{\alpha }{\mathrm{\Delta }}}}.`$ (23)
Therefore, $`R\sqrt{\alpha }`$, like for the simple perceptron MMH , but with a prefactor that depends on the mapping and the teacher.

In our model, we expect that hierarchical learning correspond to a fast increase of $`R`$ at small $`\alpha `$, mainly dominated by the contribution of $`R_u`$. As in the limit $`\alpha 0`$,

$$R\frac{R_u+R_\sigma \sqrt{\sigma ^4\mathrm{\Delta }\mathrm{\Delta }_\sigma ^{}}}{\sqrt{1+\sigma ^4\mathrm{\Delta }}\sqrt{1+\mathrm{\Delta }_\sigma ^{}}},$$
(24)
we expect hierarchical learning if $`\sigma ^4\mathrm{\Delta }1`$ and $`\mathrm{\Delta }_\sigma ^{}1`$. The first condition establishes a constraint on the mapping, which is only satisfied by the normalized one. The second condition, that ensures that $`R_\sigma <R_u`$ holds, gives the range of teachers for which this hierarchical generalization takes place. Under these conditions, $`R`$ grows fast and the contribution of $`R_\sigma `$ is negligible because it is weighted by $`\sqrt{\sigma ^4\mathrm{\Delta }\mathrm{\Delta }_\sigma ^{}}`$. The effect of hierarchical learning is more important the smaller $`\mathrm{\Delta }_\sigma ^{}`$. The most dramatic effect arises for $`Q^{}=0`$, i.e. for a quadratic SVM learning a linearly separable rule.

On the other hand, if $`\sigma ^4\mathrm{\Delta }1`$, which is the case for the non normalized mapping, both $`R_u`$ and $`R_\sigma `$ contribute to $`R`$ with comparable weights. Notice that, if the normalized mapping is used, the condition $`\mathrm{\Delta }_\sigma ^{}1`$ implies that $`Q^{}<Q_{iso}^{}\mathrm{\Delta }/(1+\mathrm{\Delta })`$, where $`Q_{iso}^{}`$ corresponds to the isotropic teacher. A straightforward calculation shows that a fraction of $`47.5\%`$ of teachers satisfies this constraint for $`N=50`$. In fact, the distribution of teachers as a function of $`Q^{}`$ has its maximum at $`Q_{iso}^{}`$. When $`N\mathrm{}`$, the distribution becomes $`\delta (Q^{}Q_{iso}^{})`$, and $`Q_{iso}^{}`$ tends to the median, meaning that in this limit, only about $`50\%`$ of the teachers give raise to hierarchical learning when using the normalized mapping.

In the limit $`\alpha \mathrm{}`$, all the generalization error curves converge to the same asymptotic value as the simple perceptron MMH learning in the feature space, namely $`ϵ_g=0.500489(1+\mathrm{\Delta })/\alpha `$, independently of $`\sigma `$ and $`Q^{}`$. Thus, $`ϵ_g`$ vanishes slower the larger the inflation factor $`\mathrm{\Delta }`$.

Finally, it is worth to point out that for $`\sigma =1`$, which would correspond to a normalizing factor $`a=\sqrt{\mathrm{\Delta }/N}`$, the pattern distribution in feature space is isotropic. Irrespective of $`Q^{}`$, the corresponding generalization error is exactly the same as that of a simple perceptron learning the MMH with isotropically distributed examples in feature space.

Since the inflation factor $`\mathrm{\Delta }`$ of the SVM feature space in our approach is a free parameter, it does not diverge in the thermodynamic limit $`N\mathrm{}`$ . As a consequence, $`ϵ_g`$ does not present any stepwise behaviour, but just a crossover between a fast decrease at small $`\alpha `$ followed by a slower decrease regime at large $`\alpha `$. The results of Dietrich et al.  for the normalized mapping, that corresponds to $`\sigma ^2\mathrm{\Delta }=1`$ in our model, can be deduced by taking appropriately the limits before solving our saddle point equations. The regime where the number of training patterns $`M=\alpha N`$ scales with $`N`$, is straightforward. It is obtained by taking the limit $`\sigma 0`$ and $`\mathrm{\Delta }\mathrm{}`$ keeping $`\sigma ^2\mathrm{\Delta }=1`$ in our equations, with $`\alpha `$ finite. The regime where the number of training patterns $`M=\alpha N`$ scales with $`N\mathrm{\Delta }`$, the number of quadratic features, obtained by keeping $`\stackrel{~}{\alpha }\alpha /(1+\mathrm{\Delta })`$ finite whilst taking, here again, the limit $`\sigma 0`$, $`\mathrm{\Delta }\mathrm{}`$ with $`\sigma ^2\mathrm{\Delta }=1`$. The corresponding curves are represented on figure 7 for the case of an isotropic teacher. In order to make the comparisons with our results at finite $`\mathrm{\Delta }`$, the regime where $`\stackrel{~}{\alpha }`$ is finite is represented as a function of $`\alpha =(1+\mathrm{\Delta })\stackrel{~}{\alpha }`$ using the value of $`\mathrm{\Delta }`$ corresponding to our numerical simulations, namely, $`\mathrm{\Delta }=25.5`$. In the same figure we represented the generalization error $`ϵ_g=(1/\pi )\mathrm{arccos}(R)`$ where $`R`$, given by eq. (19), is obtained after solving the saddle point equations with parameter values $`\sigma ^2=0.039`$ and $`\mathrm{\Delta }=25.5`$.

These results, obtained for quadratic SVMs, are easily generalizable to higher order polynomial SVMs. The corresponding saddle point equations are cumbersome, and will not be given here. We expect a cascade of hierarchical generalization behaviour, in which successively more and more compressed features are learned. This may be understood by considering the set of saddle point equations that generalize equation (13). These equations relate the teacher-student overlaps in the successive subspaces. The sequence of different feature subspaces generalized by the SVM depends on the relative complexity of the teacher and the student. This is contained in the factors $`\mathrm{\Delta }_{\sigma _m}^{}/\mathrm{\Delta }_m`$ corresponding to the $`m^{th}`$ subspace, that appear in the set of equations that generalize eq. (13).

## V Conclusion

We introduced a model that clarifies some aspects of the generalization properties of polynomial Support Vector Machines (SVMs) in high dimensional feature spaces. To this end, we focused on quadratic SVMs. The quadratic features, which are the pairwise products of input components, may be scaled by a normalizing factor. Depending on its value, the generalization error presents very different behaviours in the thermodynamic limit .

In fact, a finite size SVM may be caracterized by two parameters: $`\mathrm{\Delta }`$ and $`\sigma `$. The inflation factor $`\mathrm{\Delta }`$ is the ratio between the quadratic and the linear features dimensions. Thus, it is proportional to the input space dimension $`N`$. The variance $`\sigma `$ of the quadratic features is related to the corresponding normalizing factor. Usually, either $`\sigma 1/\sqrt{N}`$ (normalized mapping) or $`\sigma 1`$ (non normalized mapping). In previous studies, not only the input space dimension diverges in the thermodynamic limit $`N\mathrm{}`$, but also $`\mathrm{\Delta }`$ and $`\sigma `$ are correspondingly scaled.

In our model, neither the proportion of quadratic features $`\mathrm{\Delta }`$ nor their variance $`\sigma `$ are necessarily related to the input space dimension $`N`$. They are considered as parameters caracterizing the SVMs. Since we keep them constant when taking the thermodynamic limit, we can study the learning properties of actual SVMs with finite inflation ratios and normalizing factors, as a function of $`\alpha M/N`$, where $`M`$ is the number of training examples. Our theoretical results were obtained neglecting the correlations among the quadratic features. The agreement between our computer experiments with actual SVMs and the theoretical predictions is excellent. The effect of the correlations does not seem to be important, as there is almost no difference between the theoretical curves and the numerical results.

We find that the generalization error $`ϵ_g`$ depends on the type of rule to be inferred through $`Q^{}`$, the (normalized) sum of the teacher’s squared weight components in the quadratic subspace. If $`Q^{}`$ is small enough, the quadratic components need more patterns to be learned than the linear ones. However, only if the quadratic features are normalized, $`ϵ_g`$ is dominated by the high rate learning of the linear components at small $`\alpha `$. Then, on increasing $`\alpha `$, there is a crossover to a regime where the decrease of $`ϵ_g`$ becomes much slower. The crossover between these two behaviours is smoother for larger values of $`Q^{}`$, and this effect of hierarchical learning disappears for large enough $`Q^{}`$. On the other hand, if the features are not normalized, the contributions of both the linear and the quadratic components to $`ϵ_g`$ are of the same order, and there is no hierarchical learning at all.

In the case of the normalized mapping, if the limits $`\mathrm{\Delta }N\mathrm{}`$ and $`\sigma ^21/N0`$ are taken together with the thermodynamic limit, the hierarchical learning effect gives raise to the two different regimes, corresponding to $`MN`$ or $`MN^2`$, described previously .

It is worth to point out that if the rule to be learned allows for hierarchical learning, the generalization error of the normalized mapping is much smaller than that of the non normalized one. In fact, the teachers corresponding to such rules are those with $`Q^{}Q_{iso}^{}`$, where $`Q_{iso}^{}`$ corresponds to the isotropic teacher, the one having all its weights components equal. For the others, both the normalized mapping and the non normalized one present similar performances. If the weights of the teacher are selected at random on a hypersphere in feature space, the most probable teachers have precisely $`Q^{}=Q_{iso}^{}`$, and the fraction of teachers with $`Q^{}Q_{iso}^{}`$ represent of the order of $`50\%`$ of the inferable rules. Thus, from a practical point of view, without having any prior knowledge about the rule underlying a set of examples, the normalized mapping should be preferred.

## Acknowledgements

It is a pleasure to thank Arnaud Buhot for a careful reading of the manuscript, and Alex Smola for providing us the Quadratic Optimizer for Pattern Recognition program . The experimental results were obtained with the Cray-T3E computer of the CEA (project 532/1999).

SR-G acknowledges economic support from the EU-research contract ARG/B7-3011/94/97.

MBG is member of the CNRS.


format: text

text:
# Exact 𝒪(𝑔²𝛼_𝑠) top decay width from general massive two-loop integrals

## Abstract

We calculate the $`b`$-dependent self-energy of the top quark at $`𝒪(g^2\alpha _s)`$ by using a general massive two-loop algorithm proposed in a previous article. From this we derive by unitarity the $`𝒪(\alpha _s)`$ radiative corrections to the decay width of the top quark, where all effects associated with the $`b`$ quark mass are included without resorting to a mass expansion. Our results agree with the analytical results available for the $`𝒪(\alpha _s)`$ correction to the top quark width.

Precision measurements of the electroweak parameters are a powerful tool for testing the validity of the standard model and searching for new physics. The impressive accuracy attained experimentally, which is also expected to improve in the future, makes a complete two-loop analysis of the data necessary.

However, on the theoretical side, such a complete two-loop electroweak analysis is far from being a simple task. Apart from the proliferation of diagrams, a special kind of technical problem is encountered in electroweak two-loop calculations. These processes involve particles with different masses, and in general need to be evaluated at finite external momenta. It has been known already for some time that general massive two-loop Feynman diagrams, when evaluated at nonvanishing external momenta, lead to complicated and often unknown special functions. See, for instance, ref.  which relates the sunset two-point topology to the Lauricella function.

Recent progress in two-loop electroweak analyses was mainly attained by using mass or momentum expansion methods. We note that in some situations mass expansions turn out to converge well, such as the top mass expansion of the $`\alpha _s`$ corrections to the $`Zb\overline{b}`$ process , while in other cases subleading terms are known to be substantial . In certain situations it was possible to recover the exact function starting from an expansion . On the other hand, when the process under consideration involves more than one ratio of masses, recovering the exact result from an expansion would be difficult, the functions involved being complicated.

In order to circumvent the use of a mass or momentum expansion, we proposed in ref.  a general framework for treating two-loop Feynman graphs by a combination of analytical and numerical methods. Given a Feynman diagram, its tensorial structure is first reduced analytically into a combination of ten special functions $`h_i`$, which are defined by fairly simple integral representations. The result of this analytical reduction is further integrated numerically.

A subset of the ten special functions $`h_i`$ , namely $`h_1`$ and $`h_2`$, are sufficient for treating two-loop corrections in a theory involving only scalars and no derivative couplings. This is the case of radiative corrections of enhanced electroweak strength in the standard model, and in ref.  this method was used for deriving corrections of leading power in the Higgs mass. The resulting radiative corrections agree with independent calculations which use different techniques .

Going from a scalar theory to a full renormalizable theory with fermions and derivative couplings, such as the standard model, increases a lot the complexity of the problem. The analytical reduction of the tensor structure of a graph into $`h_i`$ functions along the lines of ref.  needs to be handled by a symbolic manipulation program such as FORM or Schoonschip, and results in rather lengthy expressions. Once the reduction is performed, the numerical integration methods are the same as those used in ref. , with the only difference that more computing power is needed to handle the number of $`h_i`$ functions which need to be integrated.

In ref.  we have shown how these methods can be used for calculating momentum derivatives of two-loop two-point functions around a finite, non-zero value of the external momentum. Such momentum derivatives are encountered for instance when evaluating wave function renormalization constants on-shell.

In this letter we show how this method can be used for calculating two-point functions of physical interest. We note that, due to the general nature of the tensor reduction algorithm given in ref. , the inclusion of more than two external legs proceeds in the same way. The difference is that the resulting expressions are more complicated than those stemming from two-point functions, and more computing power is needed for performing higher-dimensional numerical integration.

In figure 1 we show the two-loop Feynman graphs relevant for the $`b`$-dependent self-energy of the top quark at order $`g^2\alpha _s`$.

As for the counterterm structure, we only show in figure 1 the counterterm diagrams which are necessary for renormalizing the imaginary part of the self-energy. These are the on-shell one-loop QCD counterterms of the $`b`$-quark mass $`\delta m_b`$ and of the top wave function renormalization constant $`\delta Z_t`$:

$`\delta m_b`$ $`=`$ $`2\alpha _sN_c\pi ^2m_b\left\{{\displaystyle \frac{3}{ϵ}}+{\displaystyle \frac{3}{2}}\left[\gamma +\mathrm{log}\pi +\mathrm{log}\left({\displaystyle \frac{m_b^2}{\mu ^2}}\right)\right]2\right\}`$ (1)
$`\delta Z_t`$ $`=`$ $`\alpha _sN_c\pi ^2\left[{\displaystyle \frac{2}{ϵ}}+\gamma +\mathrm{log}\pi +\mathrm{log}\left({\displaystyle \frac{m_t^2}{\mu ^2}}\right)+2\mathrm{log}\left({\displaystyle \frac{m_t^2}{m_g^2}}\right)4\right]`$
Here, $`\mu `$ is the ’t Hooft mass scale, and $`m_g`$ is the gluon mass infrared regulator. $`N_c=4/3`$ is the color factor.

The imaginary part of the two-loop self-energy on-shell is physical. Writing the top self-energy as $`\mathrm{\Sigma }(p\gamma )=\mathrm{\Sigma }_1(p\gamma )+\gamma _5\mathrm{\Sigma }_{\gamma _5}(p\gamma )`$, the top decay width is given by $`\mathrm{\Gamma }_t=2Im\mathrm{\Sigma }_1(p\gamma =m_t)`$. This can be compared to known analytical and numerical results for the $`𝒪(\alpha _s)`$ QCD corrections to the top decay width , and thus provides a nontrivial test of the two-loop algorithm.

The algebraic reduction of the tensorial structure of the graphs was done by implementing the reduction formulae of ref.  into a symbolic manipulation program. The resulting intermediary decomposition into $`h_i`$ functions is too lengthy to reproduce here. Instead, we give the final results obtained after numerical integration over the analytical decomposition.

Here we would like to discuss briefly the treatment of infrared divergence associated with the massless gluon. To regularize the infrared singularities, we introduce a mass regulator for the gluon. We note that in higher-order QCD calculations introducing a gluon mass would affect the Slavnov-Taylor identities. In such a case, a different approach is needed to treat the infrared singular diagrams. One possibility is to try to extract analytically the IR singular part of the diagram and to treat this by dimensional regularization, while the remaining IR finite part of the diagram can be calculated by numerical integration. Whether or not this can always be done in a simple way is beyond the scope of this article. We also note that our approach is mainly designed to treat massive graphs, while other approaches are available for massless QCD calculations. However, this being an $`𝒪(g^2\alpha _s)`$ correction, the use of a gluon mass regulator is legitimate in this case. The physical quantity $`Im\mathrm{\Sigma }_1(p\gamma =m_t)`$, which is associated with the top quark width, is free of infrared singularities, and we checked numerically that indeed the result becomes independent of the gluon regulator when the mass regulator is much smaller than the $`b`$ mass.

In figure 2 we give numerical results for the ultraviolet finite part of the quantity $`\mathrm{\Sigma }_1(p\gamma =m_t)`$ derived from figure 1; the inclusion of the counterterm contributions shown in figure 1 makes only the imaginary part (which is physical) of the self-energy finite, both UV and IR. We plot the results for a range of the top mass, and assume $`G_F=1.1663710^5`$ GeV<sup>-2</sup>, $`m_W=80.41`$ GeV, $`m_b=4.7`$ GeV, and $`\alpha _s(m_t)=.108`$.

In table 1 we give numerical values for the $`𝒪(\alpha _s)`$ correction to the top decay rate $`tW+b`$ obtained from the imaginary part of the self-energy. Therefore, the $`𝒪(\alpha _s)`$ QCD correction is obtained as an inclusive quantity, integrated over the gluon spectrum. To check the size of the finite mass of the $`b`$ quark, we ran our programs in the vanishing $`b`$ mass limit. At tree level the finite $`b`$ mass amounts to a 3-4 MeV correction to the width, and in the $`𝒪(\alpha _s)`$ correction the effect is negligible. These results agree with existing calculations of $`𝒪(\alpha _s)`$ corrections to $`tW+b`$ , which provides a good test for the two-loop tensor decomposition and numerical integration algorithm.

To conclude, we have shown that the general methods proposed in ref.  can be used to calculate physical radiative corrections. We treated at two-loop order the $`b`$-mass dependent self-energy of the top quark at $`𝒪(g^2\alpha _s)`$. From the imaginary part of the self-energy we extract the $`𝒪(\alpha _s)`$ corrections to the top decay process $`tb+W`$, and find agreement with existing results for this process. The calculation is performed while respecting the actual mass and momentum kinematics of the process, and without resorting to a mass or momentum expansion of the diagrams.

Aknowledgements

The work of A.G. was supported by the US Department of Energy. The work of Y.-P. Y. was supported partly by the US Department of Energy.


text:
# An Abel ordinary differential equation class generalizing known integrable classes

## 1 Introduction

Abel-type differential equations of the first kind are of the form

$$y^{}=f_3y^3+f_2y^2+f_1y+f_0,$$
(1)
where $`yy(x)`$ and $`f_i`$ are arbitrary functions of $`x`$. Abel equations appear in the reduction of order of many second and higher order families , and hence are frequently found in the modelling of real problems in varied areas. A general “exact integration” strategy for these equations was first formulated by Liouville  and is based on the concepts of classes, invariants and the solution of the equivalence problem . Generally speaking, two Abel equations of the first kind belong to the same equivalence class if and only if one can be obtained from the other by means of a transformation of the form

$$\{x=F(t),y=P(t)u+Q(t)\},$$
(2)
where $`t`$ and $`uu(t)`$ are respectively the new independent and dependent variables, and $`F`$, $`P`$ and $`Q`$ are arbitrary functions satisfying $`F^{}0`$ and $`P0`$. By changing variables $`\{x=t,y=\left(g_1u+g_0\right)^1\}`$, where $`\{g_1,g_0\}`$ are arbitrary functions of $`t`$, Abel equations of the first kind can always be written in second-kind format

$$y^{}=\frac{\stackrel{~}{f_3}y^3+\stackrel{~}{f_2}y^2+\stackrel{~}{f_1}y+\stackrel{~}{f_0}}{g_1y+g_0}.$$
(3)
Abel equations of the second kind belong to the same class as their first-kind partners. However, due to the arbitrariness introduced when switching from first to second kind, the transformation preserving the class for the latter becomes

$$\{x=F(t),y=\frac{P_1(t)u+Q_1(t)}{P_2(t)u+Q_2(t)}\},$$
(4)
where $`P_1Q_2P_2Q_10`$. There are infinitely many Abel equation classes, and their classification is performed by means of algebraic expressions invariant under (2), the so called invariants, built with the coefficients $`f_i`$ in (1) and their derivatives.

When the invariants of a given Abel equation are constant, its integration is straightforward: the equation can be transformed into a separable equation as explained in textbooks . On the contrary, when the invariants are non-constant, the integration strategy relies on recognizing the equation as equivalent to one of a set of previously known integrable equations, and then applying the equivalence transformation to that equation. However, for non-constant-invariant Abel equations, only a few integrable classes are known. In a recent work , for instance, a classification of all integrable cases presented by Abel, Liouville and others , including examples from Kamke’s book, showed, in all, only four classes depending on one parameter and seven classes depending on no parameters.

In this work, we present a single multi-parameter Abel equation class (AIA111The acronyms AIA, AIR and AIL are explained below.) generalizing all the integrable classes collected in . In addition, AIA contains as particular case a new subclass (AIR), depending on 6 parameters, all of whose members can be systematically transformed into Riccati-type equations. This AIR class, in turn, includes a 4-parameter fully integrable subclass (AIL). Then we will see in §2.3 that many representatives of subclasses of this AIA class can be mapped into non-trivial Abel equations belonging to different subclasses. Hence, if the member being mapped is solvable, it can be used to generate a different, maybe new, solvable and non-constant-invariant class as will be shown in §3.

Due to its simplicity and the potential preparation of computer algebra routines for solving the related equivalence problem , the material being presented appears to us as a convenient starting point for a more systematic determination of exact solutions for Abel equations.

Another important differential equation problem, complementary to the one discussed in this paper, is the one where a first-order equation $`y^{}=\mathrm{\Phi }(x,y)`$, where $`\mathrm{\Phi }`$ is arbitrary, can be transformed into separable by means of linear transformations (2). This problem is a generalization of the case where $`\mathrm{\Phi }`$ is cubic in $`y`$ and the Abel equation has constant invariant. An algorithmic approach for solving this problem is presented in the precedent paper in this issue.

## 2 The AIL, AIR and AIA classes

As mentioned in the introduction, the whole collection of integrable classes presented in , consisting of four 1-parameter classes and 7 classes without parameters, can be obtained by assigning particular values to the parameters of a single multi-parameter Abel class, AIA. In turn, AIA contains a subclass, AIR, all of whose members can be transformed into Riccati-type equations, and inside AIR there is a fully integrable subclass (AIL) all of whose members can be mapped into first-order linear equations.

Since all these classes are obtained by means of the same procedure, to better illustrate the ideas we discuss first this AIL (Abel, Inverse-Linear) subclass. So, consider the general linear equation $`y^{}+g(x)y+f(x)=0`$, where $`f`$ and $`g`$ are arbitrary, after changing variables by means of the inverse transformation $`\{xy\}`$:222By $`\{xy\}`$ we mean changing variables $`\{x=u,y=t\}`$ followed by renaming $`\{uy,tx\}`$.

$$y^{}=\frac{1}{g(y)x+f(y)}.$$
(5)
An implicit solution to this equation is easily expressed in terms of quadratures as

$$C_1=x\mathrm{exp}\left(g(y)𝑑y\right)+\mathrm{exp}\left(g(y)𝑑y\right)f(y)𝑑y.$$
(6)
The key observation here is that (5) will be of second-kind Abel type for many choices of $`f`$ and $`g`$. For instance, by taking

$$f(y)=\frac{s_0y+r_0}{a_3y^3+a_2y^2+a_1y+a_0},g(y)=\frac{s_1y+r_1}{a_3y^3+a_2y^2+a_1y+a_0},$$
(7)
where $`\{s_1,s_0,r_1,r_0,a_i\}`$ are arbitrary constants, the resulting Abel family is

$$y^{}=\frac{a_3y^3+a_2y^2+a_1y+a_0}{\left(s_1x+s_0\right)y+r_1x+r_0}.$$
(8)
This equation has non-constant invariant, and can be seen as a representative of a non-trivial multi-parameter class, and from its connection with a linear equation, its general solution is obtained directly from Eqs.(6) and (7). As shown in the next subsection, four of the eight parameters in (8) are superfluous, and this number can be reduced further by splitting the class into subclasses. Even so, the class is surprisingly large, including multi-parameter subclasses we have not found elsewhere.333We noted afterwards, however, that (8) could also be obtained using a different approach; for instance by following Olver  and considering (1) as an “inappropriate reduction” of a second-order equation that has a solvable non-Abelian Lie algebra; the resulting restrictions on the coefficients in (1) surprisingly lead to a family of the same class as (8). Among others, the class represented by (8) contains as particular cases the two integrable 1-parameter classes related to Abel’s work , and most of the examples found in Kamke’s as well as in other textbooks .

We note that if instead of starting with a linear equation we were to start with a Bernoulli equation, instead of (8) we would obtain

$$y^{}=\frac{a_3y^3+a_2y^2+a_1y+a_0}{\left(s_1x+s_0x^\lambda \right)y+r_1x+r_0x^\lambda },$$
(9)
which is reducible to (8) by changing $`\{x=t^{1/(1\lambda )},y=u\}`$ followed by a redefinition of the constants $`c_i(1\lambda )a_i`$, and so it belongs to the same class as (8). However, if we start with a Riccati equation instead of a Bernoulli equation, and hence instead of (5) we consider

$$y^{}=\frac{1}{h(y)x^2+g(y)x+f(y)},$$
(10)
and then choose $`f(y)`$, $`g(y)`$ and $`h(y)`$ as in (7), we obtain a 10-parameter $`\{s_i,r_i,a_k\}`$ Abel type family (AIR, meaning Abel, Inverse-Riccati)

$$y^{}=\frac{a_3y^3+a_2y^2+a_1y+a_0}{\left(s_2x^2+s_1x+s_0\right)y+r_2x^2+r_1x+r_0},$$
(11)
which becomes a Riccati-type equation by changing $`\{xy\}`$. Four of these ten parameters are superfluous. The remaining 6-parameter class includes as particular cases the parameterized families presented by Liouville and Appell, as well as other classes depending on no parameters shown in Kamke and  and having solutions in terms of special functions. We note that (8) is also a particular case of (11).

Finally, the same ideas can be used to construct a more general Abel class (AIA, meaning Abel, Inverse-Abel) embracing the previous classes Eqs.(8) and (11) as particular cases. This AIA family is obtained by taking as starting point an Abel equation of the second kind (3) and, taking the coefficients $`f_i`$ of the form (7), arriving at

$$y^{}=\frac{\left(g_1a_3x+g_0a_3\right)y^3+\left(g_1a_2x+g_0a_2\right)y^2+\left(g_1a_1x+g_0a_1\right)y+g_1a_0x+g_0a_0}{\left(s_3x^3+s_2x^2+s_1x+s_0\right)y+r_3x^3+r_2x^2+xr_1+r_0}.$$
(12)
Taking now $`g_1`$ and $`g_0`$ as constants and redefining $`g_1a_ia_i`$, $`g_0a_ib_i`$, instead of (11) we obtain the 16-parameter Abel family

$$y^{}=\frac{\left(a_3x+b_3\right)y^3+\left(a_2x+b_2\right)y^2+\left(a_1x+b_1\right)y+a_0x+b_0}{\left(s_3x^3+s_2x^2+s_1x+s_0\right)y+r_3x^3+r_2x^2+r_1x+r_0}.$$
(13)
Remarkably, by changing $`\{xy\}`$, we arrive at an equation of exactly the same type,

$$y^{}=\frac{\left(s_3x+r_3\right)y^3+\left(s_2x+r_2\right)y^2+\left(s_1x+r_1\right)y+s_0x+r_0}{\left(a_3x^3+a_2x^2+a_1x+a_0\right)y+b_3x^3+b_2x^2+b_1x+b_0}.$$
(14)
Although we are not aware of a method for solving this AIA class for arbitrary values of its 16 parameters, by applying the change of variables $`\{xy\}`$ to representatives of solvable subclasses of AIA, one may obtain representatives of non-trivial new solvable Abel classes (see §2.3). Due to this feature and since AIA already contains as a particular case the AIR class and hence AIL (Eqs.(11) and (8)), this AIA class generalizes in one all the solvable classes presented in , collected from the literature (see §3).

### 2.1 The intrinsic parameter-dependence of the classes AIL and AIR

For the purposes of finding exact solutions to Abel equations, and in so doing, solving the equivalence problem with respect to the classes represented by (8) and to some extent (11), it is relevant to determine on how many parameters these classes intrinsically depend. In fact, the technique we have been using for tackling the equivalence problem relies heavily on the calculation of multivariate resultants, and even for classes depending on just one parameter these calculations require the use of special techniques .

Although the procedure presented below does not exhaust the possible reduction of the number of parameters found in (8) and (11), it suffices to reduce this number by four. The idea is to explore fractional linear transformations of the form

$$\{x=t,y(x)=\frac{au+b}{cu+d}\}.$$
(15)
We note that this transformation does not change the degrees in $`x`$ of the denominators of the right-hand sides of (8) and (11). Keeping that fact in mind, we rewrite the representatives of the AIL and AIR classes as

$$y^{}=\frac{\left(y\alpha _0\right)\left(y\alpha _1\right)\left(y\alpha _2\right)}{G_1(x)y+G_0(x)},$$
(16)
where $`\{G_1(x),G_0(x)\}`$ are polynomials in $`x`$ of the same degree (the ones implied by either (8) or (11)). By changing variables $`\left\{x=t,y=u(t)+\alpha _0\right\}`$, (16) becomes

$$y^{}=\frac{y\left(y\mathrm{\Delta }_{1,0}\right)\left(y\mathrm{\Delta }_{2,0}\right)}{G_1(x)y+\stackrel{~}{G_0}(x)},$$
(17)
where $`\mathrm{\Delta }_{i,0}\alpha _i\alpha _0`$ and $`\stackrel{~}{G_0}(x)`$ is a polynomial in $`x`$ of the same degree as $`G_0(x)`$. Three cases now arise, depending whether three, two or none of the roots $`\alpha _i`$ are equal.

Case 1: $`\alpha _0=\alpha _1=\alpha _2`$

In such a case, (16) already depends on four fewer parameters than (8) or (11); a further change of variables $`\{x=t,y=1/u\}`$ transforms (16) to the form

$$y^{}=\frac{1}{M_1(x)y+M_0(x)},$$
(18)
where $`\{M_1(x),M_0(x)\}`$ are polynomials in $`x`$ of the same degree as $`\{G_1(x),G_0(x)\}`$.

Case 2: $`\alpha _0=\alpha _1\alpha _2`$

In this case, $`\mathrm{\Delta }_{1,0}=0`$, so that changing variables

$$\{x=t,y=\frac{1}{u(t)+\mathrm{\Delta }_{2,0}^{}{}_{}{}^{1}}\}$$

in (17) leads to an equation of the form

$$y^{}=\frac{y}{\stackrel{~}{M_1}(x)y+\stackrel{~}{M_0}(x)},$$
(19)
where $`\{\stackrel{~}{M_1}(x),\stackrel{~}{M_0}(x)\}`$ are polynomials in $`x`$ of the same degree as $`\{G_1(x),G_0(x)\}`$. Note that (19) also depends on four fewer parameters than (8) or (11).

Case 3: $`\alpha _0\alpha _1\alpha _2`$

In this case, $`\mathrm{\Delta }_{1,0}\mathrm{\Delta }_{2,0}0`$, and by changing variables

$$\{x=t,y=\frac{1}{\left(\mathrm{\Delta }_{2,0}^{}{}_{}{}^{1}\mathrm{\Delta }_{1,0}^{}{}_{}{}^{1}\right)u(t)+\mathrm{\Delta }_{1,0}^{}{}_{}{}^{1}}\}$$

in (17) we obtain

$$y^{}=\frac{y(y1)}{\stackrel{~}{\stackrel{~}{M_1}}(x)+\stackrel{~}{\stackrel{~}{M_0}}(x)y},$$
(20)
where $`\{\stackrel{~}{\stackrel{~}{M_1}}(x),\stackrel{~}{\stackrel{~}{M_0}}(x)\}`$ are polynomials in $`x`$ of the same degree as $`\{G_1(x),G_0(x)\}`$ and so (20) too depends on four fewer parameters than either (8) or (11).

### 2.2 Splitting of AIL into cases

Since the AIL class is fully solvable, it makes sense, for the purpose of tackling its related equivalence problem, to consider the further maximal reduction in the number of parameters, and hence completely split the class into a set of non-intersecting subclasses, all of which depend intrinsically on a minimal number of parameters. With this motivation, we performed some algebraic manipulations, finally determining that the AIL class actually consists of two different subclasses, respectively depending on two and one parameters.

To arrive at convenient representatives for these two subclasses of AIL, we start by changing variables444We are assuming $`\omega r_1s_0r_0s_10`$ and $`s_10`$, which is justified, since when $`\omega =0`$ (8) has constant invariant, therefore presenting no interest, and when $`s_1=0`$ it can be transformed into an equation of the form (26) anyway by means of $`\{x=t/r_1,y=(ur_0)/s_0\}`$. in (8) according to

$$\{x=\frac{t}{s_{1}^{}{}_{}{}^{2}},y=\frac{r_0s_{1}^{}{}_{}{}^{2}+r_1u}{\left(s_1s_0+u\right)s_1}\},$$
(21)
and introducing new parameters $`\{k_0,k_1,k_2,k_3\}`$ according to

$`k_0`$ $`=`$ $`{\displaystyle \frac{\left(\left(a_2r_{0}^{}{}_{}{}^{2}+\left(s_0a_0a_1r_0\right)s_0\right)s_0a_3r_{0}^{}{}_{}{}^{3}\right)s_{1}^{}{}_{}{}^{2}}{\left(r_1s_0r_0s_1\right)^2}},`$ (22)
$`k_1`$ $`=`$ $`{\displaystyle \frac{\left(a_2s_13a_3r_1\right)r_{0}^{}{}_{}{}^{2}+\left(\left(2r_1a_22a_1s_1\right)r_0+\left(3a_0s_1a_1r_1\right)s_0\right)s_0}{\left(r_1s_0r_0s_1\right)^2}},`$ (23)
$`k_2`$ $`=`$ $`{\displaystyle \frac{\left(\left(2a_2s_13a_3r_1\right)r_1a_1s_{1}^{}{}_{}{}^{2}\right)r_0+\left(3s_{1}^{}{}_{}{}^{2}a_0+\left(r_1a_22a_1s_1\right)r_1\right)s_0}{s_{1}^{}{}_{}{}^{2}\left(r_1s_0r_0s_1\right)^2}},`$ (24)
$`k_3`$ $`=`$ $`{\displaystyle \frac{a_0s_{1}^{}{}_{}{}^{3}+\left(\left(a_2s_1a_3r_1\right)r_1a_1s_{1}^{}{}_{}{}^{2}\right)r_1}{s_{1}^{}{}_{}{}^{4}\left(r_1s_0r_0s_1\right)^2}}.`$ (25)
Thus (8) is transformed into

$$y^{}=\frac{k_3y^3+k_2y^2+k_1y+k_0}{y+x},$$
(26)
whose solution can be obtained directly from Eqs.(6) and (7) by taking $`s_1=0,r_1=1,s_0=1,r_0=0`$ and $`a_ik_i`$. For classification purposes (see §3), it is convenient to write (26) in first-kind form by changing variables $`\{x=t,y=\frac{1}{u}t\}`$, leading to

$$y^{}=(k_3x^3k_2x^2+k_1xk_0)y^3(3k_3x^22k_2x+k_1+1)y^2+(3k_3xk_2)yk_3.$$
(27)
Two different cases now arise, leading to two non-intersecting subclasses.

Case $`k_30`$

By redefining $`k_3k_{4}^{}{}_{}{}^{2}`$, then changing variables in (26) according to

$$\{x=\frac{k_2+3tk_4}{3k_{4}^{}{}_{}{}^{2}},y=\frac{1}{k_4u}+\frac{k_2+3tk_4}{3k_{4}^{}{}_{}{}^{2}}\},$$
(28)
and next redefining $`\{k_4,k_0,k_1\}`$ in terms of new parameters $`\{\alpha ,\beta ,\gamma \}`$ according to

$$k_4=\frac{\beta }{\gamma },k_0=\frac{k_{2}^{}{}_{}{}^{3}\gamma ^4}{27\beta ^4}\frac{k_2\alpha \gamma ^2}{3\beta ^2}+\gamma ,k_1=\alpha \frac{k_{2}^{}{}_{}{}^{2}\gamma ^2}{3\beta ^2},$$

we obtain a 2-parameter representative of the class, already in first kind format:

$$y^{}=\left(x\alpha \beta x^3\right)y^3+\left(3x^21\alpha \right)y^23xy+1.$$
(29)
Case $`k_3=0`$

This other branch of (26) splits into two subcases: $`k_2=0`$ and $`k_20`$. In the former case (26) becomes constant-invariant, and is thus of no interest. When $`k_20`$, by introducing a new parameter $`\alpha `$ by means of

$$\alpha =k_2k_0\frac{k_{1}^{}{}_{}{}^{2}}{4},$$
(30)
and changing variables in (26) according to

$$\{x=\frac{k_12t}{2k_2},y=\frac{2tk_1}{2k_2}\frac{1}{k_2u}\},$$
(31)
one obtains a simpler representative for this class, depending on just one parameter $`\alpha `$:

$$y^{}=\left(\alpha +x^2\right)y^3\left(2x+1\right)y^2+y.$$
(32)
### 2.3 Generating new integrable classes from solvable members of AIA

The motivation for this work was to try to unify the integrable Abel classes we have seen in the literature. This goal has been partially accomplished with the formulation of the AIL and AIR Abel classes, but there are still other integrable classes, not included in AIR, which however all have the following property: these Abel classes have representatives which could be obtained by changing variables $`\{xy\}`$ in the representative of other Abel classes. In this sense, these representatives are both of Abel and inverse-Abel types, which led us to ask the following question:

> Which Abel classes lead to other Abel classes by applying the inverse transformation $`\{xy\}`$ to one of its representatives?

We have formulated the answer to this question in terms of the following proposition and its corollary.

Proposition: If one Abel equation, $`\alpha `$, maps into another Abel equation, $`\beta `$, by means of the inverse transformation $`\{xy\}`$, then both are of the form (13).

Proof: By hypothesis $`\beta `$ is both of Abel and inverse-Abel type, that is, it is of the form $`y^{}=G(x,y)`$, where $`G`$ is both cubic over linear in $`y`$ (the Abel condition) and linear over cubic in $`x`$ (the inverse-Abel condition). Hence, $`G`$ is a rational function of $`x`$ and $`y`$, with numerator cubic in $`y`$ and linear in $`x`$, and denominator cubic in $`x`$ and linear in $`y`$, and so $`\beta `$ is in fact of the form (13).

Corollary: If a given Abel class is related through the inverse transformation to another Abel class then both classes have representatives of the form (13).

One consequence of this proposition is that AIA, (13), is in fact the most general equation which is both Abel and inverse-Abel; similarly AIR, (11), is the most general which is both Abel and inverse-Riccati; and AIL, (8), is the most general which is both Abel and inverse-linear.

It is worth mentioning here that an Abel class may have many different representatives of the form (13). Consequently, for instance, the AIL integrable class represented by (26), which naturally maps into a linear equation by means of $`\{xy\}`$, also contains members which map into non-trivial Abel equations by means of the same transformation. As an example of this, consider 151 of Kamke’s book:

$$y^{}=\frac{12xy+y^22y^3x}{x^2+1}.$$
(33)
This equation has the form (13) and, by changing variables $`\{xy\}`$, it leads to a non-trivial new Abel class and nevertheless it belongs to the AIL class (see (46)). In other words, the class represented by this equation has a representative of the form (13) different from (33), which, by means of $`\{xy\}`$, maps into a linear equation. An example where the same thing happens with a member of the AIR class is given by

$$y^{}=2xy^2+y^3.$$
(34)
This equation was presented by Liouville ; it is of the form (13) and by means of $`\{xy\}`$ leads to another non-trivial Abel class. In addition (34) belongs to the AIR class, so that it has another representative of the form (13), which by means of $`\{xy\}`$ maps into a Riccati equation. Some other examples illustrating how new solvable classes can be obtained from solvable members of AIA by changing variables $`\{xy\}`$ are shown in the next section.

## 3 AIA: a generalization of known integrable classes

As mentioned in the Introduction, all the solvable classes collected in , four depending on one parameter, labelled A, B, C and D, and another seven not depending on parameters, labelled 1 to 7, are particular members of the class represented by (13). The “classification” of these solvable classes as particular members of AIL, AIR or AIA is as follows:

Although this classification shown in Table 1 is not difficult to verify, it is worthwhile to show it explicitly. Starting with the parameterized class by Abel , shown in  as “Class A”,

$$y^{}=\left(\alpha x+\frac{1}{x}+\frac{1}{x^3}\right)y^3+y^2,$$
(35)
where $`\alpha `$ is the parameter, this equation can be obtained from (27) by taking $`\{k_3=2\alpha ,k_2=1,k_1=1/2,k_0=0\}`$ and changing variables

$$\{x=\frac{t^2}{2},y=2\frac{u+t}{t^3}\}.$$

So, (35) is a member of AIL.

Concerning the 1-parameter class by Liouville , labelled in  as class “B”,

$$y^{}=2\left(x^2\alpha \right)y^3+2\left(x+1\right)y^2,$$
(36)
this equation is obtained from (11) by taking $`\{s_2=0,s_1=0,s_0=1,r_2=1,r_1=0,r_0=0,a_3=0,a_2=0,a_0=2\alpha ,a_1=2\}`$ and changing variables

$$\{x=t,y=\frac{1}{u}t^2\}.$$

(36) is then a member of AIR.

The next solvable class, related to Abel’s work, presented in  as class “C”, is represented by

$$y^{}=\frac{\alpha \left(1x^2\right)y^3}{2x}+\left(\alpha 1\right)y^2\frac{\alpha y}{2x},$$
(37)
and this equation can be obtained from (27) by taking $`\{k_2=0,k_1=\alpha /2,k_0=0k_3=1/2\}`$ and changing variables

$$\{x=\frac{\sqrt{\alpha }}{t},y=\frac{t\left(1tu\right)}{\sqrt{\alpha }}\}.$$

(37) is then a member of AIL.

The last parameterized solvable class shown in , labelled there as class “D”, is related to Appell’s work , and is given by

$$y^{}=\frac{y^3}{x}\frac{\left(\alpha +x^2\right)y^2}{x^2}.$$
(38)
This equation is obtained from (11) by taking $`\{s_2=0,s_1=1,s_0=0,r_2=1,r_1=0,r_0=\alpha ,a_3=0,a_2=0,a_1=0,a_0=1\}`$ and changing variables

$$\{x=t,y=\frac{1}{u}+\frac{\alpha t^2}{t}\}.$$

So, (38) is a member of AIR.

Concerning the classes collected in  not depending on parameters, the first one, there labelled as “Class 1”, was presented by Halphen  in connection with doubly periodic elliptic functions:

$$y^{}=\frac{3y\left(1+y\right)4x}{x\left(8y1\right)}.$$
(39)
This equation, clearly a member of (13), can be obtained by changing $`\{xy\}`$ in

$$y^{}=\frac{y\left(8x1\right)}{3x\left(x+1\right)4y}$$

which in turn can be obtained from (37) (the solvable AIL class) by taking $`\alpha =2/3`$ and changing variables

$$\{x=\frac{(12t)\sqrt{36t}}{9t},y=\frac{12ut\sqrt{36t}}{\left(t+1\right)^2\left(4u3t^23t\right)}\}.$$

This derivation also illustrates how new solvable classes can be obtained by interchanging the roles between dependent and independent variables in solvable members of AIA.

As for the representative of Class 2, by Liouville , shown in  as

$$y^{}=y^32xy^2,$$
(40)
this equation is obtained from (11) by taking $`\{s_2=0,s_1=0,s_0=1,r_2=1,r_1=0,r_0=0,a_3=0,a_2=0,a_1=0,a_0=1\}`$ and converting the resulting equation into first-kind format by changing variables

$$\{x=t,y=\frac{1}{u}t^2\}.$$

(40) is then a member of AIR.

Also of note here, (40) is a special case of the equation presented by Appell in :

$$y^{}=\frac{y^3}{\alpha x^2+\beta x+\gamma }\frac{d}{dx}\left(\frac{ax^2+bx+c}{\alpha x^2+\beta x+\gamma }\right)y^2$$
(41)
with $`\alpha =0,\beta =0,\gamma =1,a=1,b=0,c=0`$. (41) is also seen to be a member of AIR since it can be obtained by changing

$$\{x=t,y=\frac{1}{u}\frac{at^2+bt+c}{\alpha t^2+\beta t+\gamma }\}$$

in (11) with $`s_2=\alpha ,s_1=\beta ,s_0=\gamma ,r_2=a,r_1=b,r_0=c,a_3=0,a_2=0,a_1=0,a_0=1`$.

Class 3, also by Liouville , presented in  as

$$y^{}=\frac{y^3}{4x^2}y^2$$
(42)
can be obtained from (40) by changing $`\{xy\}`$ and then converting it to first-kind format by means of

$$\{x=2t,y=\frac{1}{u}+t\}.$$

(42) is then a member of AIA. This also illustrates the derivation of a solvable class by changing $`\{xy\}`$ in solvable members of the AIR subclass.

The next class, presented in  as Class 4, collected among the Abel equation examples of Kamke’s book,

$$y^{}=y^3\frac{\left(x+1\right)y^2}{x},$$
(43)
can be obtained from (32) by taking $`\alpha =0`$ and changing variables

$$\{x=\frac{1}{t},y=t^2ut\},$$

so it belongs to the AIL subclass.

In , three new integrable classes not depending on parameters were presented too - these are classes “5”, “6” and “7”. Starting with class 5, given by

$$y^{}=\frac{\left(2x+3\right)\left(x+1\right)y^3}{2x^5}+\frac{\left(5x+8\right)y^2}{2x^3},$$
(44)
this equation can be obtained from (8) by taking $`\{s_1=0,s_0=1,r_1=1,r_0=0,a_3=6,a_2=10,a_1=4,a_0=0\}`$ and changing variables

$$\{x=\frac{1}{t},y=\frac{2t^2\left(t+1\right)u}{t\left(\left(t+1\right)u+2t\right)}\}.$$

Regarding Class 6, given by

$$y^{}=\frac{y^3}{x^2\left(x1\right)^2}+\frac{\left(1xx^2\right)y^2}{x^2\left(x1\right)^2},$$
(45)
this equation is obtained from (43) by changing $`\{xy\}`$ and then converting it to first-kind form by means of

$$\{x=t,y=\frac{t1+u}{\left(t1\right)u}\}.$$

Concluding, class 7 is given by

$$y^{}=\frac{\left(4x^4+5x^2+1\right)y^3}{2x^3}+y^2+\frac{\left(14x^2\right)y}{2x\left(x^2+1\right)}.$$
(46)
This equation can be obtained from Kamke’s first order example 151,

$$y^{}=\frac{\left(y^2+1\right)\left(12yx\right)}{x^2+1},$$
(47)
by first changing $`\{xy\}`$ and then converting the resulting equation to first kind format by means of

$$\{x=t,y=\frac{1}{2t}\frac{1}{2\left(t^2+1\right)u}\}.$$

In turn, (47) is a member of the AIL subclass and can be obtained from (37) by taking $`\alpha =4`$ and changing variables

$$\{x=\frac{i}{t},y=\frac{it\left(tu1\right)}{t^2+1}\}.$$

## 4 Conclusions

In this paper, a multi-parameter non-constant-invariant Abel class has been presented which generalizes in one the integrable cases shown in the works by Abel, Liouville, Appell and others, including all those shown in Kamke’s book and other references. This new class splits into various subclasses, many of which are fully integrable, including some not previously shown elsewhere to the best of our knowledge.

The particular subclass represented by (26) mapping non-constant-invariant Abel equations into linear first-order equations contains by itself most of the exactly integrable cases we have seen. The subclass represented by (11) appears to us to be the most general class mapping Abel equations into Riccati ones; indeed it includes the parameterized mappings presented by Liouville and Appell as particular members.

Finally, the mapping of Abel classes into other Abel classes presents a useful way of finding new integrable classes from other classes known to be solvable, as shown in §3. We are presently working on analyzing different connections between all these subclasses and expect to find reportable results in the near future.

Acknowledgments

This work was supported by the Centre of Experimental and Constructive Mathematics, Simon Fraser University, Canada, and by the Symbolic Computation Group, Faculty of Mathematics, University of Waterloo, Ontario, Canada. The authors would like to thank K. von Bülow555Symbolic Computation Group of the Theoretical Physics Department at UERJ - Brazil. for a careful reading of this paper, and T. Kolokolnikov666Department of Mathematics, University of British Columbia - Canada. and A. Wittkopf777Centre for Experimental and Constructive Mathematics, Simon Fraser University - Canada. for useful related discussions. We also acknowledge T. Kolokolnikov and one anonymous referee for their suggestion of normalizing AIR by normalizing the cubic polynomial on the numerator of (8) and (11).


text:
# 1 Introduction

## 1 Introduction

Electromagnetic probes provide a unique tool to investigate the early and interior properties of the system created in a heavy ion collision. Dileptons are of particular interest as their production cross section is related to the correlation function of the iso-vector vector current, the conserved current of the $`\mathrm{SU}_\mathrm{R}(2)\times \mathrm{SU}_\mathrm{L}(2)`$ chiral symmetry. Thus, a careful study of the dilepton spectrum may provide some evidence and signal for the existence of a chirally restored phase created in a heavy ion collisions. Of course, unambiguous proof of chiral restoration requires the measurement of both iso-vector vector and iso-vector axial-vector correlation functions. The detection of the latter, however, is rather difficult, as it does not couple to a penetrating probe but rather to pionic degrees of freedom, and thus is strongly affected by final state interactions. However, at least at low temperatures one can show , that the onset of chiral restoration goes via the mixing of the vector and axial correlators. Assuming that this mixing is the dominant mechanism for chiral restoration (for a discussion of other possibilities see e.g. ), one expects considerable changes in the vector-correlator, which is measurable via dilepton production.

Let us, however, already note at this point that the medium can not only give rise to mixing of vector and axial vector correlation functions in the iso-vector channel, but also to mixing between iso-scalar and iso-vector (‘$`\rho `$ and $`\omega `$’) correlator. This mixing is not related to chiral restoration, but certainly contributes to the dilepton spectrum, as both iso-vector and iso-scalar currents couple to the electromagnetic field. For the purpose of extracting information about possible chiral restoration, it has to be considered as ‘background’.

## 2 Chiral restoration and mixing

As already mentioned in the introduction, to leading order in the temperature the iso-vector vector correlator receives an admixture from the axial correlator . Similar arguments, although somewhat model dependent, can also be given at finite density, where one again finds this mixing of vector and axial vector correlator . In more physical terms, in both cases a pion from either the heat bath (finite temperature) or from the pion could of the nucleons (finite density) couples to the vector-current ($`\rho `$ meson) to form an axial-vector ($`a_1`$ or pion) intermediate state. This is depicted schematically on the left hand side of fig. 1. Conversely, the pions from the heatbath and/or pion cloud also give rise to a vector-admixture to the axial-vector correlator. One, therefore, can imagine that a certain temperature/density vector and axial-vector correlators are fully mixed and thus indistinguishable, which in turn means that chiral symmetry is restored.

The pions in the medium, however, also give rise to other mixing. For example a pion from the heatbath may couple to a $`\rho `$-meson to give an $`\omega `$ intermediate state. This corresponds to the mixing between the iso-vector vector and iso-scalar vector correlators. This process is depicted schematically on the right hand side of fig. 1. This mixing is not related to chiral symmetry restoration, but still affects the dilepton spectrum.

The contribution of these mixings to the dilepton spectrum can be easily understood in simple physical terms. The dilepton-production cross section is proportional to the imaginary part of these correlation functions, which corresponds to a simple cut of the diagrams depicted in fig.1. Therefore, the matrix-elements responsible for vector axial-vector mixing are identical to those for the Dalitz-decay of the $`a_1`$-meson. A similar relation hold between the the mixing of iso-vector and iso-scalar correlators to the Dalitz decay of the $`\omega `$ meson.

In other words, the contribution of the $`a_1`$-Dalitz decay to the total dilepton spectrum provides a measure for the importance of the vector axial-vector mixing, i.e. for the sensitivity to chiral restoration effects. While it might be impossible to extract this channel from the data, in a given model calculation it can be easily accessed. We will show below, that the $`a_1`$-Dalitz contribution is so small that its absence or presence would make an indistinguishable difference in the total spectrum. Consequently, the sensitivity to chiral restoration appears to be very weak.

Fig. 1. Schematic illustration of mixing induced by the medium. Left: Vector axial-vector mixing. Right: Iso-vector iso-scalar mixing.

Finally, there has been discussion about the contributions of baryons, notably the $`N(1520)`$ resonance . Again, this can be understood either as a mixing of the $`\rho `$-meson with $`N(1520)`$-hole states or simply as the Dalitz decay of the $`N(1520)`$. At present, it is unclear what the chiral properties of a $`N(1520)`$-hole state are. Therefore, one cannot say to which extent this mixing contributes to chiral restoration. However, as will be shown below, the contribution of the $`N(1520)`$-Dalitz again is so small, that it hardly affects the total Dilepton spectrum.

## 3 The N(1520)

There has been some discussion as to how important the contribution of baryons to the dilepton spectrum is. Baryon-resonances contribute chiefly through their Dalitz decay $`N^{}e^+e^{}N`$ to the dilepton spectrum. In  this contribution has been estimated to be at most 50 % of that of the $`\omega `$-Dalitz decay in the mass range of $`400500\mathrm{MeV}`$. The $`N(1520)`$ plays a special role, as it couples very strongly to the $`\rho `$ meson and thus should also contribute most to the dilepton spectrum . In order to explore this we have calculated the Dalitz decay rate and width of the $`N(1520)`$ using the model and parameters of , which is based on an analysis of pion photoproduction data. A detailed description of this calculation can be found in . We have also compared to the nonrelativistic models, which are usually employed in the calculation of the $`\rho `$-meson spectral function . In fig. 2 we show the resulting width for two relativistic parameterizations and the nonrelativistic result. They agree reasonably well, and in the transport calculations shown below we employ the nonrelativistic result. We should also mention that we found a rather weak dependence on the off-shell parameter, which is always present in the relativistic description of a spin 3/2 object.

Fig. 2. Dalitz decay width of the N(1520) based on a relativistic and nonrelativistic descriptions. The difference between the calculations labeled ‘Relativistic I’ and ‘Relativistic II’ are discussed in detail in .

## 4 Transport calculation

The results which we will present in the subsequent section are based on a transport calculation. Therefore, it is appropriate to discuss how the dilepton spectrum calculated in transport is related to that obtained from in medium correlation functions (see e.g. ). In general, the dilepton production rate is related to the imaginary part of the current-current correlation function $`C(q_0,\stackrel{}{q})`$

$`{\displaystyle \frac{dR}{d^4q}}={\displaystyle \frac{\alpha ^2}{\pi ^3q^2}}\mathrm{Im}C(q_0,\stackrel{}{q}).`$ (1)
Assuming vector dominance, and concentrating on the iso-vector vector ($`\rho `$) channel the current-current correlator is, up to constant factors, given by the in medium rho propagator $`D_\rho (q_0,\stackrel{}{q})`$

$`\mathrm{Im}C(q_0,\stackrel{}{q})\mathrm{Im}D\rho (q_0,\stackrel{}{q}){\displaystyle \frac{\mathrm{Im}(\mathrm{\Sigma })}{(m^2m_\rho ^2+\mathrm{Re}(\mathrm{\Sigma }))^2+\mathrm{Im}(\mathrm{\Sigma })^2}}`$ (2)
with

$`\mathrm{\Sigma }=\mathrm{\Sigma }_{\rho \pi \pi }+\mathrm{\Sigma }_{\pi +\rho a_1}+\mathrm{}`$ (3)
being the selfenergy of the $`\rho `$ meson. In free space $`\mathrm{\Sigma }=\mathrm{\Sigma }_{\rho \pi \pi }`$. In the medium, one gets additional contributions from all the other channels coupling to the $`\rho `$ such as the $`a_1`$, the $`\omega `$, $`N(1520)`$ etc.

In transport, on the other hand, one folds the collisions of the respective particles with the branching probability into the dilepton channel in order to generate the dilepton spectrum. At first sight, this appears to be a different approach then the one based on the imaginary part of the correlator. However, the relevant cross sections which control the collision probabilities are directly related to the imaginary part of the selfenergies. Thus the transport calculation can also be formulated in terms of an imaginary part of a correlation function which, for the results shown below has the following form

$`\mathrm{Im}C(q_0,\stackrel{}{q})_{transport}\mathrm{Im}D\rho (q_0,\stackrel{}{q})_{transport}{\displaystyle \frac{\mathrm{Im}(\mathrm{\Sigma })}{(m^2m_\rho ^2)^2+\mathrm{Im}(\mathrm{\Sigma }_{\rho \pi \pi })^2}}.`$ (4)
Note that there are two differences between the correlator (2) and that resulting from transport (4). First, in the denominator the real part of the selfenergy in neglected. This is a reasonable approximation as the real part appears to be small . The other approximation is that in the denominator only the free selfenergy enters in the imaginary part. Additional broadening due to other processes is not taken into account. This approximation becomes bad close to the resonance mass, where the imaginary part of the selfenergy dominates the denominator. There, the transport over predicts the actual dilepton yield. However, in the results below we will see that in this region the total dilepton spectrum is dominated by the decay of the $`\omega `$ meson, which happens mostly in the final state. Thus, a reduction of the $`\rho `$ peak will reduce the overall strength only by a small amount. This effect will become more visible, once the experimental mass resolution becomes sufficiently good to resolve the omega peak. Then the collisional broadening of the $`\rho `$ should be directly visible in the neighborhood of the omega. One should also note that this second approximation can in principle be overcome within a transport calculation. The transport gives the full information about the collisional width ($`Im(\mathrm{\Sigma })`$) for each spacetime point. It, however, requires a numerical tour de force to do a statistically reliable calculations. While this has not yet been attempted in the framework of nucleus-nucleus collisions, such a calculation has been carried out for dilepton production in photon-nucleus reactions .

## 5 Results for Pb+Au at 160 GeV/A and predictions for 40 GeV/A

Let us now turn to the actual results. Most of the details of the calculations can be found in . The new element here is the inclusion of the N(1520). In fig. 3 and fig. 4 we show the resulting dilepton spectra in comparison with the preliminary CERES data from 1996 . We find an overall agreement with the data, also for the low transverse momentum ones. Notice, these results have been obtained without any in-medium effects. As explained in  the initial conditions have been chosen such that the final hadronic spectra agree with experiment.

Fig. 3. Results for the dilepton invariant mass spectrum at $`160\mathrm{GeV}/\mathrm{A}`$. All transverse momenta. Calculation based on model of .

As discussed above, an indicator for the chiral mixing of axial-vector and vector correlator is the $`a_1`$ Dalitz contribution (thick full line at the bottom of the graph). In our calculation the $`a_1`$-Dalitz is at best a 2 % contribution. Considering the size of the experimental error bars, it is clear that the present data do not allow for any conclusions concerning chiral restoration. Note that our findings here are also in agreement with refs.  and . Both references find that the $`a_1`$-Dalitz rate is considerably smaller than the pion annihilation rate.

We also find a rather small contribution from the $`N(1520)`$ decay (thick dashed line at the bottom of the graph). It is about a factor of 4 below the contribution of the $`\omega `$-Dalitz decay, supporting the early estimate of . This also implies that at least for the $`160\mathrm{GeV}`$ data, the contribution of the baryons is rather small.

Fig. 4. Results for the dilepton invariant mass spectrum at $`160\mathrm{GeV}/\mathrm{A}`$. Transverse momenta below 500 MeV/c. Calculation based on model of .

Let us next turn to our prediction for the ‘low energy’ ($`40\mathrm{GeV}`$) run. In figs. 5 and 6 we show our prediction together with the 160 GeV data, which are intended as a reference. The prediction is based on central ($`N_{charge}=220`$) events generated with URQMD  (for details of the calculation see ). In our calculation, we have assumed a mass resolution of 1 %, in order to account for the substantially improved mass resolution of the CERES spectrometer. Given this mass resolution, the omega should be clearly visible and thus should put a strong constraint on the model calculations. Aside from the improved resolution, however, we do not predict a significant difference between the high energy and the low energy spectrum, if plotted in the CERES normalization i.e. if divided by the number of charged particles. The only small but visible difference is around 400 MeV, where the yield is somewhat smaller than that of the high energy run. Finally, we find that the contribution of the baryons is still comparatively small, although somewhat stronger than for the high energy collisions.

We have also carried out a calculation in the spirit of , where we have adjusted the initial hadronic state in order to reproduce the final spectra of an RQMD calculation at 40 GeV/A. The results are virtually identical to the ones using the full URQMD events.

Fig. 5. Prediction for the dilepton invariant mass spectrum for Pb+Pb at $`40\mathrm{GeV}/\mathrm{A}`$. All transverse momenta. Calculation based on URQMD . The data are those for 160 GeV/A  and shown for comparison only.

## 6 Conclusions

The imaginary part of the iso-vector vector correlator contributes significantly to the dilepton spectrum. Therefore, a dilepton measurement is in principle sensitive to in medium changes of this correlation function, and thus may be utilized to investigate effects of chiral symmetry restoration, at least indirectly. We have argued that the strength of the $`a_1`$-Dalitz decay provides a good measure for the mixing of the vector and axial-vector correlators, which is one way of restoring chiral symmetry. Our analysis of the presently available dilepton data by the CERES collaboration shows a negligible contribution from the $`a_1`$-Dalitz decay, in agreement with other calculations. Therefore, the present dilepton data do not allow for any conclusions concerning the mixing of the axial-vector and vector correlation functions, and thus about the possible restoration of chiral symmetry achieved in SPS-energy heavy ion collisions.

Fig. 6. Prediction for the dilepton invariant mass spectrum for Pb+Pb at $`40\mathrm{GeV}/\mathrm{A}`$. Transverse momenta below 500 MeV/c. Calculation based on URQMD . The data are those for 160 GeV/A  and shown for comparison only.

We have also provided a prediction for the dilepton spectrum in central Pb+Pb collisions at 40 GeV/A. If normalized by the number of charged particles we find only small differences between the low energy and the high energy spectra. Furthermore, even at the low bombarding energy, we find that the contribution from baryons is rather small.

Acknowledgments: The work of CG and AKDM was supported by the Natural Sciences and Engineering Research Council of Canada and the Fonds FCAR of the Québec Government. The work of CMK was supported in part by the National Science Foundation under Grant No. PHY-9509266 and PHY-9870038, the Welch Foundation under Grant No. A-1358, the Texas Advanced Research Program under Grant FY97-010366-068, and the Alexander v. Humboldt Foundation. MB and VK were supported by the Director, Office of Science, Office of High Energy and Nuclear Physics, Division of Nuclear Physics, and by the Office of Basic Energy Sciences, Division of Nuclear Sciences, of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098. MB was also supported by a Feodor Lynen Fellowship of the Alexander v. Humboldt Foundation, Germany.