text:
# FIRCLA, one-loop correction to 𝑒⁺𝑒⁻→𝜈𝜈̄𝐻 and basis of Feynman integrals in higher dimensions

## 1 INTRODUCTION

Electroweak SM calculations of one-loop corrections to processes with five and more external legs are quite demanding. The number of diagrams to be evaluated typically is in the several thousands, integrands of the individual diagrams are in general very lengthy and in spite of large cancellations which take place the final result usually comes out to be not as compact as required for an efficient and reliable numerical evaluation. At the present time there is no computer algebra system which allows one to perform in an efficient way complete calculations of multi-leg one-loop corrections. For example, FORM is very efficient in doing Lorentz algebra but it can not handle ratios of Gram determinants depending on several momenta and masses. To overcome this kind of problems we use a combination of several systems and exploit their most advantageous features. We are testing the effectiveness of our strategy in the evaluation of the one-loop correction to the process $`e^+e^{}\nu \overline{\nu }H`$. In the last part of this note we also will discuss improvements possible by utilizing non-standard sets of master-integrals.

## 2 FIRCLA

To achieve a better performance we use a combination of several tools like DIANA  based on C and program libraries based on FORM and MAPLE which also act as interfaces between them. Since the main tool of evaluation is based on recurrence relations we called this collection of packages FIRCLA which stands for Feynman Integral Recursive CalcuLAtor.

FIRCLA works as follows. After the process is specified to DIANA it invokes QGRAF , to generate all diagrams, then constructs input expressions suitable for use by FORM and provides additional information like types and masses of particles as well as the relations between kinematical variables. The output from DIANA is stored in a file which is used as an input by our FORM package. Additionally, DIANA may be utilized to produce pictures of all or particular diagrams as a Postscript file.

The FORM package performs the Lorentz algebra, takes traces of the Dirac $`\gamma `$-matrices, transforms products of spinors times $`\gamma `$ matrices to a chosen basis of amplitudes, and utilizing the algorithm , reduces tensor integrals to a combination of scalar integrals, some of them with shifted space-time dimension. For each diagram the FORM package creates a file with expression for further processing by our MAPLE package. The expressions are written in MAPLE format. Each of these files also provides information about scalar products of external momenta and the masses of the particles.

The values of all scalar invariants are calculated in the FORM package. Relations between the momenta carried by the lines of the diagrams ($`p_i`$) and the external momenta ($`q_j`$) are evaluated by DIANA after generating a diagram and then transfered via FORM to the MAPLE package. In this way all useful information from DIANA can be transfered to the MAPLE package. The relations between momenta are needed, for example, to find which integrals are independent or equivalent. In the FORM output files, after the kinematic information an expression in MAPLE format follows. It has the form of a sum of integrals multiplied by polynomials of scalar products, masses and spinor amplitudes. Integrals are just the names of MAPLE procedures like:

$$\mathrm{npoint}(\underset{\mathrm{j}=1}{\overset{\mathrm{N}}{}}\mathrm{P}(\mathrm{j},\mathrm{m}_\mathrm{j}^2)^{\nu _\mathrm{j}},\mathrm{p},\mathrm{data})$$
(1)
where $`P(j,m_j^2)=(k_1p_j)^2m_j^2+iϵ`$ is an inverse scalar propagator, $`\nu _j`$ its power, $`N`$ the number of different propagators, the parameter $`p`$ specifies the shift of the space-time dimension $`D=d+2p`$ and $`data`$ is a set of substitutions for the momenta and scalar invariants.

In the MAPLE package integrals are evaluated separately by using the recurrence relations algorithm for the evaluation of one-loop integrals described in  (see also ).

If the Gram determinants of an integral are different from zero then a set of three relations can be used: the relation for removing dots (a dot represents one power of momentum attached to the line) from lines reads

$`2\mathrm{\Delta }_n\nu _j𝐣^+I_n^{(d)}={\displaystyle \underset{k=1}{\overset{n}{}}}(1+\delta _{jk})\left({\displaystyle \frac{\mathrm{\Delta }_n}{Y_{jk}}}\right)`$
$`\times \left[d{\displaystyle \underset{i=1}{\overset{n}{}}}\nu _i(𝐤^{}𝐢^++1)\right]I_n^{(d)}.`$ (2)
Another relation reduces the shift of the space-time dimension and the index of the $`j`$-th line simultaneously

$`G_{n1}\nu _j𝐣^+I_n^{(d+2)}=`$
$`\left[(_j\mathrm{\Delta }_n)+{\displaystyle \underset{k=1}{\overset{n}{}}}(_j_k\mathrm{\Delta }_n)𝐤^{}\right]I_n^{(d)}.`$ (3)
Integrals with shifted space-time dimension can be expressed in terms of integrals in generic dimension by applying the formula:

$`(d{\displaystyle \underset{i=1}{\overset{n}{}}}\nu _i+1)G_{n1}I_n^{(d+2)}=`$
$`\left[2\mathrm{\Delta }_n+{\displaystyle \underset{k=1}{\overset{n}{}}}(_k\mathrm{\Delta }_n)𝐤^{}\right]I_n^{(d)}.`$ (4)
In the above formulae the shift operators $`𝐣^\pm `$ etc. shift the indices $`\nu _j\nu _j\pm 1`$,

$`\mathrm{\Delta }_n=\left|\begin{array}{cccc}Y_{11}& Y_{12}& \mathrm{}& Y_{1n}\\ Y_{12}& Y_{22}& \mathrm{}& Y_{2n}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ Y_{1n}& Y_{2n}& \mathrm{}& Y_{nn}\end{array}\right|,`$
$`Y_{ij}=(p_ip_j)^2+m_i^2+m_j^2,`$
$`p_j`$ are combinations of external momenta and $`m_j`$ is the mass of the $`j`$-th line , $`_j/m_j^2`$ and

$$G_{n1}=2^n|\begin{array}{cccc}p_1p_1& p_1p_2& \mathrm{}& p_1p_{n1}\\ p_1p_2& p_2p_2& \mathrm{}& p_2p_{n1}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ p_1p_{n1}& p_2p_{n1}& \mathrm{}& p_{n1}p_{n1}\end{array}|.$$

When one of the determinants $`G_{n1}`$ and $`\mathrm{\Delta }_n`$ is equal to zero it is possible to express integrals with $`n`$ lines as a combination of integrals with $`n1`$ lines. If $`G_{n1}=0`$ then the relation

$$I_n^{(d)}=\underset{k=1}{\overset{n}{}}\frac{(_j_k\mathrm{\Delta }_n)}{(_j\mathrm{\Delta }_n)}𝐤^{}I_n^{(d)}$$
(9)
should be applied until one of the lines of the integral will be removed. For the integral with $`n1`$ lines the Gram determinant may be different from zero and therefore the relations (2), (3) and (4) can be used if needed.

If $`\mathrm{\Delta }_n=0`$ one should apply the relation

$$(d\underset{i=1}{\overset{n}{}}\nu _i1)G_{n1}I_n^{(d)}=\underset{k=1}{\overset{n}{}}(_k\mathrm{\Delta }_n)𝐤^{}I_n^{(d2)},$$

until a line of the integral will be removed. As in the previous case if the Gram determinants for integrals with $`n1`$ lines are different from zero one can use relations (2), (3) and (4). Since the space-time dimension of $`n1`$-point integrals will be decreased by this reduction one should subsequently apply the relation

$$I_n^{(d)}=\underset{j=1}{\overset{n}{}}\nu _j𝐣^+I_n^{(d+2)},$$
(10)
in order to increase the dimension back to $`d=42\epsilon `$.

Our MAPLE package reads the FORM output diagram by diagram, evaluates each integral from the input expression, adds results for the integrals and simplifies the sum. For each diagram the result is written in a separate file. The expression stored for each diagram is a combination of master-integrals and spinor-amplitudes multiplied by ratios of polynomials in scalar products of external momenta, masses and the dimension $`d`$. Summing up the results from all diagrams is performed by a separate program.

## 3 $`O(\alpha )`$ CORRECTION TO $`e^+e^{}\overline{\nu }\nu H`$

In order to check the effectiveness of our method of evaluating Feynman diagrams we calculated the one-loop correction to Higgs production via

$$e^+(q_1)e^{}(q_2)\overline{\nu }(q_3)\nu (q_5)H(q_4).$$

This process will be important at future $`e^+e^{}`$ colliders. In the tree approximation it was considered in . At one-loop order recently in .

At $`e^+e^{}`$ linear colliders operating in the 300–800 GeV energy range, the main production mechanisms for SM-like Higgs particles are

$`e^+e^{}`$ $``$ $`(Z)Z+H`$
$`e^+e^{}`$ $``$ $`\overline{\nu }\nu (WW)\overline{\nu }\nu +H`$
$`e^+e^{}`$ $``$ $`e^+e^{}(ZZ)e^+e^{}+H`$
$`e^+e^{}`$ $``$ $`(\gamma ,Z)t\overline{t}+H`$
i.e., “Higgs-strahlung”, $`WW`$-fusion (see figure), $`ZZ`$-fusion, and “radiation off-top”, respectively. The first two processes

are the dominant ones, in particular at energies above 500 GeV.
The Higgs-strahlung cross section for large $`s`$ goes like $`1/s`$ and dominates at low energies. In contrast, the $`WW`$-fusion mechanism exhibits a cross section growing like $`\mathrm{log}(s/M_H^2)`$ and therefore dominates at high energies. At $`\sqrt{s}500`$ GeV, the two processes have approximately the same cross sections. The relevance of one-loop corrections was considered in .

In order to generate all required diagrams a ’technical SM model’ with two leptons doublets (in order to distinguish external leptons from leptons in loops) and one quarks doublet was considered. The contributions of the missing fermions have been obtained by adding up the corresponding results with the appropriate masses.

We put $`m_e=0`$ but keep all other masses $`m_\mu ,m_bm_W,m_Z,m_t,m_H`$ different from zero. The results for the pentagon diagrams depend on 5 scalar invariants. All calculations were done for arbitrary values of the gauge parameters $`\xi _W,\xi _Z,\xi _\gamma `$. In the specified model 326 diagrams contribute to the one-loop correction. Out of these 15 are pentagon type diagrams.

The results have been stored for each diagram separately in a file. Its size is $`100MB`$ ( $`5MB`$ in MAPLE format). The expression for each diagram has the form

$$D_k=\underset{i,j}{}A_{i,j}I_iO_j,$$

where $`I_i`$ are master integrals and $`O_j`$ are 24 spinor amplitudes.

After summing all diagrams we find that the result is gauge invariant as it must be. Each pentagon diagram separately gives a gauge invariant scalar 5-point integral - 15 in total. We observed a huge reduction in the number of scalar n-point (n=4,3,2) integrals $`D_0`$, $`C_0`$,$`B_0`$. For example, from 240 integrals $`C_0`$ only 79 remain in the final answer. It should be mentioned that at the initial stage before taking into account the symmetries of the integrals the total amount of master integrals $`C_0`$ with different values of momenta and masses was more than 500. A special routine has been added to the MAPLE package which recognizes the symmetries of integrals. By permuting masses and momenta into lexicographical order, taking into account momentum conservation, one obtains a set of independent standard integrals and hence a much more compact representation.

The final gauge invariant result is stored in a file of about 1.2 MB in MAPLE format. Thus, it is still rather large for further numerical evaluations. Working out a more compact representation of the result, mainly by optimizing ratios of huge polynomials depending on several variables is one of the problems.

To obtain the complete radiative correction to the observable differential cross section we still have to add the bremstrahlung correction to our result. This work is in progress.

One observation we made is that in the final result the pentagon diagrams yield terms of the form

$$\frac{1}{(d4)}I_{11111}^{(d)}.$$
(11)
In fact, second rank $`n`$-point tensor integrals produce terms of the form $`\frac{1}{(dn+1)}I_{11\mathrm{}1}^{(d)}`$ with a coefficient which is singular at $`d=4`$ for $`n=5`$. In our case, these terms originate from 2nd rank tensor integrals like

$$d^dk_1\frac{(\overline{U}\mathrm{}\widehat{k_1}\mathrm{}U)(\overline{V}\mathrm{}\widehat{k_1}\mathrm{}V)}{c_1c_2c_3c_4c_5}.$$
(12)
The tensor integral here can be written in terms of scalar integrals with shifted dimension:

$`{\displaystyle d^dk_1\frac{k_{1\mu }k_{1\nu }}{c_1c_2c_3c_4c_5}}={\displaystyle \frac{1}{2}}g_{\mu \nu }I_{11111}^{(d+2)}`$
$`+2p_{1\mu }p_{1\nu }I_{31111}^{(d+4)}+2I_{13111}^{(d+4)}p_{2\mu }p_{2\nu }`$
$`+2I_{11311}^{(d+4)}p_{3\mu }p_{3\nu }+2I_{11131}^{(d+4)}p_{4\mu }p_{4\nu }`$
$`+I_{22111}^{(d+4)}\{p_1,p_2\}_{\mu \nu }+I_{21211}^{(d+4)}\{p_1,p_3\}_{\mu \nu }`$
$`+I_{21121}^{(d+4)}\{p_1,p_4\}_{\mu \nu }+I_{12211}^{(d+4)}\{p_2,p_3\}_{\mu \nu }`$
$`+I_{12121}^{(d+4)}\{p_2,p_4\}_{\mu \nu }+I_{11221}^{(d+4)}\{p_3,p_4\}_{\mu \nu }`$
where $`\{p_1,p_2\}_{\mu \nu }=p_{1\mu }p_{2\nu }+p_{1\nu }p_{2\mu }`$.

Different tensor structures give contributions to different of the abovementioned spinor amplitudes. By using the recurrence relations one gets

$$I_{11111}^{(d+2)}=\frac{2\mathrm{\Delta }_5}{(d4)G_4}I_{11111}^{(d)}+\underset{k=1}{\overset{5}{}}\frac{_k\mathrm{\Delta }_5}{(d4)G_4}𝐤^{}I_5^{(d)}.$$

The integral at the left hand side is UV and IR finite, however, in front of the integrals at the right hand side the spurious pole $`1/\epsilon `$ appears. So it looks like we have to evaluate the $`\epsilon `$ term in the expansion of the pentagon integrals. There are no such problems for the 3- and 4-point functions.

In fact multiplying (12) by the Born term, summing over polarizations and taking the traces removes the problematic $`1/(d4)`$ terms. However, if we would attempt to calculate the amplitude numerically before squaring it, we would have to work out first it’s $`\epsilon `$-expansion. A possibility to avoid this problem is to use a different basis of master-integrals in higher dimensions.

## 4 MASTER-INTEGRALS IN HIGHER DIMENSIONS

The idea to express one-loop tensor integrals in terms of integrals with shifted dimension was proposed in  and later in . Recurrence relations which allow us to reduce any one-loop integral with shifted dimension to a combination of integrals in generic dimension were given in  and . These were extended to massless integrals with both Gram determinants zero in .

In  it was discovered that one-loop integrals in higher dimensions provide better numerical stability.

A simple formula expressing any $`n`$-point integral in terms of integrals in higher dimensions was given in . It reads

$$I_n^{(d)}=\frac{(dn+1)G_{n1}}{2\mathrm{\Delta }_n}I_n^{(d+2)}\underset{k=1}{\overset{n}{}}\frac{(_k\mathrm{\Delta }_n)}{2\mathrm{\Delta }_n}𝐤^{}I_n^{(d)}.$$

The improved stability of numerical integrations of integrals in higher dimensions can be seen from their integral representation

$$I_n^{(d)}=\mathrm{\Gamma }\left(n\frac{d}{2}\right)_0^1𝑑x_1\mathrm{}_0^1𝑑x_{n1}J_nh_n^{(d/2n)},$$

where

$`J_2`$ $`=`$ $`1,\mathrm{},J_n=x_{n2}x_{n3}^2\mathrm{}x_1^{n2},`$
$`h_2`$ $`=`$ $`p_{12}^2x_1^2+m_2^2x_1(p_{12}^2m_1^2+m_2^2),`$
$`h_3`$ $`=`$ $`x_1x_2(1x_1)p_{13}^2x_1^2x_2(1x_2)p_{12}^2`$
$`x_1(1x_1)(1x_2)p_{23}^2+x_1x_2m_1^2`$
$`+x_1(1x_2)m_2^2+(1x_1)m_3^2.`$
If we transform integrals to the dimension $`D=d+2n2`$ (assuming that $`d=42\epsilon `$) then the expansion of $`I_n^{(d+2n2)}`$ for small $`\epsilon =2d/2`$ reads

$$I_n^{(d+2n2)}=\mathrm{\Gamma }\left(1+\epsilon \right)\left[\frac{s_n}{\epsilon }s_nR_n^1+O(\epsilon )\right],$$

where

$$s_n=_0^1\mathrm{}_0^1\{dx\}J_nh_n=\frac{p_{ij}^2}{(n+1)!}+\frac{1}{n!}\underset{j=1}{\overset{n}{}}m_j^2,$$

$$R_n^1=_0^1\mathrm{}_0^1\{dx\}J_nh_n\mathrm{ln}h_n.$$

If we transform the integrals to dimension $`D=d+2n4`$ then the expansion of $`I_n^{(d+2n4)}`$ at small $`\epsilon `$ reads

$$I_n^{(d+2n4)}=\mathrm{\Gamma }\left(1+\epsilon \right)\left[\frac{1+\epsilon }{(n1)!\epsilon }R_n^0+O(\epsilon )\right],$$

where

$$R_n^0=_0^1\mathrm{}_0^1\{dx\}J_n\mathrm{ln}h_n.$$

We see that in the integrals $`R_n^{0,1}`$ there are no polynomials in the denominator and this is why these integrals in higher dimensions are more suitable for a direct numerical evaluation.

As an example we give the relations for some integrals in the logarithmic basis. In the latter for $`n>2`$ we express the integrals $`I_n^{(d)}`$ in terms of integrals $`I_n^{(d+2n4)}`$:

$`2\lambda _{123}I_3^{(d)}=(d2)g_{123}I_3^{(d+2)}_1\lambda _{123}I_2^{(d)}(23)`$
$`_2\lambda _{123}I_2^{(d)}(13)_3\lambda _{123}I_2^{(d)}(12),`$
$`4\lambda _{1234}^2\lambda _{123}\lambda _{124}\lambda _{134}\lambda _{234}I_4^{(d)}=(d1)(d3)`$
$`\times g_{123}^2\lambda _{123}\lambda _{124}\lambda _{134}\lambda _{234}I_4^{(d+4)}`$
$`+f_{1234}^{(3)}+f_{2134}^{(3)}+f_{3124}^{(3)}+f_{4123}^{(3)}+f_{1234}^{(2)}`$
$`+f_{3124}^{(2)}+f_{4123}^{(2)}+f_{3214}^{(2)}+f_{4213}^{(2)}+f_{4312}^{(2)}`$
where

$`f_{i_1i_2i_3i_4}^{(2)}=(_{i_1}\lambda _{i_1i_2i_3i_4}_{i_2}\lambda _{i_2i_3i_4}\lambda _{i_1i_3i_4}`$
$`+_{i_2}\lambda _{i_1i_2i_3i_4}_{i_1}\lambda _{i_1i_3i_4}\lambda _{i_2i_3i_4})`$
$`\times \lambda _{i_1i_2i_3i_4}\lambda _{i_1i_2i_3}\lambda _{i_1i_2i_4}I_2^{(d)}(i_3i_4)`$
$`f_{i_1i_2i_3i_4}^{(3)}=_{i_1}\lambda _{i_1i_2i_3i_4}((d2)g_{i_2i_3i_4}\lambda _{i_1i_2i_3i_4}`$
$`+(d3)g_{i_1i_2i_3i_4}\lambda _{i_2i_3i_4})`$
$`\times \lambda _{i_1i_2i_3}\lambda _{i_1i_2i_4}\lambda _{i_1i_3i_4}I_3^{(d+2)}(i_2i_3i_4).`$
where the following notation is used:

$$\lambda _{i_1\mathrm{}i_n}=\mathrm{\Delta }_n,g_{i_1\mathrm{}i_n}=G_{n1}.$$

In these expression we see multiple occurrence of Gram determinants and their derivatives. This fact may lead for some problems for kinematical regions where Gram determinants are close to zero. However outside these regions we expect good numerical stability of these integrals. In fact there is a correspondence between the integrals in higher dimensions and the method advocated in  for a direct numerical calculation of loop integrals. For the one-loop case, (3) and (4) are the explicite solutions of the relations which are derived in  by exploiting the Bernstein-Tkachov theorem.


format: text

text:
# 1 INTRODUCTION

## 1 INTRODUCTION

In refs. , a network based on the parametrical four-wave mixing process (FWM)  that is well-known in nonlinear optics was examined. Such a network is capable to hold and handle information that is encoded in the form of the phase-frequency modulation. In the network the signals propagate along interconnections in the form of quasi-monochromatic pulses at $`q`$ different frequencies

$$\{\omega _l\}_1^q\{\omega _1,\omega _2,\mathrm{},\omega _q\}.$$
$`(1)`$
The model is based on a parametrical neuron that is a cubic nonlinear element capable to transform and generate frequencies in the parametrical FWM- processes $`\omega _i\omega _j+\omega _k\omega _r`$. Schematically this model of a neuron can be assumed as a device that is composed of a summator of input signals, a set of $`q`$ ideal frequency filters $`\{\omega _l\}^q`$, a block comparing the amplitudes of the signals and $`q`$ generators of quasi-monochromatic signals $`\{\omega _l\}^q`$.

Let $`\{K^{(\mu )}\}_1^p`$ be a set of patterns each of which is a set of quasi-monochromatic pulses with frequencies defined by Eq.(1) and amplitudes equal to $`\pm 1`$:

$$K^{(\mu )}=(\kappa _1^{(\mu )},\mathrm{},\kappa _N^{(\mu )}),\text{ where }\kappa _i^{(\mu )}=\pm \mathrm{exp}(ı\omega _{l_i^{(\mu )}}t),\{\begin{array}{c}\mu =1,\mathrm{},p;\hfill \\ i=1,\mathrm{},N;\hfill \\ 1l_i^{(\mu )}q.\hfill \end{array}$$
$`(2)`$
The memory of the network is localized in interconnections $`T_{ij},i,j=1,\mathrm{},N`$, which accumulate the information about the states of $`i`$th and $`j`$th neurons in all the $`p`$ patterns. We suppose that the interconnections are dynamic ones and that they are organized according to the Hebb rule:

$$T_{ij}=(1\delta _{ij})\underset{\mu =1}{\overset{p}{}}\kappa _i^{(\mu )}\kappa _j^{(\mu )},i,j=1,\mathrm{},N.$$
$`(3)`$
The network operates as follows. A quasi-monochromatic pulse with a frequency $`\omega _{l_j}`$ that is propagating along the $`(ij)`$-th interconnection from the $`j`$th neuron to the $`i`$th one, takes part in FWM-processes with the pulses stored in the interconnection,

$$\omega _{l_i^{(\mu )}}\omega _{l_j^{(\mu )}}+\omega _{l_j}\{\omega _l\}_1^q.$$

The amplitudes $`\pm 1`$ have to be multiplied. Summing up the results of these partial transformations over all patterns, $`\mu =1,\mathrm{},p`$, we obtain a packet of quasi-monochromatic pulses, where all the frequencies from the set (1) are present. The amplitudes of the pulses are determined by the interconnection. This packet is the result of transformation of the pulse $`\omega _{l_j}`$ by the interconnection $`T_{ij}`$, and it comes to the $`i`$th neuron. All such packets are summarized in this neuron. The summarized signal propagates through $`q`$ parallel ideal frequency filters. The output signals from the filters are compared with respect to their amplitudes. The signal with the maximal amplitude activates the $`i`$-th neuron (’winner-take-all’). As a result it generates an output signal whose frequency and phase are the same as the frequency and the phase of the activating signal.

Generally, when three pulses interact, under a FWM-process always the fourth pulse appears. The frequency of this pulse is defined by the conservation laws only. However, in order that the abovementioned model works as a memory, an important condition must be add, which has to facilitate the propagation of the useful signal, and, in the same time, to suppress external noise. This condition is the principle of incommensurability of frequencies proponed in ,: no combinations $`\omega _l\omega _l^{}+\omega _{l^{\prime \prime }}`$ can belong to the set (1), when all the frequencies are different. Now we finished to describe the principle of the network operating. This network will be called the parametrical neural network (PNN).

There are arguments going in for PNN. First of all, the frequency-phase modulation is more convenient for optical processing of signals. It allows us to back down an artificial adaptation of an optical network to amplitude modulated signals. Second, when signals with $`q`$ different frequencies can propagate along one interconnection (this is an analog of the channel multiplexing), this, in fact, allows us to reduce the number of interconnections by a factor of $`q^2`$. Note, interconnections occupy nearly 98% of the area of neurochips. Third, the signal-noise analysis made with the aid of the Chebyshev-Chernov statistical method showed that the storage capacity of PNN was approximately $`q^2`$ times as much as the storage capacity of the Hopfield model. Even if $`q10`$, the gain is two orders. For computer processing of colored images the standard value is $`q=256`$. Consequently, comparing with the Hopfield model the gain is about five orders. Simultaneously with an increase of the storage capacity, the noise immunity of the network also increases. For example, we simulated PNN with the following parameters: $`N=100`$, $`q=22`$ and $`p=200`$. This network recognized any 80% noisy pattern after 100 steps (in fact, in one pass over all the neurons). The same network with the parameters $`N=100,q=25,p=1000(!)`$ recognized a 65% noisy pattern in 4-5 passes over all the neurons. We remind that some time ago the ability of the Hopfield model with $`N=400,p=30`$ to recognize a 30% noisy pattern was presented as a high mark in the patterns recognition .

In the present work we investigate the abilities of PNN. Here an important remark has to be done. Generally speaking, there are different parametrical FWM-processes complying with the principle of incommensurability of frequencies. For example in , the parametrical FWM-process of the type

$$\omega _l\omega _l^{}+\omega _{l^{\prime \prime }}=\{\begin{array}{cc}\omega _{l^{\prime \prime }},& \text{ when }l^{}=l;\hfill \\ \omega _l,& \text{ when }l^{}=l^{\prime \prime };\hfill \\ 0& \text{ in other cases.}\hfill \end{array}$$

was examined. The corresponding network will be called PNN-I. However, better results can be obtained for the parametrical FWM-process

$$\omega _l\omega _l^{}+\omega _{l^{\prime \prime }}=\{\begin{array}{cc}\omega _l,& \text{ when}l^{}=l^{\prime \prime };\hfill \\ 0& \text{in other cases.}\hfill \end{array}$$
$`(4)`$
This network will be called PNN-II.

The organization of the paper is as follows. In Section 2 we introduce a vector formalism allowing us to formulate the problem in the general form. In this section the results for PNN-II are presented. In Section 3 the vector formalism is used to examine the Potts-glass neural network. We compare it with PNN-II. Some remarks are given in Conclusions. The details of calculations are in Appendix.

## 2 PNN-II

In fact, PNN is an associative memory of the Hopfield type with neurons, which can be in more than two different states. Such models of neural networks were examined previously (see, for example, ). Usually neurons are modeled with the aid of vectors, but not scalar quantities equal to $`\pm 1`$ or $`0/1`$. The number of the representative vectors is equal to the number of different states of neurons. In the present Section for PNN-II we formulate the vector formalism and then estimate the storage capacity of such a network.

### 2.1 Vector Formalism

In order to describe the $`q`$ different states (1) of neurons we use the set of basis vectors $`𝐞_l`$ in the space $`\mathrm{R}^\mathrm{q},\mathrm{q}1`$,

$$𝐞_l=\left(\begin{array}{c}0\\ \mathrm{}\\ 1\\ \mathrm{}\\ 0\end{array}\right),l=1,\mathrm{},q.$$

The state of the $`i`$th neuron is described by a vector $`𝐱_i`$,

$$𝐱_i=x_i𝐞_{l_i},x_i=\pm 1,𝐞_{l_i}\mathrm{R}^\mathrm{q},1l_iq,i=1,\mathrm{},N.$$

The state of the network as a whole $`X`$ is determined by a set of $`N`$ $`q`$-dimensional vectors $`𝐱_i`$: $`X=(𝐱_1,\mathrm{},𝐱_N)`$. By analogy with Eq.(2) the $`p`$ stored patterns are

$$X^{(\mu )}=(𝐱_1^{(\mu )},𝐱_2^{(\mu )},\mathrm{},𝐱_N^{(\mu )}),𝐱_i^{(\mu )}=x_i^{(\mu )}𝐞_{l_i^{(\mu )}},x_i^{(\mu )}=\pm 1,\{\begin{array}{c}1l_i^{(\mu )}q,\hfill \\ \mu =1,\mathrm{},p.\hfill \end{array}$$

Since in this model neurons are vectors, the local field $`𝐡_i`$ affecting the $`i`$th neuron is a vector too. By analogy with the standard Hopfield model we write

$$𝐡_i=\underset{j=1}{\overset{N}{}}𝐓_{ij}𝐱_j$$
$`(5)`$
The $`(q\times q)`$-matrix $`𝐓_{ij}`$ describes the interconnection between the $`i`$th and the $`j`$th neurons. This matrix affects the vector $`𝐱_j\mathrm{R}^\mathrm{q}`$, converting it in a linear combination of basis vectors $`𝐞_l`$. This combination is an analog of the packet of quasi-monochromatic pulses that come from the $`j`$th neuron to the $`i`$th one after transformation in the interconnection (see Introduction). To satisfy the conditions (3) and (4), we need to take the matrices $`𝐓_{ij}`$ as

$$𝐓_{ij}=(1\delta _{ij})\underset{\mu =1}{\overset{p}{}}𝐱_i^{(\mu )}𝐱_{j}^{(\mu )}{}_{}{}^{+},i,j=1,\mathrm{},N,$$
$`(6)`$
where $`\delta _{ij}`$ is the Kronecker symbol. The elements of these matrices are

$$T_{ij}^{(kl)}=(1\delta _{ij})\underset{\mu =1}{\overset{p}{}}(𝐞_k𝐱_i^{(\mu )})(𝐱_j^{(\mu )}𝐞_l)k,l=1,\mathrm{},q.$$

Let us define the dynamics of our $`q`$-dimensional neurons. Let $`X(t)=(𝐱_1(t),\mathrm{},𝐱_N(t))`$ be the state of the system at the time $`t`$. By definition the $`i`$th neuron at the time $`t+1`$ is oriented along a direction mostly close to the local field $`𝐡_i(t)`$. Let us clarify this definition. With the aid of (6) we write Eq.(5) in the form more convenient for analysis:

$$𝐡_i(t)=\underset{l=1}{\overset{q}{}}A_l^{(i)}𝐞_l,\text{ where }A_l^{(i)}=\underset{j(i)}{\overset{N}{}}\underset{\mu =1}{\overset{p}{}}(𝐞_l𝐱_i^{(\mu )})(𝐱_j^{(\mu )}𝐱_j(t)).$$
$`(7)`$
Let $`k`$ be the index relating to the amplitude that is maximal in modulus in the series (7):

$$\begin{array}{ccc}\hfill A_k^{(i)}=& \mathrm{max}& A_l^{(i)}.\hfill \\ & 1lq& \end{array}$$

Then according to our definition

$$𝐱_i(t+1)=\mathrm{sgn}(A_k^{(i)})𝐞_k.$$
$`(8)`$
The expression (8) is identical to the ’winner-take-all’ rule of Introduction.

The evolution the system consists of consequent changes of orientations of vector-neurons according to the rule (8). We make the convention that if some of the amplitudes are maximal in modulus simultaneously, and the neuron is in one of these unimprovable states, its state does not change. Then it is easy to show that during the evolution of the network its energy $`H(t)=1/2_{i=1}^N(𝐡_i(t)𝐱_i(t))`$ decreases. In the end the system reaches a local energy minimum. In this state all the neurons $`𝐱_i`$ are oriented in an unimprovable manner, and the evolution of the system come to its end. These states are the fixed points of the system. The necessary and sufficient conditions for a configuration $`X`$ to be a fixed point is fulfillment of the set of inequalities:

$$(𝐱_i𝐡_i)(𝐞_l𝐡_i),l=1,\mathrm{},q;i=1,\mathrm{},N.$$
$`(9)`$
### 2.2 Storage capacity of PNN-II

Let us estimate the storage capacity of the network in the limit $`N1`$. Suppose that the network starts from a distorted $`m`$th pattern

$$\stackrel{~}{X}^{(m)}=(a_1\widehat{b}_1𝐱_1^{(m)},a_2\widehat{b}_2𝐱_2^{(m)},\mathrm{},a_N\widehat{b}_N𝐱_N^{(m)}).$$

Here $`\{a_i\}_1^N`$ and $`\{\widehat{b}_i\}_1^N`$ define a phase noise and a frequency noise respectively: $`a_i`$ is a random value that is equal to $`1`$ or $`+1`$ with the probabilities $`a`$ and $`1a`$ respectively; $`b`$ is the probability that the operator $`\widehat{b}_i`$ changes the state of the vector $`𝐱_i^{(m)}=x_i^{(\mu )}𝐞_{l_i^{(\mu )}}`$ , and $`1b`$ is the probability that this vector remains unchanged.

Let us examine to what extent the neural network recognizes the pattern $`X^{(m)}`$ correctly. The amplitudes $`A_l^{(i)}`$ (7) have the form

$$A_l^{(i)}=\{\begin{array}{cccc}\hfill x_i^{(m)}_{j=1}^{N1}\xi _j& +& _{r=1}^L\eta _r(l_i^{(m)}),\hfill & \text{ when }l=l_i^{(m)};\hfill \\ & & _{r=1}^L\eta _r(l),\hfill & \text{ when }ll_i^{(m)},\hfill \end{array}$$
$`(10)`$
where $`\xi _j=a_j(𝐱_j^{(m)}\widehat{b}_j𝐱_j^{(m)})`$, $`\eta _r(l)\eta _j^{(\mu )}(l)=a_j(𝐞_l𝐱_i^{(\mu )})(𝐱_j^{(\mu )}\widehat{b}_j𝐱_j^{(m)})`$, $`j=1,\mathrm{},N`$, $`ji`$, $`\mu =1,\mathrm{},p`$, $`\mu m`$. For simplicity, when writing the quantities $`\eta `$ in place of the superscript $`\mu `$ and the subscript $`j`$ we use the subscript $`r`$ which takes $`L=(N1)(p1)`$ different values $`r=1,\mathrm{},L`$.

Let us note that when the patterns $`\{X^{(\mu )}\}_1^p`$ are uncorrelated, the quantities $`\xi _j`$ and $`\eta _r`$ can be considered as independent random variables described by the probability distributions

$$\xi _j=\{\begin{array}{cccccccccccccccccccc}+1,& & & & & & & & & & & & & & & & & & & \\ 0,& & & & & & & & & & & & & & & & & & & \\ 1& & & & & & & & & & & & & & & & & & & \end{array}\begin{array}{cccccccccccccccccccc}(1a)(1b)& & & & & & & & & & & & & & & & & & & \\ b& & & & & & & & & & & & & & & & & & & \\ (1a)b& & & & & & & & & & & & & & & & & & & \end{array},\eta _r(l)=\{\begin{array}{cccccccccccccccccccc}+1,& & & & & & & & & & & & & & & & & & & \\ 0,& & & & & & & & & & & & & & & & & & & \\ 1& & & & & & & & & & & & & & & & & & & \end{array}\begin{array}{cccccccccccccccccccc}1/2q^2& & & & & & & & & & & & & & & & & & & \\ 11/q^2& & & & & & & & & & & & & & & & & & & \\ 1/2q^2& & & & & & & & & & & & & & & & & & & \end{array}.$$
$`(11)`$
Since the distributions of the quantities $`\eta _r(l)`$ are independent of $`l`$, in what follows we simply write $`\eta _r`$. According to the rule (9), the $`i`$th neuron finds itself in the state $`𝐱_i^{(m)}`$ when two conditions for the amplitudes (10) are fulfilled:

$$\mathrm{sgn}(A_{l_i^{(m)}}^{(i)})=x_i^{(m)},\underset{j=1}{\overset{N1}{}}\xi _j+x_i^{(m)}\underset{r=1}{\overset{L}{}}\eta _r\underset{r=1}{\overset{L}{}}\eta _r.$$

Otherwise there will be an error in the recognition of the vector $`𝐱_i^{(m)}`$. Since the random variable $`x_i^{(m)}\eta _r`$ has the same distribution as $`\eta _r`$, the probability of this error is

$$\underset{i}{\mathrm{Pr}}=\mathrm{Pr}\left\{\underset{j=1}{\overset{N1}{}}\xi _j+\underset{r=1}{\overset{L}{}}\eta _r<0\right\}.$$
$`(12)`$
To estimate the value of $`\mathrm{Pr}_i`$ we use the well-known Chebyshev-Chernov method , (see Appendix). As a result we obtain the expression for the probability of the error in the recognition of the pattern $`X^{(m)}`$:

$$\underset{err}{\mathrm{Pr}}=N\mathrm{exp}\left(\frac{N(12a)^2}{2p}q^2(1b)^2\right)$$
$`(13)`$
When $`N`$ increases, this probability tends to zero, if $`p`$ as function of $`N`$ increases slower than

$$p_c=\frac{N(12a)^2}{2\mathrm{ln}N}q^2(1b)^2$$
$`(14)`$
This allows us to use (14) as an asymptotically possible value of the storage capacity of PNN-II.

When $`q=1`$, Eqs.(13)-(14) transform into well-known results for the standard Hopfield model (in this case there is no frequency noise, $`b=0`$). When $`q`$ increases, the probability of the error (13) decreases exponentially, i.e. the noise immunity of PNN increases noticeably. In the same time the storage capacity of the network increases proportionally to $`q^2`$. In contrast to the Hopfield model the number of the patterns $`p`$ can be much greater than the number of neurons.

For example, let us set a constant value $`\mathrm{Pr}_{err}=0.01`$. In the Hopfield model, with this probability of the error we can recognize any of $`p=N/10`$ patterns, each of which is less then 30% noisy. In the same time, PNN-II with $`q=64`$ allows us to recognize any of $`p=5N`$ patterns with 90% noise, or any of $`p=50N`$ patterns with 65% noise.

In Fig.1 we give an example of the restoration of 90% distorted pattern ($`a=0,b=0.9`$). Here the parameters of the network are $`N=100`$, $`p=200`$, $`q=32`$. The pattern is a picture of a dog. The gray squares are noisy pixels. The states of the network after 50 and 100 steps are shown.

In Fig.2 for different values of $`q`$ we show the dependence of the probability of a pattern recognition $`𝐏_{\mathrm{rec}}=1\mathrm{Pr}_{err}`$ as function of the frequency noise $`𝐛=b100\%,b[0,1]`$, when $`\alpha =p/N=2`$ (solid line); the phase noise is equal to zero, $`a=0`$. We see that if $`q=20`$, we can recognize correctly any pattern when the noise is less than 70%, and if $`q=30`$, any pattern when the noise is less than 85%. Generally, if the noise is less than a critical value $`b_c`$,

$$b_c=1\frac{2}{q}\sqrt{\frac{p}{N}},$$
$`(15)`$
PNN can recognize a noisy pattern for sure, and if $`b>b_c`$, the probability of recognition tends to zero. Our computer simulations confirm these results.

## 3 Potts-glass neural networks

The models of associative memory with neurons that can be in more than two different states have been investigated by a lot of authors . All these models are related with the Potts model of magnetic. The last generalizes the Ising model for the case of the spin variable that takes $`q>2`$ different values ,. In all these works the authors used the same well-known approach relating the Ising model with the Hopfield model (see, for example, ). Namely, in place of the short-range interaction between two nearest spins the Hebb type interconnections between all vector-neurons were used. As a result, long-range interactions appeared. Then in the mean-field approximation it was possible to calculate the statistical sum and, consequently, to construct the phase diagram. Different regions of the phase diagram were interpreted in the terms of the ability of the network to recognize noisy patterns.

Among all the models of $`q`$-state associative memory, characteristics of the anisotropic Potts-glass neural network (APGNN) , , are most close to PNN-II. In other models the storage capacity is less than even for the Hopfield model. Below we describe APGNN in terms of our vector formalism and compare it with PNN-II.

APGNN consists of $`N`$ neurons each of which can be in $`q`$ different states. Now to describe the states of the neurons in place of the basis vectors $`𝐞_l\mathrm{R}^\mathrm{q}`$ (see Subsection 2.1) $`q`$-dimensional vectors of a special type are used. Namely, the $`l`$th state of a neuron is described by a column-vector $`𝐝_l\mathrm{R}^\mathrm{q}`$,

$$𝐝_l=\frac{1}{q}\left(\begin{array}{c}1\\ \mathrm{}\\ q1\\ \mathrm{}\\ 1\end{array}\right),l=1,\mathrm{},q.$$

The state of the $`i`$-th neuron is described by a vector $`𝐱_i=𝐝_{l_i}`$, $`1l_iq`$, $`i=1,\mathrm{},N`$. The state of the network as a whole $`X`$ is determined by a set of $`N`$ $`q`$-dimensional vectors $`𝐱_i`$: $`X=(𝐱_1,\mathrm{},𝐱_N)`$. The $`p`$ stored patterns are

$$X^{(\mu )}=(𝐱_1^{(\mu )},\mathrm{},𝐱_N^{(\mu )}),𝐱_i^{(\mu )}=𝐝_{l_i^{(\mu )}},1l_i^{(\mu )}q,\mu =1,2,\mathrm{},p.$$

The local field $`𝐡_i`$ affecting the $`i`$-th neuron is the vector $`𝐡_i=\underset{j=1}{\overset{N}{}}𝐓_{ij}𝐱_j`$, where $`(q\times q)`$-matrices $`𝐓_{ij}`$ describe the interconnections between the $`i`$-th and the $`j`$-th neurons. These matrices are

$$𝐓_{ij}=(1\delta _{ij})\underset{\mu =1}{\overset{p}{}}𝐱_i^{(\mu )}𝐱_{j}^{(\mu )}{}_{}{}^{+},i,j=1,\mathrm{},N.$$

The same as in Subsection 2.1 the dynamics of APGNN is defined as follows: the $`i`$-th neuron at the next time step $`t+1`$ is oriented along a direction mostly close to the local field $`𝐡_i(t)`$ at the time $`t`$. During the evolution of the network the energy $`H(t)=1/2_{i=1}^N(𝐡_i(t)𝐱_i(t))`$ decreases. The necessary and sufficient conditions for a configuration $`X`$ to be a fixed point is fulfillment of the set of inequalities:

$$(𝐱_i𝐡_i)(𝐝_l𝐡_i),l=1,\mathrm{},q;i=1,\mathrm{},N.$$

We see that PNN-II and APGNN are much alike. The difference between these models is that, first, in APGNN the vectors $`𝐝_l`$ are nonorthogonal, and, second, in APGNN there are no amplitudes $`\pm 1`$ relating with the vectors $`𝐝_l`$.

When $`q=2`$, APGNN is the same as the standard Hopfield model. Repeating the argumentation of Subsection 2.2. we can estimate the storage capacity of APGNN for $`N1`$. We must only take into account that there is no phase noise $`\{a_i\}_1^N`$ in this model. The distorted $`m`$th pattern has a form

$$\stackrel{~}{X}^{(m)}=(\widehat{b}_1𝐱_1^{(m)},\widehat{b}_2𝐱_2^{(m)},\mathrm{},\widehat{b}_N𝐱_N^{(m)}).$$

As above, the random operator $`\widehat{b}_i`$ with the probability $`b`$ changes the state of the vector $`𝐱_i^{(m)}`$, and with the probability $`1b`$ this vector remains unchanged. Now, the probability of the error in the recognition of the vector $`𝐱_i^{(m)}`$ is

$$\underset{i}{\mathrm{Pr}}=\mathrm{Pr}\left\{\underset{j=1}{\overset{N1}{}}\xi _j+\underset{r=1}{\overset{L}{}}\eta _r(l_I^{(m)})<\underset{r=1}{\overset{L}{}}\eta _r(l)\right\},$$

where the independent random variables $`\xi _j=(𝐱_j^{(m)}\widehat{b}_j𝐱_j^{(m)})`$, $`\eta _r=(𝐝_l𝐱_i^{(\mu )})(𝐱_j^{(\mu )}\widehat{b}_j𝐱_j^{(m)})`$ are distributed

$$\xi _j=\{\begin{array}{cc}(q1)/q,& 1b\\ 1/q,& b\end{array},\eta _r=\{\begin{array}{cc}(q1)/q,& 1/q^2\\ 1/q,& (q1)/q^2\\ 0,& (q2)/q\\ 1/q,& (q1)/q^2\\ (q1)/q,& 1/q^2\end{array}$$

Naturally, it is true for the randomized patterns $`\{X^{(\mu )}\}_1^p`$, only.

Similarly to calculations of Appendix we obtain the expression for the probability of the error in the recognition of the pattern $`X^{(m)}`$,

$$\underset{err}{\mathrm{Pr}}=N\mathrm{exp}\left(\frac{N}{2p}\frac{q(q1)}{2}(1\overline{b})^2\right),\overline{b}=\frac{q}{q1}b.$$
$`(16)`$
Then the asymptotically possible value of the storage capacity of APGNN is

$$p_c=\frac{N}{2\mathrm{ln}N}\frac{q(q1)}{2}(1\overline{b})^2.$$
$`(17)`$
When $`q=2`$, these expressions give the known estimates for the Hopfield model. For $`q>2`$ the storage capacity of APGNN is $`q(q1)/2`$ times as large as the storage capacity of the Hopfield model. In  the same factor was obtained by fitting the results of numerical calculations. Our approach allows us to obtain the same result rigorously.

## 4 Conclusions

For $`q>>1`$ the storage capacity of APGNN is two times less than the storage capacity of PNN-II (compare Eq.(17) with Eq.(14) for $`a0`$). When calculating the probability of the error in the recognition, the additional factor two appears in the exponent (see Eqs.(13),(16)). This leads to a significant decrease of a noise immunity of APGNN comparing with PNN-II. This is well seen in Fig.2, where for APGNN the dashed line shows the dependence of the pattern recognition $`𝐏_{\mathrm{rec}}`$ on the value of the frequency noise $`𝐛`$ under the same conditions as for PNN-II (the solid line). The superiority of PNN-II is easily seen, especialy in the region of not so large values of $`q10`$. For APGNN the critical value of the noise $`b_c`$ (15) is less than the analogous characteristic of PNN-II by a quantity $`\frac{0.8}{q}\sqrt{\frac{p}{N}}`$.

In conclusion, we would like to note that our approach allows us to describe not only the optical neural networks of the parametrical type, but also neural networks in which information is encoded in the form of phase delays of pulses in interconnections. It is much more easy to realize such a network in form of a device.

## 5 Appendix

The following equation is true for the probability $`\mathrm{Pr}_i`$ (12):

$$\underset{i}{\mathrm{Pr}}\mathrm{Pr}\left\{\underset{j=1}{\overset{N1}{}}\xi _j+\underset{r=1}{\overset{L}{}}\eta _r0\right\}=\mathrm{Pr}\left\{\underset{j=1}{\overset{N1}{}}\xi _j\underset{r=1}{\overset{L}{}}\eta _r0\right\}$$

Using the known approach of exponential estimates of the Chebyshev type, for any positive $`z>0`$ we obtain:

$$\underset{i}{\mathrm{Pr}}\overline{\mathrm{exp}\left(z\left(\underset{j=1}{\overset{N1}{}}\xi _j\underset{r=1}{\overset{L}{}}\eta _r\right)\right)}=\left(\overline{\mathrm{exp}(z\xi _j)}\left(\overline{\mathrm{exp}(z\eta _r)}\right)^{p1}\right)^{N1}$$

The over-line means an averaging over all possible realizations, and the last equality follows from independence of the random variables $`\xi _j`$ and $`\eta _r`$.

Taking into account the distributions (11), it is easy to obtain the averages

$$\overline{\mathrm{exp}(z\xi _j)}=(1a)(1b)e^z+b+a(1b)e^z,\overline{\mathrm{exp}(z\eta _r)}=e^z/2q^2+11/q^2+e^z/2q^2.$$

Changing the variables $`e^z=x`$ and introducing functions $`f_1(x)`$ and $`f_2(x)`$,

$$f_1(x)=a(1b)x+b+\frac{(1a)(1b)}{x},f_2(x)=\frac{1}{2q^2}\left(x+\frac{1}{x}\right)+1\frac{1}{q^2},$$

we obtain that for any positive $`x`$ the following estimate is valid:

$$\underset{i}{\mathrm{Pr}}\left(f_1(x)f_2^{p1}(x)\right)^{N1}.$$
$`(18)`$
To obtain the minimal possible value of the probability $`\mathrm{Pr}_i`$, we need to find the value of the variable $`x`$ minimizing the right-hand side of Eq.(18). This leads us to the equation

$$(p1)(x^21)+\frac{a(1b)x^2(1a)(1b)}{a(1b)x^2+bx+(1a)(1b)}(x^2+2(q^21)x+1)=0.$$

When $`p1`$, the proper root of this equation up to the terms of the order of $`1/p`$ is equal to $`x_1=1+q^2(12a)(1b)/(p1)`$. Substituting this value of $`x`$ in Eq.(18), we obtain

$$\underset{i}{\mathrm{Pr}}\left(1\frac{q^2(12a)^2(1b)^2}{2(p1)}\right)^{N1}\mathrm{exp}\left(\frac{N(12a)^2}{2p}q^2(1b)^2\right).$$

This inequality gives the estimate (13) for the probability of the error in the recognition of the pattern $`X^{(m)}`$.

List of Figures

Fig.1. The restoration of the pattern with 90% frequency noise ($`b=0.9`$), when $`N=100`$, $`p=200`$ and $`q=32`$. The pattern is a picture of a dog. The gray squares are noisy pixels. The states of the network after 50 and 100 steps are shown.

Fig.2. The probability of the pattern recognition $`𝐏_{\mathrm{rec}}=1\mathrm{Pr}_{err}`$ versus frequency noise $`𝐛=b100\%,b[0,1]`$ for different values of $`q`$ and $`\alpha =p/N=2`$ for PNN-II (solid line) and for the Potts-glass neural network (dashed line).


text:
# The scalar-isoscalar spectral function of strong matter in a large N approximation

## 1 Introduction

With the advent of the new heavy-ion experiments, the significance of the theoretical study of models at high temperature and density is growing. The study started by Hatsuda et al.  of the $`\sigma `$ particle has recently received particular interest. At nonzero temperature $`T`$ and baryon density $`n_B`$ they conjectured the narrowing of the broad $`\sigma `$ resonance: during the restoration of chiral symmetry, the mass of the $`\sigma `$ meson decreases which entails the squeezing of the phase space available for the $`\sigma 2\pi `$ decay. The gradual enhancement of the spectral function near the two-pion threshold means that in spite of the elusive nature of the $`\sigma `$ particle in the vacuum, at $`T=0`$, $`n_B=0`$, there is a chance to see it more clearly in medium.

Our aim is to study in the leading order large N limit of the linear $`\sigma `$ model the trajectory of the sigma pole and its imprint on the shape of the scalar-isoscalar spectral function. Here, we point out only the main features of the investigation, referring for a more comprehensive treatment to Ref.  for the chiral limit and Ref.  for the physical value of the pion mass.

## 2 The model and the quantities of interest

The Lagrangian of the linear $`\sigma `$ model including an explicit symmetry breaking external field $`h`$ is parameterised as:

$`L={\displaystyle \frac{1}{2}}[(\stackrel{}{\varphi })^2m^2\stackrel{}{\varphi }^2]{\displaystyle \frac{\lambda }{24N}}[\stackrel{}{\varphi }^2]^2+\sqrt{N}h\varphi _1.`$

In the broken symmetry phase one separates the vacuum expectation value of the $`\sigma `$ field by the replacement: $`\varphi ^a(\sqrt{N}\mathrm{\Phi }(T)+\varphi ^1,\varphi ^i)`$. For the $`\sigma \pi `$ system, $`N=4`$ and one has $`\mathrm{\Phi }(0)=f_\pi /2`$ where $`f_\pi =93`$ MeV. The role of the external field $`h`$ is to give a mass to the pions: $`m_\pi ^2(T)=h/\mathrm{\Phi }(T)`$.

There are several reasons that favour the use of a large N approximation. First of all since the expansion goes in powers of $`1/N`$, a strongly self-coupled theory does not pose any limitation to this expansion. This approach is insensitive to the choice of the renormalisation point. It leads to a 2$`^{\text{nd}}`$ order chiral transition and provides correct critical description near $`T_c`$. It bears also some nice features required in phenomenology: satisfies at LO the scalar-isoscalar channel unitarity and the Adler-zero condition.

The quantities of interest in our investigation are the equation of state $`\varphi ^1=0`$, the $`\sigma `$ propagator $`G_\sigma (p)`$, whose pole trajectory we are interested in, and its spectral function: $`\rho _\sigma (p_0,T)=\frac{1}{\pi }\underset{\epsilon +0}{lim}\mathrm{Im}G_\sigma (p_0+i\epsilon ,T)`$.
Using the equation of state $`m^2+\frac{\lambda }{6}\mathrm{\Phi }^2(T)+\frac{\lambda }{6N}\varphi ^a\varphi ^a=\frac{h}{\mathrm{\Phi }(T)}`$ in the pion propagator, we check the general Goldstone theorem for pions, i.e. their propagator turns out to be $`G_\pi (p)=1/(p^2m_\pi ^2(T))`$. This propagator is used when calculating the pion bubbles both in the equation of state and in the $`\sigma `$ propagator.

At leading order in the $`1/N`$ expansion $`G_\sigma ^1(p)`$ is given by
$`G_\sigma ^1(p)=p^2m^2[\text{}].`$
With the help of the equation of state, after performing a coupling constant renormalisation $`\frac{1}{\lambda }+\frac{1}{96\pi ^2}\mathrm{ln}\frac{e\mathrm{\Lambda }^2}{M_0^2}=\frac{1}{\lambda _R}`$, one can express $`G_\sigma ^1(p)`$ as

$`G_\sigma ^1(p)=p^2\frac{h}{\mathrm{\Phi }(T)}\frac{\lambda _R\mathrm{\Phi }^2(T)/3}{1\lambda _Rb_R(p)/3},\text{where}b(p)=\text{}.`$

Notice, that the imaginary part of $`G_\sigma `$ is due exclusively to the imaginary pieces in $`b_R(p)`$ ($`b_R(p)`$ is the finite part of $`b(p)`$). In what follows we restrict ourselves to the propagator with zero spatial momentum $`(𝐩=0)`$.

## 3 Analytical continuation

Since we are looking for poles of $`G_\sigma (p_0)`$ on the unphysical, second Riemann sheet, where usually resonances lie, we have to perform the analytical continuation of the expression of $`G_\sigma (p_0)`$. This actually means the continuation of the function $`b_R(p_0)`$, which initially is defined on the physical sheet. We have chosen to perform the continuation of $`b(p_0)`$ across the positive real axis where we have to deal with the discontinuity along the real axis of the original expression for the temperature dependent part of the bubble, $`b_T(p_0)`$ .

The mass $`M_\sigma `$ and the width $`\mathrm{\Gamma }`$ of the $`\sigma `$ particle are the real and the imaginary parts of the solution of $`G_\sigma ^1(p_0)=0`$, respectively. To determine the value of the coupling constant $`\lambda _R`$, we use the phenomenological value of $`M_\sigma /\mathrm{\Gamma }`$ ratio at $`T=0`$ and the restriction imposed by the existence of a tachyonic pole in the theory. The tachyon, which is a pole of $`G_\sigma (p_0)`$ on the imaginary axis of the physical sheet, is intimately related to the effective nature of the $`\lambda \mathrm{\Phi }^4`$ theory. By looking at the coupling renormalisation formula we see that, if one wants to evade triviality ($`\lambda _R=0`$ ) maintaining the stability of the theory (i.e. $`\lambda >0`$) one has to keep the value of the cutoff finite. In order to have a stable effective theory the cutoff has to be below the scale of the tachyon.

The value of the ratio $`M_\sigma /\mathrm{\Gamma }`$ approaches the phenomenologically preferred value with increasing $`\lambda _R`$, but in the same time the tachyon mass comes closer to $`M_\sigma `$. Remaining below the tachyonic scale by a factor of 3 to 5, the closest value is achieved for $`\lambda _R=310`$ ($`M_\sigma =3.5f_\pi `$, $`M_\sigma /\mathrm{\Gamma }1`$) in the limiting chiral case and $`\lambda _R=400`$ ($`M_\sigma =3.95f_\pi `$, $`M_\sigma /\mathrm{\Gamma }=1.4`$) for $`m_\pi (0)=140`$ MeV.

## 4 Results

Fig. 1 shows the basic features of the temperature driven $`\sigma `$ pole trajectory for different values of the $`T=0`$ pion mass. The imaginary part of the pole will eventually decrease with increasing $`T`$ and the pole approaches the two-pion threshold. But in a first stage the imaginary part actually increases, a feature observed also in the context of chiral perturbation theory .

One can also see that with decreasing $`m_\pi (0)`$ the pole approaches more and more the imaginary axis. Below a certain mass value actually two oppositely moving poles emerge from the collision on the imaginary axis of the pole in the 4$`^{\text{th}}`$ quadrant with its 3$`^{\text{rd}}`$ quadrant mirror. The one that moves towards the origin comes down from the imaginary axis switching over to the real axis only when colliding with one pole belonging to the infinite chain of poles which are on the imaginary axis due to the Bose-Einstein distribution function. In the chiral limit no second collision occurs, the pole goes to the origin as $`TT_c`$ along the negative imaginary axis accounting for the observed critical behaviour .

The l.h.s of Fig. 2 presents in more details the trajectory of the $`\sigma `$ pole for the physical pion mass $`m_\pi (0)=140`$ MeV. At $`T^{}0.69m_\pi (0)`$ the real part of the pole goes below the threshold while the imaginary part is still finite. Then at $`T_{\text{real}}`$ the pole collides with its 1$`^{\text{st}}`$ quadrant mirror on the real axis of the 2$`^{\text{nd}}`$ Riemann sheet, below the threshold and splits up in two poles moving along the real axis. At $`T^{}1.07m_\pi (0)`$ one of the poles reaches the threshold, climbs up to the 1$`^{\text{st}}`$ Riemann sheet and becomes physical, the other remains on the real axis of the unphysical sheet. A pole moving on the real axis of the unphysical sheet towards the threshold was recently reported in Ref.  using an alternative approach to the $`\sigma `$ model.

We observe by looking at the spectral function presented on the r.h.s of Fig. 2 that only for $`T<T^{}`$ there is trace of the $`\sigma `$ pole in the spectral function. For higher temperature the approximate Lorentzian shape is distorted and in this way no resonance can be identified. A gradual threshold enhancement occurs in the interval $`T(T^{},T^{})`$, where the real part of the pole is smaller than $`2m_\pi (T)`$. The maximum of the enhancement develops for $`T=T^{}`$, around which $`\rho _\sigma (p_0,T^{})(14m_\pi ^2(T^{})/p_0^2)^{1/2}`$. The details of how this asymptotic behaviour sets in depend crucially on the fact that the pole approaches the threshold along the real axis.

To test the generic nature of pole evolution presented above, we have investigated at $`T=0`$ the dependence of the pole trajectory on the baryon charge density $`n_B`$. Following Ref.  the decrease of the chiral order parameter with increasing values of $`n_B`$ in matter was taken into account by linearly rescaling the $`T=0`$ vacuum expectation value. Qualitatively the same pattern of the pole trajectory was obtained as $`n_B`$ varies, like the one presented in Fig. 2. This shows that the scenario is independent of the specific type of thermodynamical driving force applied to the system.

## 5 Conclusions

In the LO large N approximation of the $`\pi \sigma `$ system we presented a generic scenario for the $`\sigma `$ pole evolution. For a physical value of $`m_\pi (0)`$ the pole hits the real axis of the unphysical sheet below the threshold and produces the threshold enhancement in the spectral function when moving on the real axis towards the threshold. This represents a qualitative prediction for heavy ion collision experiments since the shape of the spectral function in a wide temperature range would lead to a high intensity narrow peak in the $`2\gamma `$ spectra coming from the $`\sigma 2\pi 4\gamma `$ decay chain.


text:
# Search for the Lepton Number Violating Decay 𝜏→𝜇𝛾

## 1 Introduction

The decay $`\tau \mu \gamma `$ is an anticipated lepton-number violating process in supersymmetric models, left-right supersymmetric models and in supersymmetric string unified models. For some ranges of model parameters, decay rates as high several parts per million are expected for this decay, even in light of the current experimental limit on the related $`\mu e\gamma `$ decay. Assuming that neutrinos have mass, the standard model branching ratio is expected at the O($`10^{34}`$) level. The current best published limit is B($`\tau \mu \gamma `$) $`<1.1\times 10^6`$ from the CLEO experiment using 12.6 million $`\tau `$ pairs.

This analysis uses $`56\times 10^6`$ $`\tau ^+\tau ^{}`$ events from the 1999-2001 BABAR data set to give a preliminary result of B($`\tau \mu \gamma `$). The general signature of the $`\tau \mu \gamma `$ signal is the presence of an isolated $`\mu `$ and $`\gamma `$ which have an invariant mass consistent with that of the $`\tau `$(1.777 GeV) and with an energy in the centre-of-mass system (CMS) consistent with the beam energy in the CMS (5.29 GeV) and the rest of the particles in the event having properties that are consistent with being produced in a generic 1-prong $`\tau `$ decay111A ‘1-prong’ decay contains a single charged particle amongst its decay products.. The analysis is performed using a blinded mass-energy region in the data which corresponds to a 3$`\sigma `$ error ellipse centred on the expected peak of the distribution.

A cut-based approach is used, in which the background in the signal-box is estimated from extrapolations using side-band data and verified with the Monte Carlo simulation. The determination of selection criteria in this analysis was based on Monte Carlo simulation of background and signal as well as on data in side-band regions further away from the signal-box.

Background sources arising from non-$`\tau ^+\tau ^{}`$ sources, the most problematic of which are radiative $`\mu ^+\mu ^{}`$ pair events, are reduced to small levels. This leaves the generic 1-prong $`\tau `$ decays as the major source of background. The dominant, and irreducible, background is from radiative muonic decays of the $`\tau `$ in which the two neutrinos have very little energy.

## 2 Monte Carlo Samples

The Monte Carlo simulation uses a complete description of the BABAR detector response employing the GEANT4 software. The most important simulation samples are those involving $`e^+e^{}\tau ^+\tau ^{}`$ decays. These include 27.6 million generic $`\tau ^+\tau ^{}`$ events, which use KORALB employing the TAUOLA decay package, and the $`\tau \mu \gamma `$ signal incorporated into TAUOLA. Forty thousand signal events were generated which form the bulk of the signal Monte Carlo used for the analysis. A second signal sample with differing beam background contributions was produced in order to provide a cross check on the sensitivity of the selection to different background conditions.

Additionally, simulated $`\mu ^+\mu ^{}`$ events, which uses the AfkQed generator, provide tools for $`\mu `$ particle identification and some guidance on $`e^+e^{}\mu ^+\mu ^{}`$ event suppression. Other potential backgrounds include $`e^+e^{}u\overline{u};d\overline{d};s\overline{s},e^+e^{}c\overline{c}`$ and $`e^+e^{}b\overline{b}`$ which were studied using the corresponding generic BABAR simulated event samples. As the selection is restricted to events in which both $`\tau `$’s decay via 1-prong modes, the $`e^+e^{}b\overline{b}`$ background is completely negligible as indicated by the simulated events. The $`e^+e^{}c\overline{c}`$ and light-quark background contributions are also found to be very small, with no events in the simulation surviving from these two sources after all cuts have been applied, and consequently have not been included further in the data-Monte Carlo comparisons. The two-photon processes are negligible backgrounds for events in which the beam energy is detected in a CMS hemisphere.

## 3 Mass and Energy Determination

Because there is no missing energy in the $`\tau \mu \gamma `$ signal, the mass of the measured decay products is the $`\tau `$ mass and the energy is the full energy of the $`\tau `$. Ignoring the effects of initial-state radiation, the energy of the final state particles is equal to the full energy of the beam in the CMS. This situation lends itself very well to the use of a mass of the $`\mu \gamma `$ system calculated from the kinematic fit employing the beam energy constraint on the energy of the system as used by ARGUS, and denoted as $`m_{EC}`$. The energy variable used is the difference between the measured $`\mu \gamma `$ energy and the beam energy, denoted by the symbol $`\mathrm{\Delta }E`$. The distributions of these variables are shown in Figure 1(a)-(c). The resolutions of the core of these distributions, which represent those events with well reconstructed photons and tracks, are $`\sigma _{\mathrm{\Delta }E}`$ = 88 MeV and $`\sigma _{m_{EC}}`$ = 19 MeV.

However, as is evident by the diagonal band of events present in Figure 1(c), when initial state radiation shifts the true energy of the $`\tau `$ from the beam energy, there a negative correlation between $`\mathrm{\Delta }E`$ and $`m_{EC}`$ and a substantial loss of resolution. The low-energy tail of the $`\mathrm{\Delta }E`$ distribution for well reconstructed $`m_{EC}`$ , evident in Figure 1(c), is populated by events with a photon reconstructed with the correct direction, but with significant energy loss.

Figure 2 shows the various regions in the $`m_{EC}`$ vs $`\mathrm{\Delta }E`$ plane used for determining the selection criteria, measuring the background, and defining the signal and blinded regions. A ‘Grand Side-Band’ is defined as the region within the $`1.0<\mathrm{\Delta }E<0.5`$ GeV and $`1.5<m_{EC}<2.1`$ GeV bounds. It is used for evaluating the reliability of the estimation of the signal efficiency and this is also shown on the figure.

The signal-box is an elliptical region centred on the peak of the two-dimensional distribution as determined by the Monte Carlo. For the positive side of $`\mathrm{\Delta }E`$, the ellipse has a 3$`\sigma `$ half-axis for both the $`\mathrm{\Delta }E`$ and $`m_{EC}`$ axes whereas for the negative side of $`\mathrm{\Delta }E`$ the ellipse has a 2$`\sigma `$ half-axis in $`\mathrm{\Delta }E`$ and a 3$`\sigma `$ half-axis for $`m_{EC}`$. The resolutions used in defining the signal-box are those of the core Gaussians obtained from the signal Monte Carlo. This asymmetric shape provides an optimal signal-box given the resolution in $`\mathrm{\Delta }E`$ and the presence of increasing background in the negative $`\mathrm{\Delta }E`$ region that is not present for positive $`\mathrm{\Delta }E`$.

Events that are ultimately selected and included in a study of the $`m_{EC}`$ and $`\mathrm{\Delta }E`$ variables must first pass a number of selection criteria. These requirements are designed primarily to reduce the number of Bhabha and $`\mu ^+\mu ^{}`$ events. The requirements of the standard BABAR background filter which the signal events pass:

* 2 charged tracks
* zero net charge
* $`\left|p_1\right|+\left|p_2\right|<9.0`$ GeV
* $`E_1+E_2<5.0`$ GeV
* $`E_1/\left|p_1\right|<0.8`$ OR $`E_2/\left|p_2\right|<0.8`$
* $`E_{CM}\left(\left|p_1\right|_{CM}+\left|p_2\right|_{CM}\right)>0`$
* $`p_{t\left(CM\right)}/\left(E_{CM}\left(\left|p_1\right|_{CM}+\left|p_2\right|_{CM}\right)\right)>0.07`$
* charged track separation $`>90^{}`$ in the CMS
* at least one $`\gamma `$
* fiducial region: $`0.775<\mathrm{cos}\theta _{track}<0.940`$

where $`p_{1(2)}`$ is the track momentum; $`E_{1(2)}`$ is the calorimeter energy associated with the track; $`p_{t(CM)}`$ is the total transverse momentum of charged tracks in the CMS; and subscript CM indicates when the quantity is in the CMS.

These requirements have an efficiency of $`37.3\pm 0.2\%`$. These efficiency losses arise largely from the fiducial acceptance of the detector, but also from the branching ratio to 1-prong modes. At this stage of the analysis, $`22\times 10^6`$ events are selected in the data. Monte Carlo studies indicate that 83% of the signal events pass these requirements for events which would have otherwise been selected.

Subsequent selections employ requirements to remove the non-$`\tau `$ background sources and the $`e^+e^{}\mu ^+\mu ^{}\gamma `$ and $`e^+e^{}e^+e^{}`$ events in particular. Cuts include a specific requirement on the observed CMS energy, $`E_{vis}^{tot}/2E_{Beam}^{CMS}<0.95`$, to exclude the aforementioned background sources in general and a $`\mu `$ veto for reduction of di-muon backgrounds in particular. The event is tagged via an electron or $`h1\pi ^0`$ on the opposing side. The missing mass on the tag side is restricted to reduce backgrounds from lost or difficlt to reconstruct tracks. To be consistent with a decay through an intermediate $`\tau `$, the tag track is required to have less than 80% of the CMS beam energy. In general, this will not be the case for the non-$`\tau `$ backgrounds.

On the signal side, we impose the requirement of an identified $`\mu `$ and an associated $`\gamma `$ with energy in excess of 400 MeV. These cuts are applied in addition to the obvious $`m_{EC}`$ and $`\mathrm{\Delta }E`$ cuts for the final selection.

The cut progression shows that once the non-$`\tau `$ background is removed by requiring an electron-vs-$`\mu \gamma `$ or $`h1\pi ^0`$-vs-$`\mu \gamma `$ events, the $`\tau `$ Monte Carlo simulation tracks the data reasonably well.

We present in Figures 3 and 4 a sample of the distributions of variables employed in the analysis after all other requirements, apart from that using the one shown, have been applied. The data distributions are represented by points and the $`\tau `$ Monte Carlo simulation by histograms in these figures. The normalizations of the Monte Carlo distributions is fixed by the luminosity.

These distributions are well described by the $`\tau ^+\tau ^{}`$ simulation in both shape and overall normalization. There are 604 events in the data in these distributions and the ratio of data to the Monte Carlo expectation is $`1.022\pm 0.069(stat)\pm 0.025(norm)`$. The first error is statistical and the second is the systematic error associated with normalization of the Monte Carlo events. This normalization uncertainty is a component of the systematic error associated with the efficiency of selecting the signal. This includes uncertainties in integrated luminosity, cross-section, radiation treatment in the generator and branching ratio uncertainties.

Because the signal has a relatively stiff momentum spectrum, muon particle identification ($`\mu `$-PID) efficiency using a sample of radiative $`\mu `$ pair events ($`e^+e^{}\mu ^+\mu ^{}\gamma `$) from the same data set used for the rest of the analysis was exploited to study the $`\mu `$-PID efficiency. Contamination from sources other than the $`\mu `$ pair, such as $`\tau ^+\tau ^{}`$ events, are negligible in this sample.

The ratio of the efficiencies as a function of LAB momentum of the data to that of the $`e^+e^{}\mu ^+\mu ^{}(\gamma )`$ simulation was studied and a correction based upon the control sample was extracted. Similarly obtained is the correction for the $`\tau ^+\tau ^{}`$ events in the Monte Carlo sample and the signal Monte Carlo. This method for applying the correction is particularly appropriate here as the polar-angle, azimuthal and momentum distributions and their correlations are very similar between the $`\mu `$ pair control sample and the signal and background. From these studies, an additional correction factor of 0.80 to 0.85 is applied to the efficiency as predicted by the Monte Carlo.

## 4 Backgrounds

There are two broad categories of background that require suppression: the background arising from non-$`\tau ^+\tau ^{}`$ sources, the most significant being $`e^+e^{}\mu ^+\mu ^{}(\gamma )`$, and those coming from standard $`\tau `$ decays. In the latter category, the $`\tau \mu \nu \nu \gamma `$ form an irreducible background.

The $`\tau `$ backgrounds are supressed by requiring the signal $`\gamma `$ to be energetic, $`E_\gamma >40MeV`$, and by imposing tight requirements on the $`\mu `$-PID of the track.

The non-$`\tau `$ background is reduced by requiring there to be a non-$`\mu `$ tag on the non-signal hemi-sphere and by removing events with measured momentum and calorimeter energy characteristic of events and hemispheres with little or no undetected energy. These events are particularly problematic as they naturally possess a value of $`\mathrm{\Delta }E`$ that overlaps the signal.

After all the selection criteria have been applied, before applying the signal-box cut, there remain background sources in the Grand Side-Band region that, from the Monte Carlo simulation, consist of $`\tau \mu \nu \overline{\nu }`$ (85.9%), $`\tau \pi \nu `$ and $`\tau K\nu `$ (10.6%) and $`\tau \rho \nu `$ (3.5%). The non-$`\tau `$ background is expected to be very small in the final sample.

The background is estimated by fitting the side-bands in $`m_{EC}`$ on a sample of data selected to have $`2\sigma _{\mathrm{\Delta }E}<\mathrm{\Delta }E<3\sigma _{\mathrm{\Delta }E}`$. As seen in Figure 6, the distribution of the background is reasonably uniform in the region of the final selection, which enables one to use a simple linear interpolation to estimate the density of background one expects in the signal-box region. A geometrical correction is applied to estimate the background contained in the elliptical signal-box.

A low statistics check of the method is made by applying it to the $`\tau ^+\tau ^{}`$ Monte Carlo sample. In this case, the luminosity scaled number of events predicted from the side-bands is $`9.1\pm 2.1`$ events, in good agreement with the $`5.5\pm 3.2`$ events observed. This estimate of $`5.5\pm 3.2`$ is the only background estimate based purely on counting Monte Carlo events.

The data side-band measurements yield estimates of $`7.8\pm 1.4`$ events.

## 5 Systematics

Using the $`\tau ^+\tau ^{}`$ Monte Carlo, the above selections yield an absolute efficiency for the signal of $`6.7\pm 0.1(stat)\%`$. After the final selection of requiring the events fall within the $`m_{EC}`$ \- $`\mathrm{\Delta }E`$ signal-box, the efficiency is $`5.2\pm 0.1(MCstat)\pm 0.5(sys)\%`$. This efficency has a number of systematic uncertainties associated with it including those arising from:

1. the trigger efficiency
2. the tracking reconstruction efficiency
3. the neutral cluster reconstruction efficiency
4. the background filter and skim selection efficiency
5. electron opposite hemisphere requirements (electron tag)
6. $`m_{opp}`$ requirements (‘$`\rho `$’ tag)
7. $`\mu `$-PID requirements
8. the photon energy scale and resolution
9. the photon direction reconstruction, scale and resolution
10. the track momentum scale and resolution
11. the track momentum direction scale and resolution
12. the beam energy scale
13. the beam energy spread

Evaluation of the efficiency done using the events in the Grand Side-Band, as these have characteristics which are very similar to those of the signal. The good agreement between the data and Monte Carlo normalization for the background sources, evident in Table 1, and the shapes of, for example the $`\mu `$ momentum, indicates that the efficiencies are well understood. The effects of systematic items (1)-(7) are incorporated into the ratio of observed to expected events: $`1.022\pm 0.069(stat)\pm 0.025(norm)`$. This yields an estimate of 7.3% on the systematic error associated with items (1)-(7).

A global estimate of tracking and calorimetry errors is provided by shifting the signal-box position and by modifying the size of the box according to the uncertainties in the resolution. Half of the observed changes are used to estimate the systematic errors. The systematic errors assessed in this manner are presented in Table 1.

Monte Carlo simulations revealed that the selection efficiency is insensitive to the uncertainties in the beam energy spread and the associated systematic error is negligible.

## 6 Limit

From the observed number of events a 90% upper limit is set on B($`\tau \mu \gamma `$) with the systematic errors included as suggested in  using the technique of .

A demonstration with the two sidebands can be made: the lower sideband had a predicted background rate of 6.1 $`\pm `$ 2.2 events, whereas 6 are observed.

When the signal-box was unblinded, we observe 13 events, shown in Figure 5 for the $`m_{EC}`$ projection. Of these, 9 had electrons identifed in the opposite hemisphere; 5 were $`h\pi ^o`$-like and one was identified as both an electron and a $`h\pi ^o`$. These tagging ratios are consistent with the expectations of $`\tau `$ decays. The distributions of the $`m_{EC}`$ and $`\mathrm{\Delta }E`$ are shown in Figure 6.

The background was estimated to be $`7.8\pm 1.4`$. This yields a limit of 11.5 event upper limit on the number of signal events @90% CL when the systematic errors are included. The 11.5 event upper limit translates into a limit: B($`\tau \mu \gamma )<2.0\times 10^6\mathrm{@}90\%`$CL. The probability of a background of $`7.8\pm 1.4`$ events up to 13 observed events in the absence of a signal is 7.6% if one includes the systematic errors.

## 7 Conclusion

The 1999-2001 BABAR data has been studied in a search for the forbidden decay $`\tau \mu \gamma `$. These studies reveal that a search with an efficiency of $`5.2\pm 0.1(MCstat)\pm 0.5(sys)\%`$ and an expected background rate of $`7.8\pm 1.4`$ results in 13 events being observed. This leads to a preliminary limit of B$`(\tau \mu \gamma )<2.0\pm 10^6\mathrm{@}90\%`$CL.


text:
# Random Vibrational Networks and Renormalization Group

\[

## Abstract

We consider the properties of vibrational dynamics on random networks, with random masses and spring constants. The localization properties of the eigenstates contrast greatly with the Laplacian case on these networks. We introduce several real-space renormalization techniques which can be used to describe this dynamics on general networks, drawing on strong disorder techniques developed for regular lattices. The renormalization group is capable of elucidating the localization properties, and provides, even for specific network instances, a fast approximation technique for determining the spectra which compares well with exact results.

\] The study of linear equations with random coefficients has a long history, dating back to the study of random matrices, systems with no spatial structure at all. Later work incorporated spatial structure, focusing on the case of regular lattices, starting with Anderson’s work studying electron localization via the Schroedinger equation. The Bethe lattice was also introduced as a soluble case, illustrating the localization transition in these systems.

Recently, the study of dynamics on more general networks has become of interest. The systems above can be included in the case of networks, as can additional systems, such as sparse random matrices, or small world networks. Networks contrast with regular lattices, as the concept of locality on networks may not be well defined. This concept of locality is essential to physics, where so many subjects rely on introducing slowly varying “hydrodynamic” fields to describe the long wavelength dynamics. Small world networks illustrate the problem: there is always a route between any two nodes in the network which involves traversing sites which are near each other in physical distance, but there may be a much shorter route which involves using some long jumps.

In this work, we consider the specific problem of a random networks of masses connected by springs. This problem is chosen for several reasons. It includes, but is richer than, the related problem of a Laplacian on a network. It also has physical applications to vibrational modes in amorphous and granular media. Finally, this problem is a good problem for studying the application of real-space renormalization group techniques, originally developed for regular lattices, to arbitrary networks, including random graphs and small-world networks. Real-space techniques were originally developed for application to statistical mechanics systems without randomness. Since then, they have been applied to a number of disordered systems, in some cases yielding exact results. The advantage of these techniques is that, without assuming any specific structure of the system as required for momentum space techniques, they are able to reveal the geometry of the network. We will see that not only are these techniques able to provide very accurate approximations, but also to find the correct notion of locality for the given network.

Random Networks, Random Masses— We consider a random network of masses connected by springs. This generalizes the problem studied by Dyson in one-dimension. The equation of motion for a particle $`i`$ of mass $`m_i`$ connected by springs with spring constants $`k_{ij}`$ to other particles $`j`$ is

$$_t^2x_i=\underset{j}{}O_{ij}x_j,$$
(1)
where the matrix $`O`$ has elements $`O_{ij}=k_{ij}/m_i`$ for $`ij`$ and $`O_{ii}=\underset{j}{}k_{ij}/m_i`$ (we define $`O`$ with these signs so its eigenvalues will be positive). By a rescaling of coordinates, $`x_im_i^{1/2}x_i`$, we can study instead the matrix $`L`$ defined by $`L_{ij}=m_i^{1/2}O_{ij}m_j^{1/2}`$. Then, $`L`$ has the same eigenvalues as $`O`$, but is a symmetric matrix.

It is possible to relate this problem to a random hopping problem, studied in quantum mechanics. Introduce a matrix $`H_{(jk),i}`$, where the index $`i`$ labels a site and the index $`(jk)`$ labels a link between sites (by definition, the links $`(jk)`$ and $`(kj)`$ are the same). Then, take $`H_{i,(ij)}=H_{(ij),i)}=\sqrt{k_{ij}/m_i}`$ and $`H_{j,(ij)}=H_{(ij),j)}=\sqrt{k_{ij}/m_j}`$. The matrix $`H`$ connects sites to links and vice-versa. The matrix $`H^2`$ is block-diagonal: it connects sites to sites and links to links. Restricting to just the block of $`H^2`$ which connects sites to sites, we find that $`L=H^2`$.

Automatically, $`O`$ has a zero mode associated with a uniform motion of all the masses. This will be essential to the renormalization procedure defined below, as we will study the behavior of $`O`$ for small eigenvalues. In the case that all $`m_i=1`$, $`L`$ is a Laplacian on the network with different couplings $`k_{ij}`$ between nodes. If instead $`m_i=1/n_i`$, with $`n_i`$ the number of nodes connected to $`i`$, and $`k_{ij}=0,1`$ depending on whether nodes $`i,j`$ are connected, then we obtain another definition of the network Laplacian. We will find that including random masses leads to, in many cases, greatly different localization properties than without; further, we will see that the renormalization procedures naturally lead to variations in the mass.

Renormalization— Inspired by real-space techniques applied to strong disorder systems in one dimension, we consider a variety of renormalization procedures. The previous work was based on a matrix such as $`H`$. Although it is not essential to the discussion below, we recall that the technique for operator $`H`$ goes as follows: pick the largest term in $`H`$. Let this term have magnitude $`\mathrm{\Lambda }`$. This term connects a site $`j`$ and a link $`(jk)`$. Let site $`j`$ also be connected to link $`(ij)`$, while link $`(jk)`$ is also connected to site $`k`$. The approximation behind the renormalization procedure consists in assuming that $`\mathrm{\Lambda }`$ is much larger than all other terms in $`H`$. Then, $`H`$ has two eigenvalues $`\pm \mathrm{\Lambda }`$. We then remove site $`j`$ and link $`(jk)`$ from the system and use second order perturbation theory to connect the link $`(ij)`$ to $`k`$. We find that that the resulting term in the renormalized $`H`$ is $`H_{(ij),k}=H_{(ij),j}H_{(jk),k)}/H_{j,(jk)}`$.

Since $`H^2=L`$, the renormalization of $`H`$ will enable us to renormalize $`L`$. Based on this fact, and the above procedure, we propose the following technique that can be applied directly to $`L`$: (1) choose the site $`i`$ with the largest $`L_{ii}`$. (2) Make the approximation that this $`L_{ii}`$ is much larger than $`|L_{ij}|`$ for any other $`j`$, and thus declare that $`L`$ has an eigenvalue equal to $`L_{ii}`$ with eigenvector concentrated on site $`i`$. (3) Remove the site $`i`$ from the network and define a new matrix $`\stackrel{~}{L}`$ connecting the remaining sites with $`\stackrel{~}{L}_{jk}=L_{jk}+L_{ji}L_{ik}/L_{ii}`$. We find that if $`i`$ has only the two neighbors, $`j`$ and $`k`$, then the masses of sites $`j,k`$ are unchanged by this procedure, while the sites are connected by a spring constant $`k_{jk}=k_{ji}k_{ik}/(k_{ji}+k_{jk})`$.

This procedure works if $`L_{ii}`$ is indeed much larger in than $`|L_{ij}|`$. In one dimension, this requires that the mass of site $`i`$ be much smaller than that of either of its neighbors. If the elements of $`H`$ are chosen randomly from a distribution with finite width, this leads to a random distribution of mass ratios of nearest neighbors in the original lattice. Then, since $`m_k/m_j`$ and $`m_j/m_i`$ are in this case uncorrelated, the distribution of the mass ratio $`m_k/m_i`$ is broader. Thus, the distribution of mass ratios broadens under renormalization, justifying the procedure.

We, however, consider a system in which the masses themselves, not the mass ratios, are chosen from a given distribution so that the mass ratios remain narrow under the renormalization and thus this procedure does not work for a one dimensional system. However, if the connectivity $`n_i`$ of site $`i`$ is large, then $`L_{ii}/|L_{ij}|n_i>>1`$. Thus for networks, this procedure can work. A one dimensional lattice remains a one dimensional lattice under this procedure, while any other lattice or network changes its topology.

Define the Green’s function $`G_L(E)=1/(EL)`$. Then, for $`E=0`$, $`G_L=G_{\stackrel{~}{L}}`$ for the sites that remain in the network, so this procedure is exact for $`E=0`$. For other $`E`$, it is possible to follow a renormalization procedure with $`\stackrel{~}{L}_{jk}=L_{jk}+L_{ji}L_{ik}/(L_{ii}E)`$.

A more powerful renormalization technique is to consider pairs of sites. This two-site technique is more accurate at each renormalization step. It also tends to increase the connectivity of the sites and randomize the masses, thus leading to a situation in which each renormalization step is more accurate than the previous. We proceed as follows: (1) choose the site $`i`$ with the largest $`L_{ii}`$, and then find its neighbor $`j`$ with the largest $`|L_{ij}|`$. (2) Consider the two-by-two submatrix of $`L`$ which involves only sites $`i,j`$. This has eigenvalues $`E^\pm =(L_{ii}+L_{jj})/2\pm \sqrt{(L_{ii}L_{jj})/2+L_{ij}^2}`$ with $`E^+>E^{}>0`$, and has eigenvectors $`v_{ij}^\pm `$. Change basis from $`i,j`$ to the basis of these eigenvectors $`v^+,v^{}`$ to diagonalize this matrix, using $`L_{kv^\pm }=L_{ki}v_i^\pm +L_{kj}v_j^\pm `$. Make the approximation that $`L_{ii},L_{ij},L_{jj}`$ are much bigger than any other term in $`L`$, and thus declare that $`L`$ has one eigenvalue equal to $`E^+`$, with corresponding eigenvector $`v^+`$. (3) Remove $`v^+`$ from the network and define a new matrix $`\stackrel{~}{L}`$ connecting the remaining sites with $`\stackrel{~}{L}_{kl}=L_{kl}+L_{kv^+}L_{v^+l}/E^+`$.

This procedure is the same as that above, but with a change of basis from $`i,j`$ to $`v^+,v^{}`$. We will see below that this makes the procedure much more accurate. In addition to the increased accuracy at each step, the procedure tends to increase the connectivity: applying this procedure to a one dimensional lattice produces second-neighbor connections. The procedure also leads to variation in the masses: the mass of the new site $`v^{}`$ is not equal to either $`m_i`$ or $`m_j`$. Thus, the procedure drives itself to a regime in which it become more accurate, due to higher connectivity and a broader distribution of masses. This procedure shows the correct idea of locality for the network: we know that, as far as the low-energy dynamics are concerned, $`i,j`$ are “close by” in the original network.

Localization Properties— We have performed numerical simulations to study localization. The specific network considered is a random graph of $`N`$ sites: for a given probability $`p`$, two sites are connected with that probability. Then, the spring constant connecting those sites is assigned from a uniform distribution between $`0`$ and $`1`$ (taking instead all $`k_{ij}=0,1`$ leads to little change in the results). The average connectivity of a site is equal to $`p(N1)`$. We consider three different mass distributions: (1) $`m_i=1`$ for all $`i`$; (2) $`1/\sqrt{m_i}`$ is chosen uniformly between $`1`$ and $`2`$, giving a smooth, bounded distribution of masses; (3) $`1/\sqrt{m_i}`$ is chosen uniformly between $`0`$ and $`1`$, so that the distribution of masses is unbounded.

The energy scales with the connectivity of the the system. After removing this scaling, we find that for a system with extensive connectivity and a bounded mass distribution, the spectrum has a single zero mode and then a gap to the next eigenvalue $`E`$. For a system with intensive connectivity or with unbounded masses, the gap is filled in by Griffiths effects.

To determine the localization properties of a normalized eigenfunction $`\psi _{i,l}`$ associated with an eigenvalue $`E_l`$, we consider the inverse participation ratio, $`_i|\psi _{i,l}|^4`$. In the delocalized phase, this quantity scales with $`1/N`$. For constant mass and large connectivity, in the center of the band the states are delocalized, while they are localized near the band edge.

We have found a very different result in the case of random masses. With large connectivity and random masses, we find that all the states, with the exception of the zero mode, become localized as shown in Fig. 1. We plot the inverse participation ratio against eigenvalue, for systems with bounded, random mass. The solid line is $`N=1000,p=1`$ averaged over $`100`$ samples while the dashed is $`N=500,p=1`$ averaged over $`1000`$ samples. In the inset, we have divided the eigenvalue by $`N`$, showing a perfect collapse of the two curves, indicating that the states are localized.

For extensive connectivity, localization holds even for very weak variation in the masses. This can be understood analytically via the renormalization. Let us divide $`L`$ by $`N`$; then the $`L_{ii}`$ are all fluctuating variables of order unity, with fluctuations, due to the random mass, which are also of order unity. The $`L_{ij}`$ are of order $`1/N`$. Using the RG to remove sites from the system at eigenvalue $`E0`$, removing a single site $`j`$ leads to corrections of order $`L_{ij}^2/(L_{jj}E)`$. The numerator is of order $`1/N^2`$, while due to the random mass the denominator does not vanish. Thus, the dynamics of a single site is only weakly affected by other sites, and the states are localized.

For low connectivity, we have considered systems with 100, 200, 400, 800, and 1600 sites, averaging over 1000 realizations for systems with up to 800 sites, and over 100 realizations for the largest systems, considering the three different mass distributions, with average connectivities of 1,2, and 4. One expects that for low connectivity, greater than 1, the network approximates an infinite dimensional system, or Bethe lattice. Then, at zero eigenvalue there is an extended state, and for low eigenvalues the variations in masses and spring constants only weakly scatter the vibrational waves. Thus, one expects a delocalized phase for small eigenvalues for average connectivity greater than $`1`$.

Numerically, the results for low connectivity are less clear. For average connectivity equal to $`1`$, the inverse participation ratio is found to be independent of $`N`$, and thus the system appears localized. For average connectivity greater than $`1`$, the inverse participation ratio decreases with increasing $`N`$; however, we do not observe a clear $`1/N`$ scaling of the inverse participation ratio for an average connectivity of $`2`$. For average connectivity $`4`$, we do see approximate $`1/N`$ scaling for low $`E`$, indicating delocalized states.

Comparison of Eigenspectra— We have then tested the renormalization procedure against exact results by applying the procedure to single realizations of the system. We have found that in the correct regimes, the procedures are highly accurate even for the details of specific realizations. The procedure work best for high connectivity, and for wide distributions of the masses, with the advantage of the two-site technique being that it drives the system to the regime in which it works accurately.

In Fig. 2a, we show the results of the single-site technique for a system with 100 sites, $`p=1`$, and the unbounded distribution of masses. The sorted eigenvalues are plotted, with the $`n`$-th eigenvalue $`E_n`$ plotted at position $`(n,E_n)`$. The dashed line gives the results of the renormalization procedure, while the solid line is the exact result. For low eigenvalues, the two are indistinguishable. In Fig. 2b, we show the same for a system with 1000 sites; the two lines cannot be distinguished on the figure.

At lower connectivity, with less randomness in the mass, the two-site procedure becomes necessary. In Fig. 3a, we show a system with $`N=100,p=.01,`$ and the bounded, random distribution of masses. The solid line is the exact result, the upper dashed line is the two-site result, and the lower is the single-site result.

As one measure of the accuracy of the procedure, we sort the eigenvalues from highest to lowest, and divide the root-mean-squared error in the $`n`$-th eigenvalue by the root-mean-squared value of that eigenvalue. The result is quite accurate, with the relative error averaging, for example, $`1\%`$ for $`N=100,p=1`$ and bounded, random mass, becoming even more accurate for larger systems. As a more stringent test of the accuracy, if we instead divide the root-mean-square error by the root-mean-square sample-to-sample fluctuations in that eigenvalue, the relative error averaged over all $`n`$ is $`15\%`$ for the same system. For lower connectivity, the relative error (compared to sample-to-sample fluctuations) is worse, $`28\%`$ for $`N=100,p=.01`$. However, for the smaller eigenvalues, the accuracy increases: for $`N=100,p=.01`$, the error relative to fluctuations in the 80 lowest eigenvalues averages $`11\%`$.

Finally, we test the renormalization on a periodic one dimensional chain with constant masses and springs constants, combining both low connectivity and no randomness. If, in the two-site procedure, when faced with a choice between different sites with the same $`L_{ii}`$, we choose randomly, the low energy properties are reasonable, while if we deliberately make poor choices, keeping the connectivity low, the results are much worse, as shown in Fig. 3b. Using a network improves on the results obtained on a one dimensional chain with a similar renormalization procedure, which preserves the one dimensional structure. In the inset to Fig. 3b we consider a small-world network. We take a periodic, one dimensional chain with $`500`$ sites and constant mass, and with probability $`p=0.0002`$ we connect pairs of sites which are not already connected by the chain.

Conclusion— In conclusion, we have introduced a real-space renormalization procedure for studying certain linear operators on networks. The procedure is exact at low energy, and in many cases leads to very good results for the spectra even at much higher energy. This procedure can be used to obtain analytic results on the localization properties. It can also be used to obtain results for dynamics on a specific graph more rapidly than would be possible with matrix diagonalization routines. While for a fully connected graph it still requires $`O(N^3)`$ time, the prefactor is much smaller than that required for matrix diagonalization. For less connected networks, the procedure runs much more rapidly. This difference will be important when studying large social networks which may have billions of nodes.

We have also found that the introduction of random masses leads to many more localized states than constant masses. For amorphous systems, the random mass case is the relevant case.

The procedure naturally helps simplify the network by finding a simpler network with a similar spectrum. Although we have considered only linear problems, we hope that this kind of real-space procedure will be an important tool in studying nonlinear problems in future work, such as spin-glass or optimization problems.

Acknowledgments— I thank Z. Toroczkai for inspiration and discussions, and C. Jarzynski for suggesting the relative error measurement. This work was supported by DOE W-7405-ENG-36.