format: text

text:
# NIKHEF-00-032 New features of FORM

## 1 Introduction

It is already more than 10 years ago that the symbolic manipulation program FORM was introduced. At the time it was constructed more or less as a successor to SCHOONSCHIP  which, due to the limitations of its language of implementation, became harder and harder to use. In addition the continuous development of methods of computation in many fields of science required a rapid expansion of the possibilities of the symbolic programs. Admittedly programs like Reduce, Maple and Mathematica offer a host of possibilities combined with extensive libraries, but the speed of SCHOONSCHIP was always something that set it apart from the other programs. Hence the task was to create a program that would be at least as fast as SCHOONSCHIP and be fully portable. In particular it should be usable on the fastest computers that are suitable for this kind of work (eg. vector processors will rarely give great benefits to symbolic manipulation in general. As we will see this is completely different for parallel processing). This resulted in version 1 of FORM. Due to more modern algorithms and the benefits of a second generation approach it is actually somewhat faster than SCHOONSCHIP when the problems remain limited in size and much faster when the expressions become large. Compared to the more popular computer algebra systems the difference in speed and use of memory depends largely on the type of the problem. This is due to a completely different philosopy of design and use. Typical factors are between 10 and 100, but there have been rare cases for which the factor was considerably less and cases for which the factor was considerably larger. Such a factor translates into a timeshift with respect to when a given calculation can be done in practise with other programs. Over the past years the computer industry has given us on average every three years a factor 4 in increase in capabilities. Hence for FORM this translates into a shift of more than 5 years. As a consequence most of the calculations in quantum field theory for which the available computer resources are a consideration have been using FORM over the past years.

Of course the development of FORM has continued. Version 2.0 was introduced in 1991. It contained a number of improvements that were inspired by the use of the program in a number of actual research projects. Since then regular small improvements were made, dictated by more and more complicated projects. After a number of years it became clear that a number of fundamental internal changes would enhance the capabilities very much and hence work on version 3 was started. Most of the code of version 2 was replaced, introducing better structured ways of dealing with internal objects. Even though the language of implementation is still regular C, the methods used are in many cases those of object oriented programming. It was of course important to do this in such a way that the speed and compactness of the expressions would not suffer from this. This has been achieved. After this revamping the new facilities were added. Some of them are completely new in the field of symbolic manipulation. This does not set the claim that such things are impossible with other programs. It just means that the FORM approach to some things is far more direct and often easier to program, even if some people may not like the language.

In this paper the most important new features are described shortly. In most cases some short examples of code are given with the corresponding output. It should be understood that this paper is not a manual. Because a number of packages comes with the FORM distribution, the examples to illustrate the features come mainly from these packages. In addition there are examples from some incomplete future packages. To improve the understanding of the above, first a very short description of FORM is given.

There is also a short section on a parallel version of FORM. It looks like this will be a major new feature that is currently still under development but results have been obtained already with one of the packages. The actual target is to run nearly unmodified FORM programs (the mincer package that has been used already was not modified at all) on a number of processors in such a way that to the user it looks like a much faster sequential machine. This is very useful for –among others– the development cycle of big programs.

The final section discusses the availability of FORM. In the appendix we give a short description of the packages with some examples of their use.

## 2 A short overview of FORM

Here a very brief overview of FORM is given. For a complete introduction one can consult the manual which contains a complete tutorial and a reference section.

FORM is not an interactive system. It runs in batch. In UNIX terms one would say that FORM is a filter. For large programs this is the natural mode anyway. It also allows the user to prepare programs with a familiar editor. The editor used by the author is stedi which has a folding capability that is recognized by FORM. The sources of this editor are available in the FORM distribution as well. If everan interactive version will be constructed, it will be based on this editor.

The language of FORM uses strong typing. This means that all objects that will be used in a program will have to be declared. This allows FORM to use a number of properties of such objcets and this adds to the speed of the program. Algebraic objects in FORM are expressions, symbols, functions, vectors and indices. In addition there are $-expressions, sets, numbers, preprocessor variables and argument field wildcards. The functions are subdivided in categories of commuting and noncommuting functions. Additionally there are the classes of tensors, tables and regular functions.

The basic modus operandi of FORM is that variables are declared. Then one or more expressions are defined by giving them initial contents. After that operations are applied to one or more of these expressions.

A FORM program is subdivided in modules. Each module is an individual entity of which the contents are compiled when it is the turn of the module to be executed. After the module has been executed, the compiled version is forgotten again, but depending on the type of the module, declarations and expressions will be kept. Hence at the level of complete modules FORM acts as an interpreter. The contents of expressions as they were at the start of the module will be forgotten, unless the user takes some special action. This is needed to prevent the disk from filling up too quickly. At the beginning of each module extra variables can be declared and new expressions can be defined. If necessary old expressions can be deleted. In principle the operations of the module are applied to all old and new expressions.

The basic object for manipulation is the term. Expressions consist of terms, separated by addition operators (or subtraction of course). This addition operator is commutative and associative. Function arguments can consist of subexpressions, i.e. smaller expressions which are limited in size. Such subexpressions also consist of terms. Terms consist of subterms. The subterms are separated by the multiplication operator. This multiplication operator is associative and usually commutative. Only in the case of noncommuting functions and still unsubstituted objects that could contain noncommutative functions is the multiplication order sensitive. The number zero is not considered to be a term. An expression that is zero contains no terms.

FORM works its way through an expression in a sequential manner. It starts with the first term, applies all statements of a module to it and puts the results away for sorting (normal ordering of the complete expression is called sorting here. It does involve adding coefficients of terms that are identical with the exception of their prefactor). Partial sorting may occur when buffers become too full. After all terms of an expression have been treated this way the sorting is completed. Once all expressions have been treated this way the new versions of the expressions replace the old versions.

The actions inside a module can be either predefined operations or substitutions that may involve rather complicated patterns involving wildcarding. Wildcarding is the use of generic variables that represent a whole category of objects, both in the left hand side (the pattern) and in the right hand side (the replacement expression).

A number of special functions is available to make the use of the program easier and often much faster. These involve amoung others a sum function, a number of combinatoric functions and a number of numerical functions. In addition there are special tensors like the Kronecker delta and the Levi Civita tensor. And there is much more.

For the use in particle physics there are Dirac gamma matrices and trace operations. It should be noted however that this is nowadays a rather unimportant part of the program as there exist many facilities to emulate these objects and operations (and similar operations that can be much more complicated).

Results of a calculation can be presented in many different ways. One way is to print them to regular output. They can be written to file in internal notation for use in future FORM programs. Or they can be written to special files, etc.

There is a powerful preprocessor that can manipulate the input at the character level. It is fully adapted to the requirements of symbolic manipulation. There are now also some mechanisms for the feedback of algebraic results to the preprocessor. This facilitates the possibility to have FORM programs write parts of themselves.

A number of buffers in FORM have to be allocated at startup. The user can influence these allocations with a special setup file. Ordinarily this is considered an advanced feature.

## 3 New features

For the description of the new features it is assumed that the user has some experience with the older versions. The distribution contains a good tutorial. The interested reader without experience could look at the many examples in it. Here we concentrate on the new concepts. For a complete description and the exact syntax and for the many additional new features we refer to the new manual.

### 3.1 Variable administration

The complete internal adminstration of variables has been made dynamic. This means that the old restrictions on the number of variables have been lifted. The setup file is no longer needed when more than 100 variables of the same type are asked for. The limit lies now at 6000 for 32-bit systems and $`10^8`$ for 64-bit systems.

Similarly the compiler buffer, the input buffer and the various preprocessor buffers are arranged dynamically. The reason that in the older versions these things were set up in a static way was due to the size of the computer memories in the late 80’s.

### 3.2 Preprocessor variables

Preprocessor variables can be used ‘recursively’. The old notation to enclose a preprocessor variable in quotes has been changed into a notation in which the variable is enclosed in a backquote-quote pair as if these are brackets. This way they can be nested as in

```
#do itabs = 1,‘size’
#ifndef ‘STABLE‘itabs’HFILE’
#call table‘itabs’
#endif
#enddo

```

In the summer package there are 8 files table‘itabs’.prc all having something like:

```
#procedure table5
#ifndef ‘STABLE5HFILE’
#define  STABLE5HFILE "1"
    lots of tables
#endif
#endprocedure

```

This way the procedures are not called twice. Of course the contents would not be used twice anyway, but some of them are rather lengthy and FORM would have to read the whole procedure to find the $`\mathrm{\#}`$endif. That is avoided this way.

To make the preprocessor variables even more powerful there are also postincrement and postdecrement operators for preprocessor variables as in the C language:

```
     #define k "1"
     Local F‘k++’ = ...;
     Local F‘k++’ = ...;
     Local F‘k++’ = ...;
     Local F‘k++’ = ...;

```

would result in the definitions of `F1` to `F4`. In the end `k` would have the value 5. This is particularly handy when the program is generating expressions and we do not know in advance how many.

### 3.3 Preprocessor abbreviations

The notation has been cleaned up a bit by introducing a new operator for the compilation stages. This operator is indicated by three dots as in

```
    Symbols x1,...,x100;
    Local Fac10 = 1*...*10;

```

Actually there is an even more general notation involving $`<>`$:

```
    id  f(<p1?,m1?,i1?>,...,<p6?,m6?,i6?>) =

```

Here the $`<>`$ enclose a pattern. FORM will look for the differences in the begin pattern and the end pattern, calculate their numerical differences and if they all have the same absolute value we can make a ‘running pattern’ and make the substitution. Hence

```
    id  f(<p1,m4>,...,<p4,m1>) =

```

would evaluate into

```
    id  f(p1,m4,p2,m3,p3,m2,p4,m1) =

```

This facility cannot be used recursively but already it can clean up the code considerably when compared to how one would have to do this with do loops in the old versions of FORM. It becomes really powerful in combination with preprocessor variables:

```
    #do j = 1,‘MAX’-1
        id  e_(i1?,...,i‘j’?) =
                sum_(k,1,‘MAX’,e_(i1,...,i‘j’,k)*f(‘j’+1,k));
        id  f(‘j’,k?) = tab(‘j’,k);
        .sort:step ‘j’;
    #enddo
    id  e_(1,...,‘MAX’) = 1;

```

This would be the central step in computing the determinant of a matrix that has been defined in the two dimensional table named tab. We use an intermediate function f to make sure that first the Levi-Civita tensors with two indentical indices are killed before table substitution is attempted. This is a very efficient algorithm for sparse matrices, especially when the matrices have been ordered in such a way that the first columns have the largest number of zeroes.

### 3.4 $variables

These variables are at the same time small algebraic expressions and character strings. They can be used as either. When they are used and they are enclosed between a backquote-quote pair, the compiler will make their character representation and substitute this like a preprocessor variable. This use as a preprocessor variable even includes the possibility of using the postincrement or postdecrement operators. If they are not enclosed between the backquote-quote pair the contents are seen as an algebraic expression to be used during execution time. This can lead to very flexible code. If they are to be given an initial value, this can be done with the character $`\mathrm{\#}`$ before their name to have the preprocessor provide the value (hence this is done once during compilation time), or without the $`\mathrm{\#}`$ to indicate that this value should be provided during execution, each time a term passes this position in the code. There are additional methods to provide the $-variables with a value that is the match of a wildcard variable. Let us have a look at an example:

```
    #$max = -100;
    if ( count(x,1) > $max ) $max = count_(x,1);
    .sort
    #do i = 1,‘$max’
        ...
    #enddo

```

Here we give `$max` the initial value -100. During execution we look for the maximum power of the variable `x`. Once the module has finished this maximum should be indeed the value of `$max`, provided it is at least -100. At this point we can use this value as a parameter in a do loop. Note that the .sort is essential, because if we were to omit it, the do loop would have the bounds 1,-100 because the do loop instruction would be translated before any algebraic manipulation would have taken place (first the complete module is translated, then execution takes place).

It is also possible to manipulate the contents of the $-variables, provided these do make sense.

The contents of these $-variables are kept in memory. This means that the user should not try to see them as a substitute for regular expressions and have some of them grow to a large size. It would work, but performance might suffer when the physical memory is exhausted in the same way that Mathematica and Maple become rather disagreeable when they need more memory than is available.

### 3.5 Hide

When there are several expressions active at the same time, sometimes there is the need to temporarily put some of them ‘out of the way’. FORM has several mechanisms for this. The Skip statement would just skip the indicated expression(s) by copying them directly from the input to the output of the module. This is not always efficient. In addition one has to specify this in each module seperately. A more basic solution is to copy the expression(s) to a special location where they remain suspended until called back. This is done with the Hide and Unhide statements. The hidden expressions can still be used in the right hand side of expressions in the same way as fully active expressions can be used. This is a faster mechanism than the one that is used with stored expressions.

```
    .sort
    Hide F1,F2;
    Local F4 = F1+f3;
    .sort
         more modules
    .sort
    Unhide F1,F2;
    id  x = y+1;

```

Here F1 and F2 are not operated upon for a few modules. Then they are called back and they are immediately active again. This means that the substitution of x also holds for the terms of F1 and F2.

There is another mechanism when one or a few statements should be applied only to a limited number of expressions. This concerns a new option of the if statement:

```
    if ( expression(G) == 0 );
       statements
    endif;

```

In this case the statements are applied to the terms of all active expressions except for the terms of the expression G.

### 3.6 $`\mathrm{\#}`$write

This is a preprocessor controled instruction that can write either to regular output, to the log file or to any named file. It has a format string as in the printf function of the C language. It specifies what should be printed. There are some examples in which complete procedure files are constructed this way for later use in the same program.

### 3.7 Print ””

The new feature of the print statement is a mode like printf in the C language that can print information during execution. It can handle the printing of messages and even individual terms or $-variables. This gives it great value as a debugging tool. A good example is given with the SplitArg statement in 3.11.

### 3.8 Table manipulations

This facility involves a statement, a special function and a special print statement. The table\_ function can be used to load the contents of a table into an expression. The table indices are represented by powers of a given set of symbols. The fillexpression statement does the opposite. It takes an expression that is bracketed in a number of symbols and loads the contents into a table. The powers of the symbols indicate the table indices and the contents of the brackets form the corresponding table elements. The printtable statement writes a table to file in the form of fill statements. This file can be used in a future program to fill the table. Together these facilities allow the dynamic extension of tables in a program. If a table element that is needed has not yet been defined, one could compute it on the spot and add it to the table. The results can be written in such a way that in future programs the same element does not have to be recalculated.

### 3.9 Symmetry properties

Functions can have properties of symmetry. These properties hold for the set of all arguments. The functions or tensors can be symmetric, antisymmetric, cyclesymmetric and reversecyclesymmetric. This last means that the function is invariant both under cyclic permutations and a complete reversal of all arguments. These properties are put in the declaration of the function or tensor. Each time some arguments of the function are changed, FORM will check whether permutations have to be made. There are some restrictions. FORM will not accept argument field wildcards (wildcards that match zero, one or more function arguments) in patterns in which such a wildcard occurs inside a symmetric or antisymmetric function. The pattern matching would become extremely slow when no match occurs (all possibilities have to be searched and that involves n! possibilities for n arguments). There are tricks to circumvent this in most cases.

### 3.10 Topology

In combination with the symmetry properties FORM has now a command to look for loop structures in argument contractions. What this means is that if we have a network in which v is a vertex function and its arguments are indices and connections are represented by two vertices having identical indices, FORM can find the smallest loop in such a structure and make replacements based on this. There are some examples in the Color package in which the vertex function is the antisymmetric tensor f with three indices. In the program

```
#-
Indices i1,...,i9;
CF f(antisymmetric),ff;
Local F = f(i1,i2,i3)*f(i2,i4,i5)*f(i3,i5,i6)*
          f(i4,i7,i8)*f(i6,i7,i9)*f(i1,i8,i9);
ReplaceLoop,f,arguments=3,loopsize=all,outfun=ff;
Print +f;
.end

```

the replaceloop command takes the smallest loop, removes the loop indices and puts the remainder in the output function ff which should be at least cyclesymmetric (at least the action we take here implies it). If there is more than one possibility FORM takes the first one it encounters. The result of this program is

```
   F = f(i1,i8,i9)*f(i4,i7,i8)*f(i6,i7,i9)*ff(i1,i6,i4);

```

Only a single loop has been replaced. If one would like to replace all loops one should put a repeat/endrepeat combination around the replaceloop statement. This would result in

```
   F = - ff(i1,i4,i6)*ff(i1,i6,i4);

```

If such a command would not exist, the color package would be much more restricted because it would need a statement for each loop (loops with one vertex, two vertices, three vertices etc.) Also the corresponding pattern matching would be much slower. In the crucial phases of the program this made a factor 10 difference in speed and it also made the code much shorter.

### 3.11 SplitArg

This is a rather useful tool for investigating arguments of functions. In its basic form it splits a multiterm argument into a number of arguments, each with a single term. There are variations for faster operation when only a single term should be taken from the multiterm argument. Take for instance the following example in which we want to select all functions den that contain the symbol `j1` in its argument.

```
    SplitArg,((j1)),den; * take all terms in j1 separately
    id  den(j1) = den1(0,j1); * a single term does not get split
    *                           we may need more of these statements
    id  den(x?,x1?) = den1(x/x1*j1,j1)/x1*j1; * normalize the j1 argument.
    *
    * all den with two arguments have now just j1 without factor.
    *
    repeat id den1(x?,j1)*j1 = 1-x*den1(x,j1);
    repeat id den1(x1?!{x2?},j1)*den1(x2?!{x1?},j1) =
               (den1(x1,j1)-den1(x2,j1))*den(x2-x1);
    id  den(x?number_) = 1/x;

```

Code of this nature occurs frequently in the summer package. It splits the fractions, being careful not to divide by zero (in some cases one has to be far more careful and involve theta\_ and delta\_ functions). Let us see what this code does. We use the following program (this is also a demonstration of the new print statement capabilities):

```
#-
S   j1,x,x1,x2;
CF  den,den1;
Off Statistics;
L   F = den(j1)*den(2+j1)*den(3-2*j1);
  Print +f "<1> %t";
SplitArg,((j1)),den;
  Print +f "<2> %t";
id  den(j1) = den1(0,j1);
  Print +f "<3> %t";
id  den(x?,x1?) = den1(x/x1*j1,j1)/x1*j1;
  Print +f "<4> %t";
repeat id den1(x?,j1)*j1 = 1-x*den1(x,j1);
  Print +f "<5> %t";
repeat;
    id den1(x1?!{x2?},j1)*den1(x2?!{x1?},j1) =
         (den1(x1,j1)-den1(x2,j1))*den(x2-x1);
    Print +f "<6> %t";
endrepeat;
id  den(x?number_) = 1/x;
  Print +f "<7> %t";
Print +f;
.end

```

which results in

```
    #-
<1>  + den(j1)*den(2 + j1)*den(3 - 2*j1)
<2>  + den(2,j1)*den(3, - 2*j1)*den(j1)
<3>  + den(2,j1)*den(3, - 2*j1)*den1(0,j1)
<4>  - 1/2*den1(0,j1)*den1(2,j1)*den1( - 3/2,j1)
<5>  - 1/2*den1(0,j1)*den1(2,j1)*den1( - 3/2,j1)
<6>  - 1/2*den(2)*den1(0,j1)*den1( - 3/2,j1)
<6>  - 1/2*den(2)*den( - 3/2)*den1(0,j1)
<6>  - 1/2*den(2)*den( - 3/2)*den1(0,j1)
<7>  + 1/6*den1(0,j1)
<6>  + 1/2*den(2)*den( - 3/2)*den1( - 3/2,j1)
<6>  + 1/2*den(2)*den( - 3/2)*den1( - 3/2,j1)
<7>  - 1/6*den1( - 3/2,j1)
<6>  + 1/2*den(2)*den1(2,j1)*den1( - 3/2,j1)
<6>  + 1/2*den(2)*den( - 7/2)*den1(2,j1)
<6>  + 1/2*den(2)*den( - 7/2)*den1(2,j1)
<7>  - 1/14*den1(2,j1)
<6>  - 1/2*den(2)*den( - 7/2)*den1( - 3/2,j1)
<6>  - 1/2*den(2)*den( - 7/2)*den1( - 3/2,j1)
<7>  + 1/14*den1( - 3/2,j1)

   F =
       - 2/21*den1( - 3/2,j1) + 1/6*den1(0,j1) - 1/14*den1(2,j1);

```

Note that the termwise option of the print statement gives us a trace of what happens and in what order.

### 3.12 The term environment

Sometimes one needs to work things out for just a few terms, but an intermediate sort for the results of these manipulations is needed in order to have things cancel. Or one would like to do even fancier things like trying whether the contents of certain brackets can be factorized by a given expression. If so, one would like to make that factorization. The term environment defines the current term temporarily as a small expression. Because this is done during the processing of the terms, sorting cannot be done with the .sort instruction which works for the whole expression. This means that we need a separate sort statement for all these small expressions. It should however be noted that the complete term environment consisting of the term statement, the endterm statement, and the intermediate regular statements and sort statements are all part of the same module and hence are all compiled before execution starts. Basically we use the same mechanism of operation in the argument environment. There function arguments can be operated upon and when the endargument statement is encountered the results are sorted and put back inside the function. The sort statement is hence equivalent to an endargument statement followed by a new argument statement. This last construction has the limitation that the result has ‘to fit’ and hence is limited to MaxTermSize or less. The Term environment does not suffer from this restriction because we do not have to keep the results inside a function argument and also in the intermediate stages even the disk can be used and the expression may be quite large.

```
#-
S   x,y,[x+1],[x+2],x1,x2;
CF  acc;
L   F = y*(x+1)^2*(x+2)+y^2*(x+1)*x+y^3*(x+2)*x+y^4*x^2;
B   y;
Print +f;
.sort
Collect,acc;
Term;
    $min1 = 1000;
    $min2 = 1000;
    id acc(x?) = x;
    id  x = x1-1;
    sort;
    if ( count(x1,1) < $min1 ) $min1 = count_(x1,1);
    sort;
    Multiply ([x+1]/x1)^$min1;
    id  x1 = x+1;
    id  x = x2-2;
    sort;
    if ( count(x2,1) < $min2 ) $min2 = count_(x2,1);
    sort;
    Multiply ([x+2]/x2)^$min2;
    id  x2 = x+2;
EndTerm;
Print +f;
.end

```

which gives

```
    #-

   F =
       + y * ( 2 + 5*x + 4*x^2 + x^3 )
       + y^2 * ( x + x^2 )
       + y^3 * ( 2*x + x^2 )
       + y^4 * ( x^2 );

   F =
      x*y^2*[x+1] + x*y^3*[x+2] + x^2*y^4 + y*[x+1]^2*[x+2];

```

Of course this works only if one has a suspicion which factors are involved and things would go much better if FORM would be able to deal with factorizations by itself. This is considered for the future. For now this provides at least a partial solution.

### 3.13 Bracket indexing

In the old versions FORM would look for the contents of brackets inside an expression in a sequential manner. For relatively small expressions in computers with not much memory this was still not optimal, but also not a complete disaster. With the growing of the computers and hence the expressions the inefficiency of this method became a major factor as it is an algorithm that has a linear dependence of the size of the expression. Hence now the brackets can be indexed if the user asks for it. The index will of course take space of its own, and this space can approach (in principle) the size of the expression if the brackets have very little contents. This could result in bad problems again. Therefore there is a maximum to the size of the index (can be changed in the setup file) and when this maximum is reached FORM will start skipping brackets in the index. If such a skipped bracket is asked for, FORM will do a linear search, starting from the nearest bracket in front of it. The result is still a speedup by a big factor. The brackets in the index are located by a tree-like binary search which gives a logarithmic dependence of the number of brackets in the index. The degrading factor when brackets are skipped is proportional to the size of the expression divided by the number of brackets in the index. Hence in that case the improvement is proportional to a factor $`B/\mathrm{ln}(B)`$ with B the number of brackets in the index. Usually however there should be no problem in fitting the entire index inside the memory. Because this use of brackets is not very common and the making of the index does cost resources, the default is to not make an index. Here we have a simple example

```
#-
Symbols  x1,...,x10;
Local    F = (x1+...+x10)^10;
Bracket+ x1;
.sort
Drop F;
#do i = 0,10
Local    F‘i’ = F[x1^‘i’];
#enddo
.end

```

When we run this program the second module takes 1.59 sec on a pentium 300. Without the + in Bracket+ this step takes 2.84 seconds and this is not a very big expression yet. The program that solves the equations for the values in $`N\mathrm{}`$ in the summer package had an improvement factor of more than 100.

### 3.14 Error messages

The error messages have been changed completely. The most noticeable change is the indication of where the error occurs when the input listing is off. This mode is standard for procedures. There are additional warnings and there is a way to let FORM print its interpretation of the input. Together with the already existing printing of the name lists this provides extra tools for the debugging of FORM programs. In this example we have made two errors. It is the program of the bracket index. The name of the file is ex1.frm. Because we had turned off the listing, the line number information is rather crucial:

```
#-
Symbols  x1,...,x10;
Local    F = (x1+...+x10)^^10;
Bracket+ x1;
.sort
Drop F;
#do i = 0,10
Local    F‘i’ = FF[x1^‘i’];
#enddo
.end

```

When this program is run we obtain:

```
    #-
ex1.frm Line 3 --> Illegal position for operator: ^10
ex1.frm Line 8 --> Undeclared variable FF
++++Errors in Loop

```

Because nearly all error messages have been changed and actually the complete compiler has been replaced, it is possible that some error messages either do not make complete sense or are too vague. The author is open for suggestions.

## 4 The parallel version

Because of the way FORM processes expressions it should be possible to achieve a high degree of parallelization for the complete program. After an initial study at the special computer at FNAL in 1991 this assumption was verified. It took however till the late 90’s before a more serious attempt was made. This was made possible in the framework of a large phenomenology project centered at the university of Karlsruhe (DFG Contract Ku 502/8-1). Two people worked part time on this (Albert Rétey and Denny Fliegner) and after struggling for a long time with the hardware and the parallelization library software they did indeed manage to obtain a working version of FORM that can run the Mincer package on an eight processor machine with a speed that ranges from 3 to 5 times the speed of a single processor. In addition many tests were made with other architectures .

It should be noted that thus far the parallelization has been rather crude. The terms of the input expression in a module are distributed over the various processors. These do their work, and they sort their results. After this the outputs are collected by the master processor for a final merge and the output is written to disk. It is possible to be far more refined about this and make a much better load balancing by interfering, if needed, in the main expansion tree. The main problem however is the speed of transfer of operation. Thus far use has been made of the MPI protocol (and in tests of much more of course). In the future it is expected that the use of threads will increase the data transfer rate considerably. It is of course this transfer rate that puts a limit to the eventual speed that can be obtained. A given transfer speed in combination with slow processors will lead to saturation for a much larger number of processors than in the case that fast processors are used. In a shared memory machine with threads however the transfer rate would be mostly a function of the processor speed and hence much is expected of this model. Unfortunately the implementation is not yet complete.

At the moment the only readily available version is the one for a Compaq Alpha server (with 8 Alpha 21264 processors). It is expected that more will follow. A time prognosis cannot be given yet.

Of course it is not always possible to run code in parallel. Especially the new $-variables could cause problems because their eventual contents depend on the order in which terms passed. If each processor keeps its own version of a $-variable the processors may not have identical values at the end of the module. Take the simple code

```
#$max = -100;
if ( count(x,1) > $max ) $max = count_(x,1);

```

which tries to find the maximum power of x. The various processors may have a different value in the end. Hence FORM should not parallelize modules in which $-variables obtain values during running time, unless the user can provide some extra information. This can be done in the moduleoption statement. In that statement one can tell FORM that the final value of the $-variable is either a sum, a maximum, a minimum or completely unimportant. In such cases FORM knows what to do with the results and the module can be run in parallel after all. It should be realized that this is all very new. Hence it could well be that many new facilities will be developed here and that there may be more ways in which future versions of FORM can decide about parallelization.

Any statement that concerns the parallelization will be completely ignored by a sequential version of FORM. This means that when a program is being prepared one could or should include already some of the indications for the parallel version. This way the upgrading to a parallel machine would be rather smooth.

## 5 Availability

Whereas version 2 of FORM has been commercial, it was judged that the scientific community as well as the reputation of the author would benefit more from a free distribution. Hence, starting version 3, FORM will be freely available again. For the moment this will involve a number of binary executables of the program. Amoung these will be at least executables for LINUX on PC architectures and for alpha processors running UNIX. More binaries will become available in the near future. The manual will be available both as a postscript file and as a .pdf file. This will allow the average user to obtain the system and print out its manual as well as view it on the computerscreen.

Of course it is hoped that this development will not only improve the visibility and popularity of FORM, but will also lead to more possibilities to improve FORM itself. The potential for this depends of course on people in charge of jobs knowing about the importance of FORM in the science community. Hence it is expected that people who use FORM for scientific publications will refer to the current paper. This way the use of FORM can be ‘measured’. The author hopes that the free distribution in combination with a proven popularity will lead to more people being involved in the development of FORM itself. There is still a whole list of potential improvements waiting for implementation. Undoubtedly many users will have their own suggestions. Also more packages would be appreciated.

The FORM distribution contains:

* Executables for a variety of computers.
* The conv2to3 program with its sources. This program converts old FORM code into new FORM code.
* A postscript and a pdf version of the manual and a tutorial by André Heck.
* A number of packages. Currently these are ‘summer’, ‘harmpol’, ‘meltran’, ‘color’ and ‘mincer’.
* The sources of the stedi editor. This editor works on BSD 4.0 systems and systems that still have an old compatibility mode. It also works under LINUX. There is a manual with it. Originally it was an editor for the Atari ST computers, made to work as well under MS-DOS. Then it was ported to BSD UNIX where it did run inside terminal windows. Under LINUX and X-windows it has recovered most of its old ease of use.
* The sources of the minos database program for organizing large calculations. There is a small manual with it. The program should work on nearly any system. The database files should be system independent.
* The axodraw system for including simple drawings and figures inside files. This includes the manual.

The files can be obtained from the www pages in http://www.nikhef.nl/$``$form.

Disclaimer: The files are provided without guarantee whatsoever. If some commands do not work as mentioned in the manual or as the user thinks they should work, the user should submit a bug report. This may or may not be reacted to depending on circumstances.

Acknowledgements: The author is thankful for many comments during the development of FORM. Many people contributed in the form of bug reports and useful suggestions. Some people however stand out. Denny Fliegner and Albert Rétey for developing the parallel version of FORM (with support of DFG Contract Ku 502/8-1 (DFG-Forschergruppe ‘Quantenfeldtheorie, Computeralgebra und Monte Carlo Simulationen’), Geert Jan van Oldenborgh for being the main guinea pig during the early stages of the project, André Heck for writing a good quality tutorial, Walter Hoogland for allowing me to spend so much time on the project in its early stages when success had not been proven yet.

## Appendix A The packages

A few packages are provided with the FORM distribution. These packages are useful not only for their field of applicability. They also serve as examples of FORM packages and programming techniques. In this paper we took most of the examples from them. Here we give a very short description of what these packages are good for. More information is in a number of papers and in the future there may be some manuals.

### A.1 Summer

This packages deals with harmonic sums. A full description of these sums is given in ref. . We give a short overview. The harmonic series is defined by

$`S_m(n)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{1}{i^m}}`$ (1)
$`S_m(n)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{(1)^i}{i^m}}`$ (2)
in which $`m>0`$. One can define higher harmonic series by

$`S_{m,j_1,\mathrm{},j_p}(n)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{1}{i^m}}S_{j_1,\mathrm{},j_p}(i)`$ (3)
$`S_{m,j_1,\mathrm{},j_p}(n)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{(1)^i}{i^m}}S_{j_1,\mathrm{},j_p}(i)`$ (4)
with the same conditions on $`m`$. The $`m`$ and the $`j_i`$ are referred to as the indices of the harmonic series. In the program we will use $`S_{j_1,\mathrm{},j_p}(n)`$ `S(R(j1,...,jp),n)`.

There is an alternative notation in which the indices are one of the values 0,1,-1. A zero indicates that actually one should be added to the absolute value of the nonzero index to the right of the zero as in:

$`S_{0,1,0,0,1,1}(n)`$ $`=`$ $`S_{2,3,1}(n)`$ (5)
The weight of a sum is sum of the absolute value of its indices. In the notation with the 0,1,-1 it is equal to the number of indices. The harmonic sums with the same argument form a (weight preserving) algebra: The product of two sums with weights $`w_1`$ and $`w_2`$ respectively is a sum of terms, each with a single sum of weight $`w_1+w_2`$. The routine ‘basis’ in the package implements this property. Example:

```
#-
#include summer6.h
.global
Local F = S(R(2,3),N)*S(R(-1,2),N);
#call basis(S)
Print +f;
.end

```

This program gives the output

```
   F = - S(R(-3,2,3),N) - S(R(-3,3,2),N) + S(R(-3,5),N)
       + 2*S(R(-1,2,2,3),N) + S(R(-1,2,3,2),N) - S(R(-1,2,5),N)
       - S(R(-1,4,3),N) - S(R(2,-4,2),N) + S(R(2,-1,2,3),N)
       + S(R(2,-1,3,2),N) - S(R(2,-1,5),N) + S(R(2,3,-1,2),N);

```

This defines a standard basis. A large number of sums involving harmonic sums and denominators can be done to any level of complexity. We give here a relatively simple example:

```
#-
#include summer6.h
.global
Local F = sum(j,1,N)*S(R(2,3),0,j)*S(R(-1,2),N,-j)*den(0,j);
#call summit()
Print +f;
.end

```

which gives the output

```
Time =       0.15 sec    Generated terms =        104
                F        Terms in output =         78
                         Bytes used      =       2938

   F = - S(R(-1,-3,2,3),N) + S(R(-1,-2,-2,-1,3),N) + S(R(-1,-2,-1,-2,3),N)
       - S(R(-1,-2,2,-4),N) + S(R(-1,-2,2,-3,-1),N) + S(R(-1,-2,2,-2,-2),N)
                ....
       + S(R(1,2,-1,-3,-2),N) + S(R(1,2,-1,2,-3),N) + 2*S(R(1,2,-1,3,-2),N)
       + 3*S(R(1,2,-1,4,-1),N) - 4*S(R(1,2,-1,5),N) + S(R(1,2,3,-1,2),N);

```

When $`N\mathrm{}`$ there are relations in addition to the regular algebraic relations. They allow all the sums in infinity to be expressed in terms of a rather small number of trancendental numbers of the Euler Zagier type . In order to obtain these expressions one has to solve for each weight a number of equations. This number grows exponentially with the weight. The system has been solved exactly up to weight 9 and the file with the solutions is part of the distribution. It is about 20 Mbytes. In the file ‘summer6.h’ only the solutions up to weight 6 have been included. The programs that solve these systems of equations are also part of the distribution. In the case of weight 9 there were more than 40000 equations in 13122 objects. The divergent sums were part of this because they can all be expressed properly into combinations of the basic divergence $`S_1(\mathrm{})`$. This database of sums is rather important for whole categories of integrals as we will see in the harmpol and meltran packages. It should be noted that at the moment of writing, no constructive algerithms are known to determine these constants one by one. The only ‘algorithm’ that comes close is the numerical evualuation of each individual sum, making an ansatz about which constants form a minimal set and then fitting all sums with a special numerical program . This has been done for $`w=10`$ and in a restricted way for some higher values of $`w`$ .

### A.2 Harmpol

This is a package for Harmonic Polylogarithms as defined in ref. . A short description is in order here.

The weight one harmonic polylogarithms are given by:

$`\text{H}(0;x)`$ $`=`$ $`\mathrm{ln}x,`$
$`\text{H}(1;x)`$ $`=`$ $`{\displaystyle _0^x}{\displaystyle \frac{dx^{}}{1x^{}}}=\mathrm{ln}(1x),`$
$`\text{H}(1;x)`$ $`=`$ $`{\displaystyle _0^x}{\displaystyle \frac{dx^{}}{1+x^{}}}=\mathrm{ln}(1+x).`$ (6)
For the higher weights we first define

$`f(0;x)`$ $`=`$ $`{\displaystyle \frac{1}{x}},`$
$`f(1;x)`$ $`=`$ $`{\displaystyle \frac{1}{1x}},`$
$`f(1;x)`$ $`=`$ $`{\displaystyle \frac{1}{1+x}}.`$ (7)
Then

$$\text{H}(\stackrel{}{0}_w;x)=\frac{1}{w!}\mathrm{ln}^wx,$$
(8)
with $`\stackrel{}{0}_w`$ an array with w zeroes while if the sequence $`\stackrel{}{m}_w`$ of w elements 0,1,-1 is not equal to $`\stackrel{}{0}_w`$

$$\text{H}(\stackrel{}{m}_w;x)=_0^x𝑑x^{}f(m_w;x^{})\text{H}(\stackrel{}{m}_{w1};x^{}).$$
(9)
The weight of the harmonic polylogarithms is the number of indices. The argument $`x`$ can in principle take any complex value, but depending on the indices there may be a cut for real values $`1`$ or $`0`$. These cuts are due to powers of $`\mathrm{ln}(1x)`$ or $`\mathrm{ln}(x)`$. Such powers can be extracted, resulting in combinations of these logarithms and H-functions without these cuts. This will be needed for the Mellin transforms. It should be noted that we can use also the notation of the harmonic sums for the indices (with values beyond the set 1,0,-1).

Also the harmonic polylogarithms form a weight preserving algebra. In addition there are extra relations for $`x=1`$. One has to be rather careful with the divergent elements at $`x=1`$. The powerseries expansion of H contains harmonic sums. The algebraic relations for general $`x`$ are related to the $`N\mathrm{}`$ relations of the harmonic sums and the extra relations in $`x=1`$ are related to the general algebraic relations for the harmonic sums.

The transformations $`x1x`$, $`x1/x`$ and $`x\frac{1x}{1+x}`$ give some interesting relations that allow for a reasonably fast numerical evaluation of the function over most of the complex plane (excepting the regions around $`\pm i`$).

The values in $`x=1`$ can be expressed in terms of the harmonic sums in $`N\mathrm{}`$. Because a weight w H-function is an integral over a denominator and a weight $`w1`$ H-function this gives us whole categories of integrals that were till now extremely hard to solve. The results will be in terms of a limited number of trancendental numbers, but the relations for $`x1x`$ at $`x=\frac{1}{2}`$ will allow them to be expressed in terms of relatively fast converging series. Hence these numbers can be evaluated to any precision when the need arises.

Because the harmpol package uses a number of the routines of the summer package it will automatically load this package. Hence once harmpol.h is included, so is summer6.h. Example:

```
#-
#include harmpol.h
Off statistics;
.global
Local F = H(R(1,0,1),x)*H(R(-1,1,-1),x);
#call hbasis(H,x)
repeat id H(R(?a,n?!{1,0,-1},?b),x?) = H(R(?a,0,n-sig_(n),?b),x);
.sort
On statistics;
Print +f;
.end

```

gives:

```
    #-

Time =       0.78 sec    Generated terms =         14
                F        Terms in output =         14
                         Bytes used      =        468

   F =
      H(R(-1,1,-1,1,0,1),x) + H(R(-1,1,0,1,-1,1),x) + 2*H(R(-1,1,0,1,1,-1),x)
       + 2*H(R(-1,1,1,-1,0,1),x) + 2*H(R(-1,1,1,0,-1,1),x) + 2*H(R(-1,1,1,0,1,
      -1),x) + H(R(1,-1,0,1,-1,1),x) + 2*H(R(1,-1,0,1,1,-1),x) + H(R(1,-1,1,-1
      ,0,1),x) + H(R(1,-1,1,0,-1,1),x) + H(R(1,-1,1,0,1,-1),x) + H(R(1,0,-1,1,
      -1,1),x) + 2*H(R(1,0,-1,1,1,-1),x) + H(R(1,0,1,-1,1,-1),x);

```

The conversion of combinations of denominators, logarithms and polylogarithms in terms of H-functions is usually a straightforward exercise in integration by parts, provided that it can be done at all. We give one example of an integral

```
#-
#include harmpol.h
#call htables(6)
CF  int;
Off Statistics;
.global
Local F = ln_(1+x)^2*ln_(1-x)^2*ln_(x)/x*int(x,0,1);
id  ln_(1+x) =  H(R(-1),x);
id  ln_(1-x) = -H(R(1),x);
id  ln_(x)   =  H(R(0),x);
#call hbasis(H,x)
id  H(R(?a),x)/x*int(x,0,1) = H(R(0,?a),1)-H(R(0,?a),0);
id  H(R(?a,n?{1,-1},?b),0) = 0;  * Provided there is a nonzero index
#do i = 1,6
id  H(R(n1?,...,n‘i’?),1) = Htab‘i’(n1,...,n‘i’);
#enddo
.sort
On Statistics;
Print +f;
.end

```

which gives

```
    #-

Time =       0.61 sec    Generated terms =         10
                F        Terms in output =         10
                         Bytes used      =        172

   F = - 1/2*z2*ln2^4 - 129/140*z2^3 + 7/6*z3*ln2^3 - 37/16*z3^2
       - 31/8*z5*ln2 + 8*ln2*li5half + 4*ln2^2*li4half
       + 1/9*ln2^6 + 8*li6half + 2*s6;

```

The constants are explained both in the paper and the htable files. These integrals are not easy by conventional techniques.

### A.3 Meltran

This package does the Mellin transforms and inverse Mellin transforms between the harmonic polylogarithms and the harmonic sums. It allows calculations to be done in either space and then converted to the other, whatever is easiest. The definition of a Mellin transform is:

$`M(f(x),N)`$ $`=`$ $`{\displaystyle _0^1}𝑑xx^Nf(x)`$
$`M({\displaystyle \frac{f(x)}{(1x)_+}},N)`$ $`=`$ $`{\displaystyle _0^1}𝑑x{\displaystyle \frac{x^Nf(x)f(1)}{1x}}`$
$`M({\displaystyle \frac{f(x)\mathrm{ln}^p(1x)}{(1x)_+}},N)`$ $`=`$ $`{\displaystyle _0^1}𝑑x{\displaystyle \frac{\left(x^Nf(x)f(1)\right)\mathrm{ln}^p(1x)}{1x}}`$ (10)
in which the function $`f`$ is supposed to be finite for $`x=1`$ when the factor $`1/(1x)_+`$ is present. Hence we see that the first step in the Mellin transformation of a H-function is the extraction of potential powers of $`\mathrm{ln}(1x)`$. The Mellin transform of a H-function, divided by either $`1x`$ or $`1+x`$ can be solved recursively and the result is a combination of harmonic sums in the Mellin parameter $`N`$. There is a one to one correspondence between the Mellin transform of H-functions of weight $`w`$, divided by either $`1x`$ or $`1+x`$ and harmonic sums of weight $`w+1`$ (the Mellin transform can also contain sums of a lower weight, but of the highest weight there is only a single term). Hence the inverse Mellin transform can also be obtained rather easily. The only complicating factor is that there will be terms in $`x=1`$ or $`N\mathrm{}`$. Currently these can only be combined into a ‘minimal’ representation if the weights are at most 9.

The mellin package is part of the harmpol package. Hence it is sufficient to include the file harmpol.h. Example:

```
#-
#include harmpol.h
Off statistics;
.global
L   F = S(R(1,2,1,2),N);
#call invmel(S,N,H,x)
.sort
On Statistics;
Print +f;
.end

```

results in

```
    #-

Time =       2.72 sec    Generated terms =          4
                F        Terms in output =          4
                         Bytes used      =        144

   F =
       - 38/35*z2^3 + 1/2*z3^2 - 1/2*H(R(1),x)*[1-x]^-1*z2^2
       - H(R(1,0,1,1,0),x)*[1-x]^-1;

```

Of course it should be clear that with these techniques one can do many integrals involving combinations of powers of $`x`$, denominators with $`1\pm x`$, $`\mathrm{ln}(x)`$, $`\mathrm{ln}(1\pm x)`$ and polylogarithms with a large variety of arguments. This is usually a soft spot in most computer algebra systems.

### A.4 Color

This package handles group invariants. The theory behind this is given in ref. . All internal variables in this package have a name that starts with the characters cOl in order to avoid potential name conflicts. With this package one can, for instance, compute the group invariants that are needed in complicated Feynman diagram calculations in such a way that the values for a specific representation of a specific group can be substituted at a later stage. This makes calculations more general. A simple example of the use of this package is the calculation of a diagram with 14 vertices each having three legs, connected in such a way that the smallest loop contains 6 vertices. This diagram is unique. If the vertices all belong to the adjoint representation the program would be:

```
#include color.h
Off Statistics;
.global
G   g14 = cOlf(cOli1,cOli2,cOli3)*cOlf(cOli1,cOli4,cOli5)
         *cOlf(cOli2,cOli6,cOli7)*cOlf(cOli3,cOli8,cOli9)
         *cOlf(cOli4,cOli10,cOli11)*cOlf(cOli5,cOli12,cOli13)
         *cOlf(cOli6,cOli14,cOli15)*cOlf(cOli7,cOli16,cOli17)
         *cOlf(cOli8,cOli18,cOli19)*cOlf(cOli9,cOli20,cOli21)
         *cOlf(cOli10,cOli21,cOli15)*cOlf(cOli13,cOli19,cOli14)
         *cOlf(cOli17,cOli11,cOli18)*cOlf(cOli12,cOli16,cOli20);
sum cOli1,...,cOli21;
.sort
#call docolor
.sort
On Statistics;
Print +f +s;
.end

```

and its output is

```
Time =       1.93 sec    Generated terms =          3
              g14        Terms in output =          3
                         Bytes used      =         82

   g14 =
       + 1/648*cOlNA*cOlcA^7
       - 8/15*cOld444(cOlpA1,cOlpA2,cOlpA3)*cOlcA
       + 16/9*cOld644(cOlpA1,cOlpA2,cOlpA3)
      ;

```

This was run on a pentium 300. Note that the customary f has to be given as cOlf etc. Similarly the regular $`C_A`$ is given as cOlcA. The values for the invariants for a given group can be found in the paper. As an example we give here the value of E8:

```
*
*   Numbers for exceptional algebras
*
S   [NA+2],eta;
id  cOld644(cOlpA1,cOlpA2,cOlpA3) = 175/48*cOlcA^7*cOlNA/[NA+2]^2;
id  cOld444(cOlpA1,cOlpA2,cOlpA3) = cOlcA^6*cOlNA/[NA+2]^2
                                    *(125/27+125/216*cOlNA);
*
*   And for E8
*
id  cOlNA = 248;
id  1/[NA+2] = 1/250;
id  cOlcA = 30*eta;
Print +f +s;
.end

```

which gives

```
   g14 =
       + 2075760000*eta^7
      ;

```

with the value of $`\eta `$ usually being two.

### A.5 Mincer

This is a pure particle physics package. It concerns the evaluation of massless propagator type integrals with up to three loops. In addition there are routines for the expansion in (fixed) Mellin moments. Computation time for higher Mellin moments can be rather large. In addition the disk requirements increase exponentially with the number of the moment. One needs to specify the topology of the diagram and a few other options. A simple example is given by

```
#-
#define SCHEME "0"
#define TOPO "la"
V   p1,...,p8,Q;
#include mincer.h
Off Statistics;
.global
L   Dia = 4*p7.p8^2/p1.p1/p2.p2/p3.p3/p4.p4^2/p5.p5/p6.p6/p7.p7/p8.p8*Q.Q;
.sort
#call integral(‘TOPO’)
.sort
On Statistics;
Multiply ep^3;
id  ep = 0;
print;
.end

```

which results in

```
    #-
~~~Answer in MS-bar

Time =       0.85 sec    Generated terms =          5
              Dia        Terms in output =          5
                         Bytes used      =         64

   Dia =
       - 139/2 - 5/2*ep^-3 - 19/2*ep^-2 - 55/2*ep^-1 + 5/2*z3;

```

For more details one should have a look at the paper that comes with the package. This mincer package has been in use by the author for 9 years now. Over the years it has steadily been improved and optimized. The code is basically code for version 2 of FORM. It is this package that was used to test the parallel version.


text:
# A Neutron-Star Identification for the High-Energy Gamma-ray Source 3EG J1835+5918 Detected in the ROSAT All-Sky Survey

## 1. Introduction

The nature of the persistent high-energy ($`>100`$ MeV) $`\gamma `$-ray sources in the Galaxy remains an enigma three decades after their discovery. Ever since the identification of the mysterious $`\gamma `$-ray source Geminga as the first radio quiet but otherwise ordinary pulsar (see review by Bignami & Caraveo 1996), it has been argued that rotation-powered pulsars likely dominate the Galactic $`\gamma `$-ray source population (e.g., Yadigaroglu & Romani 1997), and that many of them will be radio quiet. That Geminga is a bright EGRET source at a distance of only $`160`$ pc (Carveo et al. 1996) begs for it not to be unique.

Several excellent neutron-star candidates for EGRET sources have emerged from X-ray observations in recent years. One is present in the supernova remnant CTA1 which is close to 3EG J0010+7309 (Brazier et al. 1998). Another is in the $`\gamma `$-Cygni supernova remnant which is coincident with 3EG J2020+4017 (Brazier et al. 1996). An intriguing 34 ms X-ray pulsar with a Be star companion lies in the error circle of 3EG J0634+0521 (Kaaret et al. 2000). X-ray nebulae that are inferred to contain pulsars have been detected in the error circles of the EGRET sources 2EG J1811–2339 (Oka et al. 1999) and 3EG J1420–6038 (Roberts & Romani 1998; Roberts et al. 1999). Radio-quiet neutron stars have been found which are not apparently EGRET sources, but they are important as members of the growing class of cooling neutron stars (RX J185635–3754, Walter 2000; PKS 1209–52, Zavlin et al. 2000). Neutron-star candidates discovered by ROSAT were reviewed by Caraveo, Bignami, & Trümper (1996) and by Motch (2000).

As the brightest of the as-yet unidentified intermediate-latitude EGRET sources at $`(\mathrm{},b)=(89^{},25^{})`$, 3EG J1835+5918 is a good a priori target for the next such identification. We recently completed an exhaustive search for a counterpart of 3EG J1835+5918. The observations included deep radio, X-ray, and optical surveys, as well as optical spectroscopic classification of every active object within or close to its 99% confidence error ellipse (Mirabal et al. 2000, hereafter Paper 1). In summary, we identified optically all but one of the ROSAT and ASCA sources in the region of 3EG J1835+5918 to a flux limit of $`5\times 10^{14}`$ erg cm-2 s-1, which is $`10^4`$ of the $`\gamma `$-ray flux, without finding any suggestive evidence for a possible counterpart among the identified sources. We also proposed that the one unidentified X-ray source, RX J1836.2+5925, is the most promising candidate for identification with 3EG J1835+5918 principally because the absence of an optical counterpart for it is the standard signature of an isolated neutron star, or perhaps a more exotic compact object.

The recent release of the decade-old ROSAT All-Sky Survey (RASS) data (Voges et al. 2000) provides an important new tool for detecting neutron stars by means of their soft, thermal X-ray emission, and thus for furthering the search for the counterparts of nearby EGRET sources. In this Letter, we report new results from the RASS on the previous candidate RX J1836.2+5925 which greatly enhance its credentials as the likely pulsar counterpart of 3EG J1835+5918. We also present a detailed study of the X-ray and optical objects in its immediate neighborhood which establishes that RX J1836.2+5925 is indeed undetected to a limit of $`V>25.2`$, as would be expected for a neutron-star counterpart of 3EG J1835+5918.

## 2. X-ray Observations of RX J1836.2+5925

A total of four X-ray observations that cover the entire 99% error ellipse of 3EG J1835+5918 have been made

FIG. 1.— A photon map of the $`1^{}\times 1^{}`$ region around RX J1836.2+5925 from the ROSAT All-Sky Survey. Only events with pulse heights in the range 0.1–0.35 keV are shown. The net, background subtracted counts from RX J1836.2+5925 number $`22.3\pm 6.3`$.

to date, one by the RASS with the PSPC instrument, two by the ROSAT High Resolution Imager (HRI), and one by the ASCA GIS. The dates and exposure times of these observations are listed in Table 1. The results of the two pointed ROSAT observations and the one ASCA observation were presented in Paper 1, and they will not be repeated here except where relevant to the properties of RX J1836.2+5925.

In the RASS Faint Source Catalog, we find an entry at the position (J2000) $`18^\mathrm{h}36^\mathrm{m}13^\mathrm{s}.59,+59^{}25^{}30^{\prime \prime }.5`$ with an uncertainty of $`26^{\prime \prime }`$ and a count rate of $`0.0146\pm 0.0041`$ s-1. The effective exposure time at this location was 1532 s. The hardness ratio of the source is given as $`1.0`$, meaning that all of its photons fall at energies below 0.4 keV. This source is consistent in position and flux with the later observed ROSAT HRI source RX J1836.2+5925 as described in Paper 1, and since it is by far the brightest source within a region of radius $`10^{}`$ there is no reason to suspect confusion with the many fainter HRI sources in its vicinity. An image of the source and its surrounding field in the RASS is shown in Figure 1.

A deep X-ray image of the field around RX J1836.2+5925, taken from the longer ROSAT HRI observation, is shown in Figure 2. RX J1836.2+5925 is the brightest source in this Figure. Six of the fainter X-ray sources were identified using multicolor CCD imaging and optical spectroscopy. Their positions and counterparts are given in Table 1 of Paper 1. We use these sources to tie the X-ray positions to the optical astrometric reference frame employed by the USNO–A2.0 catalog (Monet et al. 1996). Figure 3 shows the offsets between the X-ray positions and the optical positions before and after a zero-point shift was applied to the X-ray positions in order to correct a small systematic error which is hardly apparent in the individual identifications. This offset is only $`1^{\prime \prime }.3`$ in right ascension and $`0^{\prime \prime }.8`$ in declination. The final error circle

FIG. 2.— A $`10^{}\times 10^{}`$ portion of the ROSAT HRI image centered on the brightest, unidentified X-ray source RX J1836.2+5925 (large circle). As refined with the help of Figure 3, the position of RX J1836.2+5925 is (J2000) $`18^\mathrm{h}36^\mathrm{m}13^\mathrm{s}.77,+59^{}25^{}30^{\prime \prime }.4`$.

radius of $`3^{\prime \prime }`$ includes all of the X-ray source identifications in Figure 3, and it is undoubtedly a conservative estimate of the positional error of the bright unidentified source RX J1836.2+5925. The corrected position for the latter is (J2000) $`18^\mathrm{h}36^\mathrm{m}13^\mathrm{s}.77,+59^{}25^{}30^{\prime \prime }.4`$. We adopt this as the best X-ray position. It differs by only $`1^{\prime \prime }.5`$ from the position given in Paper 1, and it agrees very well with the position of the corresponding RASS source.

We further examined the nature of RX J1836.2+5925 in the RASS by extracting source photons from within a radius of $`3^{}`$, and obtaining a background estimate from an annulus of inner radius $`4^{}`$ and outer radius $`8^{}`$. Although this source contains only 22 net photons it is clearly real; Voges et al. (2000) give a probability of $`3\times 10^7`$ that it is spurious. By binning the photons into pulse-height intervals of 0.1 keV, we see that there is no evidence for any emission above 0.3 keV (see Table 2). If fitted by a blackbody model, such a pulse-height distribution is consistent with $`T5\times 10^5`$ K, but it is also dependent upon the unknown intervening column density. If we assume $`1\times 10^{20}<N_\mathrm{H}<3\times 10^{20}`$ cm-2, the bolometric flux corresponding to an assumed $`T=5\times 10^5`$ K is in the range $`(1.55.7)\times 10^{13}`$ ergs cm-2 s-1. The soft spectrum of this source also explains why it was not detected in the ASCA observation, which had a 1–10 keV detection limit of $`1\times 10^{13}`$ ergs cm-2 s-1.

The net count rate of RX J1836.2+5925 in the ROSAT HRI observation of 1998 is $`(2.25\pm 0.23)\times 10^3`$ s-1. For the spectral parameters assumed above, the HRI flux agrees with the PSPC flux measured 7.5 yr earlier to better than 10%. Although this is consistent with a thermal neutron star interpretation, the source was possibly not detected in the shorter HRI observation in 1995. The formal count rate in the latter observation was only $`(0.90\pm 0.51)\times 10^3`$ s-1. Because of the small number of photons detected, the significance of the variability should be evaluated using Poisson statistics. In this analysis we

FIG. 3.— The result of using six optically identified, fainter HRI sources from Figure 2 to tie the X-ray positions to the optical astrometric reference frame. Top: The uncorrected position offsets. Bottom: The final offset between X-ray and optical position after applying a zero-point shift of $`1^{\prime \prime }.3`$ in right ascension and $`0^{\prime \prime }.8`$ in declination. Error bars are the 90% X-ray statistical uncertainties. We adopt a conservative radius of $`3^{\prime \prime }`$ for the uncertainty in position of RX J1836.2+5925.

use pulse-height channels 1–9, as there are virtually no X-ray photons in higher channels in the HRI. If we assume a constant source given by the count rate in the longer HRI observation, then we should expect to find in the source detection circle in 1995 a total of 24 photons (20 source and 4 background), whereas a total of only 12 were detected. The Poisson probability of obtaining 12 or fewer events when 24 are expected is $`5.4\times 10^3`$. Thus, the source is variable at the 99.5% confidence level. However, the existence of a ROSAT PSPC flux which is consistent with the latest HRI detection may be regarded as evidence contradictory to this indication of variability, and we would not be surprised if the one weak HRI detection turns out to be a statistical anomaly. Furthermore, the source was persistent during the 1 month span of the HRI observation in which it was detected. The existence or not of variability is crucial to the detailed interpretation of the physics of this source, as we discuss below, but not to its identification with 3EG J1835+5918.

## 3. Optical Observations of RX J1836.2+5925

Armed with the precise X-ray position of RX J1836.2+5925, we reexamined the deepest optical images which we obtained on the MDM Observatory 2.4m telescope. Figure 4 is a reproduction of the $`V`$-band image from Paper 1, a total of 2 hr of exposure obtained on 2000 July 24. The adopted $`3^{\prime \prime }`$ radius error circle superposed. The circle is blank to a $`3\sigma `$ limit of $`V>25.2`$. A limit of $`R>24.5`$ was also obtained on 2000 July 15 (see Paper 1).

In order to guard against the possibility of an anomalous systematic error in the X-ray position, we obtained optical spectra of the three nearest objects to the west of the error circle on the Palomar 5m telescope. These objects are a late K star of magnitude $`V=20.7`$, and two faint galaxies which show no evidence of activity in their optical spectra. They are therefore not viable candidates for identification with RX J1836.2+5925. We have also determined that no optical object in this field shows proper motion which could account for its positional discrepancy with the X-ray source. Thus, RX J1836.2+5925 remains undetected optically to a limit of $`V>25.2`$, even allowing for a conservative uncertainty on its position. This upper limit implies that the ratio of X-ray-to-optical flux $`f_X/f_V`$ is greater than 300, an extreme which is seen only among neutron stars.

FIG. 4.— A deep $`V`$-band image at the location of the unidentified X-ray source RX J1836.2+5925. This 2 hr summed exposure from the MDM 2.4m telescope in seeing of $`0^{\prime \prime }.8`$ yielded a $`3\sigma `$ detection limit of $`V=25.2`$. North is up, and east is to the left. The field is $`70^{\prime \prime }`$ across, and the ROSAT HRI error circle is drawn as a conservative $`3^{\prime \prime }`$ in radius as derived from Figure 3. The best position is (J2000) $`18^\mathrm{h}36^\mathrm{m}13^\mathrm{s}.77,+59^{}25^{}30^{\prime \prime }.4`$.

## 4. Interpretation and Conclusions

The latest analysis of the EGRET observations of 3EG J1835+5918 leads to the conclusion that it shows no evidence for long-term variability (Reimer et al. 2000). Its spectrum can be fitted by a relatively flat power law of photon index –1.7 from 70 MeV to 4 GeV, with a turndown above 4 GeV. Such temporal and spectral behavior are more consistent with a rotation-powered pulsar than a blazar, which is the other major class of EGRET source. The soft X-ray source RX J1836.2+5925 located in the error box of the 3EG J1835+5918 presents properties that resemble a neutron star, and could be the second case of a nearby, radio-quiet $`\gamma `$-ray pulsar.

The detection of RX J1836.2+5925 in the RASS 7.5 yr prior to the final HRI observation eliminates an alternative hypothesis mentioned in Paper 1, that it might have been a luminous soft X-ray transient caused by the tidal disruption of a star by a supermassive black hole in a distant galaxy, and unrelated to 3EG J1835+5918. Such events are expected (and observed) to last only a few months (e.g., Komossa & Bade 1999; Komossa & Greiner 1999), while this source has persisted for at least 7.5 yr. Furthermore, our deep optical imaging failed to find a candidate galaxy at this location.

The available data from radio through $`\gamma `$-rays are thus consistent with the hypothesis that RX J1836.2+5925 is the actual counterpart of 3EG J1835+5918 and a more distant or older cousin of the Geminga pulsar. The implications of a neutron-star interpretation of RX J1836.2+5925 were discussed at length in Paper 1. Although the intrinsic flux of the source is highly uncertain because of the unknown column density, the latter is probably not greater than $`3\times 10^{20}`$ cm-2 or we would not be seeing such a soft spectrum. Then RX J1836.2+5925 is evidently at least 10–40 times fainter than Geminga. Thus it is either more distant than Geminga ($`d>160`$ pc, Caraveo et al. 1996), or cooler ($`T<5\times 10^5`$ K, Halpern & Wang 1997). But a cooling neutron star of age $`<10^6`$ yr should not be too distant, because to place it $`>400`$ pc above the Galactic plane at $`b=25^{}`$ would require a kick velocity at birth of $`>500`$ km s-1. A reasonable upper limit on the distance is therefore 1 kpc, which we note implies an isotropic $`\gamma `$-ray luminosity of $`6\times 10^{34}(d/1\mathrm{kpc})^2`$ ergs s-1, rather more than the spin-down power $`I\mathrm{\Omega }\dot{\mathrm{\Omega }}`$ of Geminga ($`3.3\times 10^{34}`$ ergs s-1). If the actual spin-down power of 3EG J1835+5918 is determined, it will place an additional constraint on the distance and beaming factor.

The only observation that possibly complicates this interpretation is the marginally significant long-term (and non-monotonic) variability by about a factor of two. Such behavior is not consistent with surface thermal radiation of the heat of formation from the interior. However, it is possible that thermal emission from older neutron stars is dominated by the heated polar caps impacted by inflowing relativistic particles from the acceleration regions which are responsible for their pulsar action. Indeed, such reheating has been shown to be the source of the X-rays from the nearby recycled millisecond pulsar PSR J0437–4715 (Zavlin & Pavlov 1998). If under such conditions the particle acceleration regions in the pulsar were intermittent, the thermal X-rays resulting from this external heating would also be variable. In principle, 3EG J1835+5918 could be as powerful as Geminga at $`\gamma `$-ray energies, but less conspicuous at soft X-ray wavelengths because its surface has cooled except for the $`1`$ km radius polar caps which are being reheated directly by the accelerator. The most efficient $`\gamma `$-ray pulsars are expected to be the ones like Geminga, operating closest to their death lines (Chen & Ruderman 1993). A millisecond pulsar typically possesses spin-down power and magnetospheric gap voltage similar to Geminga and other ordinary pulsars. Since there is also evidence that the magnetospheric X-rays from Geminga are variable (Halpern & Wang 1997), either an exotic middle-aged pulsar or a maximally efficient recycled pulsar are plausible scenarios for 3EG J1835+5918.

We thank Peter Mao and Fiona Harrison for obtaining the optical spectra mentioned in §3. This work was supported by NASA grant NAG 5–9095.


text:
# A Measurement of the Branching Ratio of 𝐾_𝐿→𝑒⁺𝑒⁻𝛾𝛾

## Abstract

We report on a study of the decay $`K_Le^+e^{}\gamma \gamma `$ carried out as a part of the KTeV/E799 experiment at Fermilab. The 1997 data yielded a sample of 1543 events, including an expected background of $`56\pm 8`$ events. An effective form factor was determined from the observed distribution of the $`e^+e^{}`$ invariant mass. Using this form factor in the calculation of the detector acceptance, the branching ratio was measured to be $`(K_Le^+e^{}\gamma \gamma ,E_\gamma ^{}>5\text{MeV})=(5.84\pm 0.15(\mathrm{stat})\pm 0.32(\mathrm{sys}))\times 10^7`$.

The process $`K_Le^+e^{}\gamma \gamma `$ occurs mainly through kaon Dalitz decay ($`K_Le^+e^{}\gamma `$) with an internally radiated photon. As such, measurements of the radiative decay can test radiative correction calculations and probe the form factor of the $`K_L\gamma \gamma ^{}`$ vertex. The study of $`K_Le^+e^{}\gamma \gamma `$ is also important because it is a significant background to searches for $`K_L\pi ^0e^+e^{}`$, a mode which has an important direct $`CP`$-violating component and was a primary focus of the KTeV/E799 rare-decay program .

The best existing measurements of the branching ratio $`(K_Le^+e^{}\gamma \gamma ,E_\gamma ^{}>5\text{MeV})`$ were made by the E799-I and NA31 experiments, which measured $`(6.5\pm 1.2(\mathrm{stat})\pm 0.6(\mathrm{sys}))\times 10^7`$  and $`(8.0\pm 1.5(\mathrm{stat})_{1.2}^{+1.4}(\mathrm{sys}))\times 10^7`$ , respectively, based on signals of 58 and 40 events. This Letter presents a branching ratio measurement with much higher precision, using data collected by the KTeV/E799 experiment at Fermilab in 1997.

Figure 1 shows a plan view of the KTeV/E799 detector. The Tevatron at Fermilab made an 800 GeV proton beam, typically delivering $`3.5\times 10^{12}`$ protons per minute. The beam was incident on a BeO target. Photons from the target were converted by a lead absorber immediately downstream. Charged particles were then removed with magnetic sweeping. Collimators defined two neutral beams that entered the KTeV/E799 apparatus 94 m downstream of the target. A 65 m vacuum decay region extended to the first drift chamber. The charged-particle spectrometer consisted of a dipole magnet and four drift chambers (two upstream and two downstream of the magnet) with $`100\mu \mathrm{m}`$ position resolution in both horizontal and vertical directions. The magnetic field imparted a $`205\text{MeV/c}`$ horizontal momentum kick. The spectrometer had a momentum resolution of $`\sigma (P)/P0.38\%0.016\%P`$, where $`P`$ is in GeV/c.

Photons were detected using an electromagnetic calorimeter, which consisted of 3100 pure CsI crystals, each 50 cm long . Crystals in the central section of the calorimeter had a cross-sectional area of $`2.5\times 2.5\text{cm}^2`$, and those in the outer region had a $`5\times 5\text{cm}^2`$ area. The calorimeter’s energy resolution for photons was $`\sigma (E)/E0.45\%2\%/\sqrt{E}`$, where $`E`$ is in GeV. The position resolution was about 1 mm. By comparing spectrometer momentum to calorimeter energy, $`e^\pm `$ could be distinguished from $`\pi ^\pm `$. Nine photon veto assemblies (lead scintillator sandwiches) detected particles leaving the fiducial volume.

Additional $`e^\pm /\pi ^\pm `$ separation was provided by eight transition radiation detector (TRD) stations located downstream of the last drift chambers . Each TRD station consisted of a polypropylene felt radiator followed by two planes of multiwire proportional chamber (MWPC) filled with an 80/20% mixture of Xe/CO2 gas. The pulse heights from the 16 MWPC planes, compared with the pion pulse height distributions from a sample of $`K_L\pi ^\pm e^{}\nu `$ decays, were used to calculate the confidence level that a given track was a pion. Pion rejection factors over 200:1 were possible with a 90% electron efficiency.

The trigger system required hits in two scintillation hodoscopes just upstream of the CsI, as well as drift chamber hits consistent with two coincident charged particles passing through the detector. The trigger used in this study counted the number of isolated clusters of in-time energy over 1 GeV with a special processor ; at least four such clusters were required. The total energy deposited in the calorimeter had to be greater than 28 GeV. The scintillator hodoscope behind a lead wall downstream of the calorimeter vetoed events with hadronic showers in the calorimeter, and events with activity in the photon veto system were similarly discarded. Events which passed the hardware trigger requirements were read out into on-line CPUs, which performed basic event reconstruction and applied a few loose event-topology and particle-identification cuts to select events to be written to tape.

During offline event reconstruction, each event was required to contain two tracks from oppositely-charged particles originating from a common vertex. The reconstructed vertex had to be between 96 m and 158 m from the target. The tracks had to be at least 1 cm apart at the most upstream drift chamber, with an opening angle between the tracks of at least 2.24 mrad. To select electrons, each track was required to point to a calorimeter cluster with an energy equal, within $`\pm 5\%`$, to the track momentum measured by the spectrometer. The probability that a track was caused by a pion, formed by combining TRD chamber signals, was required to be less than $`4\%`$ for each track. The total momentum of all decay products had to be under 216GeV/c and the total energy greater than 33 $`\mathrm{GeV}`$. Events with a total of four clusters were analyzed as $`K_Le^+e^{}\gamma \gamma `$ candidates. Events with a total of five clusters were analyzed as $`K\pi ^0\pi _D^0`$ decays, where $`\pi _D^0`$ denotes a pion decaying to $`e^+e^{}\gamma `$. The latter mode was used both to check the simulation and to measure the total number of $`K_L`$ decays in the data sample.

In addition to the $`K_Le^+e^{}\gamma \gamma `$ signal mode, several potential background processes were considered. The $`K_Le^+e^{}\gamma `$ mode was a background when an additional photon was present. This photon could have been radiated by an $`e^\pm `$ while passing through matter in the detector (“external” radiation), or it could have been reconstructed from coincidental activity in the calorimeter not associated with the decay (“accidental” activity). $`K\pi ^0\pi _D^0`$ was a background when one photon was missed. Similarly, $`K_L\pi ^0\pi ^0\pi _D^0`$ could have been a background but was suppressed because missing three photons generally resulted in a reconstructed mass well below the $`K_L`$ mass. Background from $`K_L\pi ^\pm e^{}\nu `$ plus two additional accidental photons with the $`\pi ^\pm `$ misidentified as $`e^\pm `$ by the calorimeter, was eliminated using the TRDs.

Monte Carlo (MC) simulations of the detector were performed to calculate the acceptances for $`K_Le^+e^{}\gamma \gamma `$ and $`K\pi ^0\pi _D^0`$ and the misidentification rates for $`K_Le^+e^{}\gamma `$, $`K\pi ^0\pi _D^0`$ and $`K_L\pi ^\pm e^{}\nu `$. The effect of accidental clusters was simulated by overlaying on the MC events with data events taken with a random trigger that had a rate proportional to the beam intensity. The finite angle of external radiation was simulated. Simulations of Dalitz decays included radiative corrections to $`𝒪(\alpha _{EM}^2)`$ based on the work of Mikaelian and Smith . The radiative correction simulation also quantitatively predicted the ratios $`\mathrm{\Gamma }(K_Le^+e^{}\gamma )/\mathrm{\Gamma }(K_L\gamma \gamma )`$ and $`\mathrm{\Gamma }(K_Le^+e^{}\gamma \gamma )/\mathrm{\Gamma }(K_L\gamma \gamma )`$ and kinematic distributions such as photon energy in the $`K_L`$ center-of-mass, $`E_\gamma ^{}`$. The simulations of $`K_Le^+e^{}\gamma \gamma `$ and $`K_Le^+e^{}\gamma `$ used the Bergström, Massó, and Singer (BMS) form factor , which includes a parameter $`\alpha _K^{}`$ describing the relative strengths of the intermediate vector and pseudoscalar meson amplitudes in the $`K_L\gamma \gamma ^{}`$ vertex.

The number of $`K_L`$ decays was measured by the $`K\pi ^0\pi _D^0`$ decay mode, as in the $`K_L\pi ^0e^+e^{}`$ analysis . There were $`(2.642\pm 0.015(\mathrm{stat})\pm 0.024(\mathrm{sys})\pm 0.091(\mathrm{BR}))\times 10^{11}`$ decays in the fiducial region, where the third term is from uncertainties on the branching ratios of $`K_L\pi ^0\pi ^0`$ and $`\pi ^0e^+e^{}\gamma `$ .

Several cuts on reconstructed quantities were made on events with exactly four calorimeter clusters to identify $`K_Le^+e^{}\gamma \gamma `$ signal candidates. First, to ensure that candidate events did not include radiated photons below the nominal infrared cutoff of 5 MeV, $`E_\gamma ^{}`$ was calculated for each photon and was required to be at least 8 MeV (to allow for finite detector resolution). Since the MC sample was subjected to the same requirement, the acceptance correction procedure yielded a branching ratio valid for $`E_\gamma ^{}>5\text{MeV}`$, allowing direct comparison to theoretical predictions as well as to other published measurements.

In order to reduce backgrounds in which one or more particles are missing or mismeasured in the detector, or which involve accidental activity, a momentum balance cut was imposed. The square of the component of the total momentum of the decay perpendicular to a line drawn from the target to the vertex was required to be less than $`300(\text{MeV/c})^2`$ for $`K_Le^+e^{}\gamma \gamma `$ candidates.

The invariant mass of the decay, $`m_{ee\gamma \gamma }`$, was calculated assuming that the tracks were electrons. This mass was required to be within $`11\text{MeV/c}^2`$ of the neutral-kaon mass. Figure 2 shows the distributions of $`m_{ee\gamma \gamma }`$, after all other cuts, for data, signal, and the two main backgrounds. Values of $`m_{ee\gamma \gamma }`$ for $`K\pi ^0\pi _D^0`$ decays tended to be less than the kaon mass because one photon was missing. The excess in data at very low $`m_{ee\gamma \gamma }`$ was from $`K_L\pi ^0\pi ^0\pi _D^0`$ decays with three missing photons. The $`K_Le^+e^{}\gamma `$ decays often had a mass above the kaon mass because most such decays had an extra, accidental photon. The peak at the kaon mass in the $`K_Le^+e^{}\gamma `$ distribution is caused by events with external radiation.

The distributions of minimum angle, $`\theta _{min}`$, between any photon and any electron in the center of mass of the decay are plotted in Fig. 3 after all other cuts. The small angle events were dominated by $`K_Le^+e^{}\gamma `$ with an externally radiated photon. The large angle events were dominated by $`K\pi ^0\pi _D^0`$, where the electron directions were relatively uncorrelated with photon directions. Signal events were required to have $`0.05<\theta _{min}<0.63\text{rad}`$.

There were 1,543 $`K_Le^+e^{}\gamma \gamma `$ candidate events that satisfied all requirements. The distribution of $`e^+e^{}`$ mass in these events after background subtraction was used to determine the form factor. However, since the MC did not include radiative corrections to $`K_Le^+e^{}\gamma \gamma `$, an effective parameter, $`\alpha _{\mathrm{Eff}}`$, was measured rather than the true $`\alpha _K^{}`$. We simulated the effect of different effective form factors in the MC by reweighting $`K_Le^+e^{}\gamma \gamma `$ and $`K_Le^+e^{}\gamma `$ MC events with a ratio of the BMS form factor squared over the generated form factor (where $`\alpha _{\mathrm{Eff}}=0.28`$) squared. The best fit to the data was found at $`\alpha _{\mathrm{Eff}}=+0.016`$ with $`\chi ^2=4.3`$ for 9 degrees of freedom. Fig. 4 shows the ratio of data to MC for this value and for $`\alpha _{\mathrm{Eff}}=0.30`$.

The statistical uncertainty was found by using 22 subsets of MC events, each the same size as the data. The spread of different $`\alpha _{\mathrm{Eff}}`$ found by using each subset in place of the data was taken as the statistical uncertainty, 0.083. A sample of 12 million $`K_L\pi ^0\pi ^0\pi _D^0`$ decays was used to check the MC simulation of the acceptance as a function of $`m_{ee}`$. A ratio of the data over MC results in a slope of $`0.038\pm 0.096`$ $`(\mathrm{GeV}/\mathrm{c}^2)^1`$ and a $`\pm 1\sigma `$ range of this slope was fed back into the $`K_Le^+e^{}\gamma \gamma `$ simulation, resulting in a systematic uncertainty of $`\pm 0.032`$ on $`\alpha _{\mathrm{Eff}}`$. Uncertainty in our ability to reconstruct events with $`K_L\pi ^0\pi ^0\pi _D^0`$ decays, where the two (usually closely spaced) tracks both passed through one of the beam regions of the first drift chamber led us to assign an uncertainty to $`\alpha _{\mathrm{Eff}}`$ of $`\pm 0.027`$. Other systematic uncertainties, similar to those assessed below for the branching ratio, were evaluated, but the results were negligible. These systematic uncertainties are combined to give a total uncertainty of $`\pm 0.042`$ on $`\alpha _{\mathrm{Eff}}`$.

The acceptance for $`K_Le^+e^{}\gamma \gamma `$ with generated $`E_\gamma ^{}>5\text{MeV}`$ was found to be $`0.984\%`$ using MC reweighted to have $`\alpha _{\mathrm{Eff}}=+0.016`$. A Monte Carlo simulation that was normalized to the number of $`K_L`$ decays in the detector predicted 31 background events from $`K_Le^+e^{}\gamma `$ (and $`K_Le^+e^{}\gamma \gamma `$ with true $`E_\gamma ^{}<5\text{MeV}`$) and 25 background events from $`K\pi ^0\pi _D^0`$. Thus, the branching ratio of $`K_Le^+e^{}\gamma \gamma ,E_\gamma ^{}>5\text{MeV}`$ is $`(5.84\pm 0.15(\mathrm{stat}))\times 10^7`$.

Several sources of systematic uncertainty for the branching ratio measurement were considered. The normalization calculation contributed uncertainties as mentioned previously. A systematic uncertainty of $`\pm 0.49\%`$ was assigned in order for the MC to match the relative excess of simulated events at small angles (Fig. 3). Another uncertainty was assigned to account for uncertainty in the inefficiency for tracks in the neutral-beam region of the drift chambers (the beam region was $`0.965\pm 0.028`$ efficient), by reweighting MC events (signal, backgrounds, and normalization) based on the number of tracks in the beam region in each event; the resulting $`0.49\%`$ drop in branching ratio was taken as a symmetric uncertainty. Varying $`\alpha _{\mathrm{Eff}}`$ by our combined uncertainty of $`\pm 0.093`$ changed the acceptance and $`K_Le^+e^{}\gamma `$ background, and the resulting $`\pm 2.25\%`$ change in branching ratio was taken as an uncertainty. Finally, uncertainties due to inaccuracies in the MC simulation of the detector were estimated by varying one cut at a time over a reasonable range simultaneously for data, signal, backgrounds, and normalization. The most significant of these uncertainties were from the transverse-momentum cut, the cut on the opening angles between the the two tracks in the lab frame, and the $`m_{ee\gamma \gamma }`$ cut. When added in quadrature, all of the cut variations contributed a systematic uncertainty of $`\pm 3.92\%`$. The combined systematic uncertainty on $`(K_Le^+e^{}\gamma \gamma ,E_\gamma ^{}>5\text{MeV})`$ was $`\pm 5.44\%`$.

MC numerical integration of the tree-level QED partial width, using the BMS form factor with $`\alpha _K^{}=+0.016`$ and $`(K_L\gamma \gamma )`$ = $`(5.86\pm 0.15)\times 10^4`$ , predicts $`(K_Le^+e^{}\gamma \gamma ,E_\gamma ^{}>5\text{MeV})`$ = $`(5.70\pm 0.14)\times 10^7`$, based on $`(K_Le^+e^{}\gamma )`$ = $`(9.1\pm 0.5)\times 10^6`$, in agreement with the measurement.

In summary, we have determined the branching ratio of $`K_Le^+e^{}\gamma \gamma `$, with an infrared threshold of $`E_\gamma ^{}>5\text{MeV}`$, to be $`(5.84\pm 0.15(\mathrm{stat})\pm 0.32(\mathrm{sys}))\times 10^7`$. This calculation uses an effective form factor, which does not include effects from radiative corrections, with $`\alpha _{\mathrm{Eff}}=+0.016\pm 0.083(\mathrm{stat})\pm 0.042(\mathrm{sys})`$. This measurement is an improvement in precision over previous measurements, and agrees with the predicted value at tree level.

We gratefully acknowledge the support and effort of the Fermilab staff and the technical staffs of the participating institutions for their vital contributions. This work was supported in part by the U.S. Department of Energy, The National Science Foundation and The Ministry of Education and Science of Japan.


text:
# CLASSICAL ORDER PARAMETER DYNAMICS AND THE DECAY OF A METASTABLE VACUUM STATE11footnote 1Presented at Strong and Electroweak Matter (SEWM2000), Marseille, France, 14-17 June, 2000. Based on an article appeared in Phys. Rev. D62 (2000) 085013 written in collaboration with Szabolcs Borsányi, András Patkós, János Polonyi.

## 1 The model and the time-history of the order parameter (OP)

We study numerically the dynamics of a classical $`\mathrm{\Phi }^4`$ theory in $`2+1`$ dimensions governed by the discretized field equation of motion $`\mathrm{\Phi }_{\widehat{n}}(t+a_t)+\mathrm{\Phi }_{\widehat{n}}(ta_t)2\mathrm{\Phi }_{\widehat{n}}(t)+a_t^2(\mathrm{\Phi }_{\widehat{n}}+\mathrm{\Phi }_{\widehat{n}}^3h)\frac{a_t^2}{a^2}_i(\mathrm{\Phi }_{\widehat{n}+\widehat{i}}(t)+\mathrm{\Phi }_{\widehat{n}\widehat{i}}(t)2\mathrm{\Phi }_{\widehat{n}}(t))=0,`$ with initial conditions : $`\dot{\mathrm{\Phi }}_{\widehat{𝐧}}(t=0)=0`$, $`\mathrm{\Phi }_{\widehat{𝐧}}(t=0)=\mathrm{\Phi }_0+\xi _{\widehat{n}}\mathrm{\Phi }_1`$ where $`\xi _{\widehat{n}}`$ is an evenly distributed white noise in the range $`(1/2,1/2)`$ (all quantities are expressed in units of mass). The corresponding initial kinetic power spectrum is $`E_k(𝐤)\omega ^2(𝐤)=1+4(\mathrm{sin}^2(k_xa/2)+\mathrm{sin}^2(k_ya/2))/a^2`$. We have chosen $`\mathrm{\Phi }_0=0.815`$ and $`\mathrm{\Phi }_1=4/\sqrt{6}`$ which corresponds to a temperature value $`T_i=0.57`$ assuring that the system is in the broken symmetry phase. Our goal is to describe the evolution of the system only in terms of an effective equation for the OP $`\mathrm{\Phi }(t)=\frac{1}{V}_𝐧\mathrm{\Phi }_𝐧(t)`$.

Fig. 1 shows a typical OP-history together with the later time history of the OP mean square (MS)-fluctuation ($`\mathrm{\Phi }^2\mathrm{\Phi }^2`$) and of its third moment ($`(\mathrm{\Phi }\mathrm{\Phi })^3`$). In general one can distinguish five qualitatively distinct parts of the OP-trajectory that starts with large amplitude damped oscillations corresponding to the excitation of resonating modes, followed by a rather slow relaxation to a metastable state characterised by $`\mathrm{\Phi }0.72`$. The (quasi)thermal motion in the metastable state is followed by the metastable $``$ stable transition induced by the external field $`h`$, during which we can observe characteristic variations of the second and third moments. Quantitative interpretation of this variation will be presented in the following section. The last portion of the trajectory represents thermal motion in the true ground state.

## 2 The effective OP-theory

For the description of (quasi)thermal motion near a (meta)stable point we assume the validity of the ergodicity hypothesis for a system which consist of a single degree of freedom, the OP of the system. We describe its local time evolution by an effective Newton type equation: $`\ddot{\mathrm{\Phi }}(t)+\eta (\mathrm{\Phi })\dot{\mathrm{\Phi }}(t)+f(\mathrm{\Phi })=\zeta (t)`$. $`\eta (\mathrm{\Phi })`$, $`f(\mathrm{\Phi })`$, $`\zeta (t)`$ are obtained by a fitting procedure attesting in this way also the presence of a term violating time-reversal invariance $`\eta (\mathrm{\Phi })`$. $`\zeta (t)`$ is the “error” of the best global fit to the homogeneous equation at time $`t`$. The force felt by the OP, $`f(\mathrm{\Phi })`$ agrees with the force derived from the equilibrium two-loop effective potential as shown in Fig. 2. To probe this agreement in a relatively wide region one has to measure the force by shifting the center of motion to different values of $`\mathrm{\Phi }`$ by applying appropriate $`h`$ fields to the system. The coefficient $`\eta `$ is well-defined and positive, but its value depends quite substantially on the time resolution. The analysis of the relationship of time averaging to the nonzero value of $`\eta `$ is left for future investigations.

The thermalization “history” of the system is shown in Fig. 3. The displayed temperature variations correspond to the kinetic energies of the soft and hard modes. During the time evolution they approach each other, manifesting relaxation to a local thermal equilibrium. Because of our “white noise” initial condition the equilibration process is characterized by an energy transfer towards the low $`\underset{¯}{𝐤}`$ modes.

The fluctuation moments depicted in Fig. 1 along with the measured force shown on the right of Fig. 2 tell about how the transition proceeds. The increased values of the moments indicate the enhanced importance of the soft interactions. Preceding directly the transition towards the direction of negative $`\mathrm{\Phi }`$ values, the fitted force bends down and its average becomes a small positive constant. The vanishing of the force implies the flatness of the effective potential along the motion of the OP (the mode $`𝐤=0`$ in momentum space) indicating the dynamical realization of the Maxwell-cut.

Taking into account the existence of a mixed phase during the transition period we can construct a model that reproduces exactly the shape of the moments and shows the vanishing of the force felt by OP. Upon space averaging the microscopic field equation of motion of section 1 we find the following equation of motion for the OP: $`0=\ddot{\mathrm{\Phi }}\mathrm{\Phi }+\mathrm{\Phi }^3+3\overline{\phi ^2}^V\mathrm{\Phi }+\overline{\phi ^3}^Vh\ddot{\mathrm{\Phi }}+dV_{inst}/d\mathrm{\Phi }`$ where the symbol $`\overline{\phi ^n}^V`$ means the space average of $`\phi ^n`$, ($`\overline{\phi }^V=0,\mathrm{\Phi }(t,𝐱)=\mathrm{\Phi }(t)+\phi (𝐱,t)`$). The instant potential $`V_{inst}`$, which is a fluctuating quantity, contains as a deterministic piece the sum of the tree-level potential and of the slowly varying part of the second and third moments $`\overline{\phi ^n}^V(t)=(\overline{\phi ^n}^V)_{det}+\zeta _n(t)`$, $`n=2,3`$. We assume that the space can be splitted into sharp domains (neglecting the thickness of the walls in between), where the field is the sum of the constant background values $`\mathrm{\Phi }_{0\pm }`$ and the fluctuations $`\phi _\pm `$ around it, $`\mathrm{\Phi }_\pm (𝐱,t)=\mathrm{\Phi }_{0\pm }+\phi _\pm (𝐱,t).`$ Based on the smooth evolution of the temperature as displayed in Fig. 3 we assume also local equilibrium in both phases.

The actual value of the OP is determined by the surface ratio $`p(t)`$ occupied by the stable phase: $`\mathrm{\Phi }(t)=p(t)\mathrm{\Phi }_0+(1p(t))\mathrm{\Phi }_{0+}.`$ Simple calculation yields

$`(\overline{\phi ^2}^V)_{det}(t)=\frac{\mathrm{\Phi }_{0+}\mathrm{\Phi }(t)}{\mathrm{\Phi }_{0+}\mathrm{\Phi }_0}\left(\overline{\mathrm{\Phi }_{}^2(𝐱,t)}^V\overline{\mathrm{\Phi }_+^2(𝐱,t)}^V\right)+\mathrm{\Phi }_{0+}^2\mathrm{\Phi }^2(t)+\overline{\phi _{+}^{}{}_{}{}^{2}}^V,`$

$`(\overline{\phi ^3}^V)_{det}(t)=\frac{\mathrm{\Phi }_{0+}\mathrm{\Phi }(t)}{\mathrm{\Phi }_{0+}\mathrm{\Phi }_0}\left(\overline{\mathrm{\Phi }_{}^3(𝐱,t)}^V\overline{\mathrm{\Phi }_+^3(𝐱,t)}^V\right)+\overline{\mathrm{\Phi }_+^3(𝐱,t)}^V3\mathrm{\Phi }(t)(\overline{\phi ^2}^V)_{det}(t)\mathrm{\Phi }^3(t)`$.

If one takes the values of $`\mathrm{\Phi }_{0\pm },\overline{\phi _\pm ^n}^V`$ from the respective equilibria a quite accurate description of the shape of the two fluctuation moments arises in the whole transition region and its close neighbourhood using the measured $`\mathrm{\Phi }(t)`$ to parametrize their $`t`$-dependence. Substituting the expressions of the moments into the equation of the OP one finds for the deterministic part of the force, $`f\left(\mathrm{\Phi }\right)=\frac{\overline{\mathrm{\Phi }_+^3(𝐱,t)}^V\overline{\mathrm{\Phi }_{}^3(𝐱,t)}^V}{\mathrm{\Phi }_{0+}\mathrm{\Phi }_0}\mathrm{\Phi }\mathrm{\Phi }+\frac{\mathrm{\Phi }_{0+}\overline{\mathrm{\Phi }_{}^3(𝐱,t)}^V\mathrm{\Phi }_0\overline{\mathrm{\Phi }_+^3(𝐱,t)}^V}{\mathrm{\Phi }_{0+}\mathrm{\Phi }_0}h.`$ The average of the equations of motion in the respective equilibria, $`\mathrm{\Phi }_\pm ^3(𝐱,t)\mathrm{\Phi }_{0\pm }h=0`$ implies the vanishing of the deterministic force, when exploiting the equality of averaging over the volume and the statistical ensemble. The equation of motion for OP goes into $`\ddot{\mathrm{\Phi }}(t)+\zeta _3(t)+3\zeta _2(t)\mathrm{\Phi }(t)=0`$ which reflects the dynamical realization of the Maxwell construction when assuming the validity of local equilibrium in both phases.

## 3 The nucleation picture

The statistics of the release time, the time necessary for the system to escape from the metastable region by nucleating a growing bubble of the stable state, is of the form $`P(t)\text{exp}(t\mathrm{\Gamma }(h)L^2)`$ where $`\mathrm{\Gamma }(h)=A\text{exp}(S_2(h)/T)`$ is the nucleation rate and $`L`$ is the lattice size . We extract the exponent $`S_2`$ by fitting the measured rate to the expression of $`\mathrm{\Gamma }(h)`$ assuming the $`h`$-independence of $`A`$. Comparing with the bounce action of the nucleation theory we obtain: $`S_{\text{thin wall}}=\frac{4\pi }{9}hT15S_{\text{measured}}`$. Using the $`T0`$ effective potential instead of the bare one one can achieve a much better agreement, but still a factor of 2 discrepancy remains with the mean field approach.

## 4 Nucleation vs. spinodal instabilities

Our results offer a “dualistic” resolution of the competition between the nucleation and the spinodal phase separation mechanisms in establishing the true equilibrium. We find that the statistical features of the decay of the false vacuum agree with the results of thermal nucleation. Alternatively, the effective OP-theory displays the presence of soft modes and produces dynamically a Maxwell cut when the time dependence of the transition trajectory is described in the effective OP theory. During the transition the OP travels trough a narrow but flat valley around the $`𝐤=0`$ mode along which we expect the effective potential to be flat. This flatness is reflected also in the decrease of the kink-like action $`S_{measured}`$ relative to $`S_{thinwall}`$. The larger is the system the smaller is the external field which is able to produce the instability.

## References