CHIRON: a package for ChPT numerical results at two loops

This document describes the package CHIRON which includes two libraries, chiron itself and jbnumlib. chiron is a set of routines useful for two-loop numerical results in chiral perturbation theory (ChPT). It includes programs for the needed one- and two-loop integrals as well as routines to deal with the ChPT parameters. The present version includes everything needed for the masses, decay constants and quark-antiquark vacuum-expectation-values. An added routine calculates consistent values for the masses and decay constants when the pion and kaon masses are varied. In addition a number of finite volume results are included: one-loop tadpole integrals, two-loop sunset integrals and the results for masses and decay constants. The numerical routine library jbnumlib contains the numerical routines used in chiron. Many are to a large extent simple C++ versions of routines in the CERNLIB numerical library. Notable exceptions are the dilogarithm and the Jacobi theta function implementations. This paper describes what is included in CHIRON v0.50.


Introduction
Chiral perturbation theory (ChPT) is the low-energy effective field theory of QCD. It was introduced by Weinberg, Gasser and Leutwyler [1][2][3] and the present state of the art are calculations performed at two-loop level. A review is [4] but many more exists. The long term goal of this project is to make available all these calculations with a consistent interface in C++. Many of the original programs were written in FORTRAN77 and are available on request from the authors, but they are not always consistent in the interfaces and the use of common blocks for moving parameters around has occasionally lead to difficult to find errors.
A general knowledge of C++ is assumed throughout this paper. The routines are at present not guaranteed to be threadsafe, some global variables inside the various files are used for the loop functions and integration routines. These are however never used for setting outside the files, there are always functions provided for this. It is recommended to always use these. The routines return double precision types if not indicated directly.
Kheiron, X ειρων, or Chiron, was the eldest and wisest of the Centaurs, half-horse men of Greek mythology. His name is derived from the Greek word for hand (Kheir) which also formed the basis of the word chiral which is why his name was chosen for this package [5].
The license chosen is the General Public License v2 or later from the Free Software Foundation [6].
The chiron routines have mainly been tested against the FORTRAN codes of the original publications. These were in turn implemented in at least two independent versions originally. The jbnumlib routines have their output compared with the original CERNLIB routines in case they were simple translations to C++. In the other cases, the tests are described in the relevant sections.
This paper describes what is included in CHIRON v0. 50. The package itself is available from [7]. The included files and how to install it is described in Sect. 2. The numerical routines included in jbnumlib are described in Sect. 3. Some short comments on ChPT notation are in Sect. 4. The main part describing the contents of chiron are Sects. 5 to 8. Section 5 contains the objects implemented to deal with input data for the ChPT calculations. A large part of the work is in implementing the relavant loop integrals, especially those at finite volume. This is the content of Sect. 6. The simplest quantities are masses, decay constants and vacuumexpectation-values. The functions for these are discussed in Sect. 7 and the finite volume extensions for masses and decay constants in Sect. 8. Some comments about errors, some warnings about the use of the routines and definitions in ChPT as well as a number of planned/possible extensions are discussed in Sect. 9. A short summary is given in the final section.

Files and setup
The package is delivered as a gzipped tarred file (chironvvvv.tar.gz), where vvvv is version information. Untarring it creates a directory chironvvvv, which is referred to as the root directory below.
The package has a number of subdirectories when delivered. The root directory contains a Makefile and files COPYING, INSTRUCTIONS and GUIDELINES. The main instructions to produce CHIRON are to do first "make libjbnumlib.a" to produce the numerical library and then "make libchiron.a" to produce the main library or simply "make" to do both. The latter also puts the newly produced library versions in the subdirectory lib. To install, put the two library files and the content of the include directory somewhere where they can be linked to and included. Linking should be indicated on the compile line with the options "-lchiron -ljbnumlib".
The subdirectory doc contains this manuscript, a filelist and possibly more files in future versions. The subdirectory src contains the source files and include the various header files. The subdirectory test contains a number of testing programs where the names "testyyy.cc" indicates a program testing the code in the source file "yyy.cc". A testing program can be compiled using "make testyyy" in the root directory. The program "a.out" should then produce the output as shown in the file "testyyy.dat" in the subdirectory testoutputs. The test subdirectory contains in addition the file LiCiBE14.dat with the latest determination of the LECs [8].

jbnumlib
The functions in this section are included to make the program collection self-contained. They are mainly implementations of well known programs in C++ and in particular many of the routines are a port to C++ from the CERNLIB [9] FORTRAN routines. Some, as mentioned in the respective texts, are fully original. The definitions are all in jbnumlib.h and contained in libjbnumlib.a. The implementations are in the files mentioned in each subsection below. In order to avoid conflict with other implementations all routines in this section have names starting with jb. The exact interface is best checked by looking in the include file jbnumlib.h.

Special functions
The way the vertex integrals are implemented requires a Spence or Li 2 function which returns complex values for all possible complex inputs. The routine implemented uses the algorithm given in [10], Appendix A up to Bernouilly number B 28 . Defined in jbdli2.cc. For real numbers the output has been compared to that of the CERNLIB routine DDILOG. It has also been checked that the function satisfies a number of the relations between values with different arguments that were not used in its evaluation.

Bessel functions
The modified Bessel functions I 0 , I 1 , K 0 , K 1 with real arguments are available as jbdbesi0, jbdbesi1, jbdbesk0 and jbdbesk1. These are implemented in jbdbesio which is a simple port to C++ of the CERNLIB routines dbesi0,…. In addition the modified Bessel functions K 2 , K 3 , K 4 are available as jbdbesk2, jbdbesk3 and jbdbesk4. These are evaluated using the recursion relations from K 0 and K 1 .

Theta functions
The functions defined are related to the Jacobi theta functions. jbdtheta30(q) returns the function It uses the idea behind the CERNLIB routines DTHETA. For small q it simply sums the series (1) and for larger q it uses the modular invariance and a series in the changed variable instead. The accuracy has been checked by running both series too much higher orders and comparing the two results. Implemented in jbdtheta30.cc. A function without the 1 which is needed to keep accuracy for small q is available as jbdtheta30m1. Implemented in jbdtheta30m1.cc. jbdtheta32(q) returns the function For small q it simply sums the series (2) and for larger q it uses the modular invariance and the derivative of the series in the changed variable The accuracy has been checked by running both series too much higher orders and comparing the two results. Defined in jbdtheta32.cc. A function with n 4 multiplying q (n 2 ) is available as jbdtheta34 and implemented in jbdtheta34.cc.

Higher dimensional theta functions
There are higher dimensional generalizations of the theta functions. These satisfy a more general modular invariance which can be used to get much faster convergence series. The functions defined below use some of these possible optimizations.
The basic function is the 2-dimensional generalization This function is implemented as jbdtheta2d0 with arguments α, β, γ in jbdtheta2d0.cc. In addition the function with the one removed is also available as jbd theta2d0m1 implemented in jbdtheta2d0m1.cc.
The function with the exponential multiplied by n 2 1 is called jbdtheta2d02 and implemented in jbdthet a2d02.cc.

Integration routines
An adaptive gaussian quadrature routine jbdgauss and an adaptive integration routine jbdcauch, that integrates symmetrically around a singularity, are included. These are ports of the CERNLIB routines dgauss and dcauch respectively. Implementation is in jbdgauss.cc and jbdcauch.cc. The complex equivalent is jbwgauss in jbwgauss.cc.
The higher dimensional integration CERNLIB routine DADMUL based on [11] has been ported to C++ and a simple interface for two and three-dimensional integration implemented as jbdad2 and jbdad3. The code can be found in jbdadmul.cc.

ChPT notation
The notation used in ChPT is not fully unique. The notation used in this program collection is the main one used by the author and his collaborators. The main point to be observed is that chiron uses a normalization for the decay constants with F π ≈ 92 MeV.
For the low-energy-constants, we use the conventions of [3,12,13] with dimensionless renormalized couplings L r i and C r i . The lowest order couplings are denoted by F 0 and B 0 . The quark masses arem = m u = m d , note that we work in the isospin limit, and m s . The lowest order masses 1 are given by

Data structures
This section discusses the structures available for dealing with masses, F π , the L r i and C r i , and the subtraction constant μ. Note that μ is present in all three data structures and the user should make sure that their use is consistent. The defaultvalue mechanism in C++ has been used to define the values when the constructors are called with less than the full data needed.
These data structures are implemented as classes. Note that we assume dimensional units in GeV, but if all dimensional inputs are scaled accordingly the routines give the correct answer. However, typical precisions set are assuming ChPT applications with dimensional units in powers of GeV.

physmass: Masses
The physical masses are defined in a class physmass defined in inputs.h and implemented in inputs.cc.
These are the physical pion, kaon and eta mass, m π , m K , m η the physical pion decay constant, F π , and the subtraction point μ. The default constructor puts them all at some reasonable values but they can be created with any number of the inputs specified, starting from the left. In addition there are functions void setmpi(double) etc., defined that set one of the values only. These use the same default values.
The values can be obtained from the function void out( mpi, mk, meta, fpi, mu) that returns all of them via referenced doubles. Functions double getmpi(void) etc. are defined that return one of the values.
The output/input stream format is defined as well so cout << mass1 and cin >> mass1 make sense. The input stream should have the same format as the output stream produces. This works for all streams, not just the standard cout and cin.
A test for equality is defined which checks that all data members agree to 7 significant digits. So expressions like if(mass1 == mass2) can be used. This is relative precision, so it is assumed no mass is zero here. The reason for not using exact equality is that calculated masses might not be exactly the same using double precision variables.

Classes for the NLO LECS: Li
The class for dealing with the next-to-leading-order (NLO) low-energy constants (LECs) defined in [3] is named Li. This is implemented in Li.cc and defined in Li.h.
The Li class has 13 double precision variables to store the LECs L r 1 , . . . , L r 10 , H r 1 , H r 2 and the subtraction scale μ. It also contains a string with a name for the set of constants. The LECs default to zero, the scale to 0.77 and the name to "nameless Li." When the LECs are referred to with numbers, 11,12 correspond to H r 1 , H r 2 respectively. Typical declarations are: Li Li1; Li Lifitall(0.88e-3,0.61e-3,-3. 04e-3,0.75e-3,0.58e-3,0.29e-3, -0.11e-3,-0.18e-3,5.93e-3,0.,0.,0.,0.77, "fit All"); Operations defined on the Li: overloaded operators are defined such that sets of Li can be added or subtracted and multiplied by a number (double). The output/input stream format is defined as well so cout << Li1 and cin >> Li1 make sense. The input stream should have the same format as the output stream produces.
Member functions that can be used to set the parameters are setli which takes an integer and a double (in either order) as argument to set the corresponding LEC to the double, setmu which sets the scale 2 and setname which changes the name of the set of L r i . Output member functions exist to obtain a single LEC, out(int), or the 10 L r i , the 10 L r i and μ, the 12 LECs, the 12 LECs and μ, and the 12 LECs, μ and the name. These are all called out and return the results via a reference to 10, 11, 12, 13 double or 13 double and a string variable.
The member function changescale changes the scale μ and changes the L r i and H r i according to the scale dependence as obtained first in [3].
There are also three functions defined that return a set of random NLO LECs. These are Lirandom which gives each LEC a random value between ±1/(16π 2 ). LirandomlargeNc does the same but leaves L r 4 , L r 6 and L r 7 zero. Finally, LirandomlargeNc2 does the same but L r 4 , L r 6 and L r 7 get a random value between ±(1/3)/(16π 2 ). Note that 1/(16π 2 ) ≈ 0.0063 so the ranges include the values of the fitted L r i . The random numbers are generated using the system generator rand() so we recommend initializing using something like srand(time (0)). These latter functions were used in the random walks in the L r i in [14].

Classes for the NNLO LECS: Ci
The class for dealing with the next-to-next-to-leading-order (NNLO) low-energy constants (LECs) defined in [12] is named Ci. This is implemented in Ci.cc and defined in Ci.h. Note that this set of routines uses the convention where the C r i are dimensionless. The parameters in the Lagrangian have dimension mass (−2) but the definition of the subtracted C r i in [13] is dimensionless. Going from one-convention to the other is with the appropriate power 3  An additional constructor is provided that has as input the resonance parameters where the resonance model is the simple version described in Sect. 5 of [15].
Operations defined on the Ci are: overloaded operators are defined such that sets of Ci can be added or subtracted and multiplied by a number (double). The output/input stream format is defined as well so cout << Ci1 and cin >> Ci1 make sense. The input stream should have the same format as the output stream produces.
Member functions that can be used to set the parameters are setci which takes an integer and a double (in either order) as argument to set the corresponding LEC to the double, setmu which sets the scale 4  The member functions changescale(double,Li) and changescale(Li,double) change the scale μ and changes the L r i , H r i and C r i according to the scale dependence as obtained first in [13]. Note that the Li set here has also the scale and values changed accordingly, not only the C r i . There are also three functions defined that return a set of random NLO LECs. These are Cirandom which gives each LEC a random value between ±1/(16π 2 ) 2 . CirandomlargeNc does the same but leaves the large-N c suppressed constants zero. Finally, CirandomlargeNc2 does the same but the large-N c suppressed constants get a value between ±(1/3)/(16π 2 ) 2 . The random numbers are generated using the system generator rand() so we recommend initializing using something like srand(time (0)). Typically these values of the C r i are somewhat on the large side when fitting data.

Loop integrals
Most of the integrals used have been treated in many places. I refer only to the papers where our particular notation has been defined and/or the method used to evaluate them was developed. When comparing with other packages, keep in mind the differences in subtraction and/or differences in defining the integrals. 6.1 One-loop integrals 6. 1

.1 Tadpoles
These integrals have been defined in [16] and correspond to the finite parts of the integral After the subtraction and renormalization as usual in ChPT, we are left with the finite four-dimensional partĀ(n, m 2 ) which is implemented as Ab(n,msq,mu2) and the n = 1, 2, 3 as Ab,Bb,Cb respectively with arguments msq,mu2.
These functions are defined in oneloopintegrals.h and implemented in oneloopintegrals.cc.

Bubble integrals
These have been defined in [16,17] Again one needs to do the subtraction and renormalization with ChPT convention. The analytical values can be obtained using the methods of [10] for the B integral and the others can be reduced to it using the methods of [18]. All functions have been implemented via a method that does the integration over the Feynman parameter x numerically. These have real arguments m1sq,m2sq,psq,mu2 and are called Bbnum,B1bnum,B21bnum,B22bnum,B31bnum,B32 bnum and return a complex value. The analytical evaluation has been implemented in Bb,B1b, B21b,B22b with the same arguments and a complex return value. The simpler analytical expression for the case of the two masses equal has been implemented analytically for Bb and B22b called with argument msq,psq,mu2 For the cases with numerical integrations, the precision can be set using setprecisiononeloopintegrals(double) and obtained by getprecisiononeloopintegrals(void).
All functions are defined in oneloopintegrals.h and implemented in oneloopintegrals.cc.

Sunset integrals
These give the sunsetintegral loop integral functions H F , H F 1 , H F 21 defined in App. C of [16]. The definition in finite volume is given in (9) below. H F 31 is the function multiplying the p μ p ν p ρ part of the integral with r μ r ν r ρ . These exist in a real version valid below threshold and a complex version valid everywhere. The method used is derived in [16].
The derivative w.r.t. p 2 is included for the real version of H , H 1 , H 21 . Functions defined in sunsetintegrals.h and implemented in sunsetintegrals.cc.
Input arguments are real and are m 2 1 , m 2 2 , m 2 3 , p 2 , μ 2 . Naming conventions are hh,hh1, hh21,hh31 for the real versions valid below threshold and the complex versions valid everywhere zhh,zhh1,zhh21,zhh31. The real versions are normally faster when applicable. In addition, wave-function renormalization requires some derivatives w.r.t. the external momentum p 2 . These are encoded in hhd,hh1d,hh21d and zhhd,zhh1d,zhh21d for the real and complex case respectively.
The precision of the numerical integrations can be set using setprecisionsunsetintegrals(double) and obtained by getprecisionsunsetintegrals(void).
These functions are defined in sunsetintegrals.h and implemented in sunsetintegrals.cc.

One-loop finite volume integrals
The methods used for these are derived in detail in [19], references to earlier literature can be found there. The integrals used here are given in the Minkowski conventions as defined in [20]. All of the integrals are available with two different methods, one using a summation over Bessel functions and the other an integral over a Jacobi theta function. The versions included at present are using periodic boundary conditions, all three spatial sizes of the same length L and the time direction of infinite extent.

Tadpoles
The tadpole integrals A and A μν are defined as The B tadpole integrals are the same but with a doubled propagator.
The subscript V on the integral indicates that the integral is a discrete sum over the three spatial components and an integral over the remainder. At finite volume, there are more Lorentz-structures possible. The tensor t μν , the spatial part of the Minkowski metric g μν , is needed for these. The functions for A μν are In infinite volume A 22 is related to A and A 23 vanishes. We denote the finite volume part by a superscript V and one should remember that for the full integrals, the infinite volume results of Sect. 6.1.1 need to be added. The functions are defined as AbVt(msq,L), BbVt (msq,L), AbVb(msq,L), BbVb(msq,L). The last letter indicates whether they are computed with the theta function or Bessel function method. The results were checked by comparing against each other and by comparing with the independent Bessel function implementation done in [21].
The functions A22bVt(msq,L), A22bVt(msq,L), and A23bVb(msq,L), A23bVb(msq,L) are available as well. setprecisionfinitevolumeoneloopt (Abacc,Bbacc,printout) and setprecisionfinitevolumeoneloopt(maxsum, Bbacc,printout) set the precision. The last variable printout is a logical variable which can be set to true or false, default is false. The first and second argument give the (mainly absolute) precision of the numerical integration for the tadpole and bubble integral numerical integrations. maxsum indicates how far the sum over Bessel function is taken. Maximum at present is 400.
These functions are defined in finitevolumeoneloop integrals.h and implemented in finitevolumeone loopintegrals.cc.

Bubble integrals
Not implemented in this version.

Sunset finite volume integrals
Sunset integrals are defined as The subscript V indicates that the spatial dimensions are a discrete sum rather than an integral. The conventions correspond to those in infinite volume of [16]. Integrals with the other momentum s in the numerator are related using the trick shown in [16] which remains valid at finite volume in the cms frame [19].
In the cms frame we define the functions 5 The arguments of all functions in the cms frame are (m 2 1 , m 2 2 , m 2 3 , p 2 ). These functions satisfy the relations, valid in finite volume [19], The arguments of the sunset functions in the second relation are all (m 2 1 , m 2 2 , m 2 3 , p 2 , L , μ 2 ). (L only for the finite volume part).
We split the functions in an infinite volume part,H i , and a finite volume correction, The infinite volume part has been discussed above. For the finite volume parts we definẽ The finite parts are defined differently from the infinite volume case in [16]. The parts with A V are removed here as well.
The functions H V i can be computed with the methods of [19]. They correspond to adding the parts labeled with G and H in Sect. 4.3 and the part of Sect. 4.4 in [19].
They are implemented as functions hhVt,hh1Vt, hh21Vt,hh22Vt,hh27Vt with arguments m1sq,m2sq, m3sq,psq,L,mu2. The derivatives w.r.t. p 2 exist as hhdVt,hh1dVt,hh21dVt,hh22dVt,hh27dVt. These are the functions using the theta function method. Those using the Bessel function method are implemented with a b instead of t as last letter in the name. The arguments are the same.
For all cases discussed we have done checks that both methods, via Bessel or (generalized) Jacobi theta functions, give the same results. In addition the derivatives w.r.t. p 2 for all the integrals are compared with taking a numerical derivative.
Note that the sunset functions at finite volume call the tadpole integrals evaluated with the same method. Do not forget to set precision for those as well. The precision for the sunset integrals can be set with the functions setprecisionfinitevolumesunsett(racc, rsacc,printout) and setprecisionfinite volumesunsetb(maxsum1,mxsum2,racc,rsacc, printout). The bool variable printout defaults to true and sets whether the setting is printed. The double values sunsetracc and sunsetrsacc set the accuracies of the numerical integration needed when one or two loop-momenta "feel" the finite volume. Default values are 1e-5 and 1e-4 respectively. The integers maxsum1 and maxsum2 give how far the sum over Bessel functions is used for the same two cases. The first is maximum 400, the second maximum 40 in the present implementation. In the latter case a triple sum is needed, hence the much lower upper bound. For most applications it makes sense to have a higher precision for the case with one loop momentum quantized, i.e. racc smaller than rsacc.

Masses
The masses of the pion, kaon and eta at two-loops in three flavour ChPT were calculated in [16]. The pion and eta mass were done earlier with a different subtraction scheme and a different way to perform the sunset integrals in [22]. The expressions for the physical masses for a = π, K , η are given by The superscripts indicate the order of the diagrams in p that each contribution comes from. The lowest order masses are given in (4). The expressions can be found in [16]. In addition the contributions themselves are split in the parts depending on the NLO LECs L r i , on the NNLO LECs C r i and the remainder as All the parts in (14) are implemented as the functions mpi4(physmass,Li), mpi4L(physmass,Li), mpi4R(physmass), mpi6(physmass,Li,Ci), mpi6L(physmass,Li), mpi6C(physmass,Ci) and mpi6R(physmass) . The equivalent functions also exist for the kaon, with pi to k, and eta, with pi to eta.
The functions are defined in massesdecayvev.h and implemented in massesdecayvev.cc.

Decay constants
The decay constants of the pion, kaon and eta at two-loops in three flavour ChPT were calculated in [16]. The pion and eta decay constants were done earlier with a different subtraction scheme and a different way to perform the sunset integrals in [22].
The expressions for the decay constants for a = π, K , η are given by The superscripts indicate the order of the diagrams in p that each contribution comes from. F 0 denotes the decay constant in the three-flavour chiral limit. The expressions were originally derived in [16], but note the description in the erratum of [15]. The expressions corrected for the error can be found in the website [23]. In addition the contributions themselves are split in the parts depending on the NLO LECs L r i , on the NNLO LECs C r i and the remainder as a L + F (6) a C + F (6) a R . (16) All the parts in (16) are implemented as the functions fpi4(physmass,Li), fpi4L(physmass,Li),fpi 4R(physmass), fpi6(physmass,Li,Ci), fpi6L (physmass,Li), fpi6C(physmass,Ci) and fpi6R (physmass) . The equivalent functions also exist for the kaon, with pi to k, and eta, with pi to eta. For the η the decay constant has been defined with the octet axial-vector current.
The functions are defined in massesdecayvev.h and implemented in massesdecayvev.cc.

getfpimeta
A problem that occurs in trying to compare to lattice QCD is that the present routines are written in terms of the physical pion decay constant and masses. In particular, the eta mass is treated as physical. One thus needs a consistent eta mass and pion decay constant when varying the input pion and kaon mass. This assumes we have fitted the LECs L r i and C r i with a known set of m π , m K , m η , F π .
The functions getfpimeta6(mpiin,mkin,massin, Li,Ci) and getfpimeta4(mpiin,mkin,massin, Li) return a physmass with a consistent set of F π and m η for input values of the pion and kaon mass. The other input is the physmass massin, the Li and Ci that are used as input. The formulas used are (14) and (16) up to order p 6 and p 4 respectively. The solution is obtained by iteration and stops when six digits of precision are reached. This method was used in [20] to obtain the consistent set of masses and decay constants used there.

Vacuum-expectation-values
The corrections to the vacuum expectation values (vevs) 0|qq|0 for up, down and strange quarks in the isospin limit were calculated at two-loops in three flavour ChPT in [15]. The expression for the up and down quark vev are identical since we are in the isospin limit.
We write the expressions in a form analoguous to the decay constant treatment: The superscripts indicate the order of the diagrams in p that each contribution comes from. The lowest order values are −F 2 0 B 0 . Note that the vevs are not directly measurable quantities. They depend on exactly the way the scalar densities are defined in QCD. ChPT can be used for them when a massindependent, chiral symmetry respecting subtraction scheme is used. M S in QCD satisfies this, but there are other possibilities. Even within a scheme, B 0 and the quark masses depend on the QCD subtraction scale μ QCD in such a way that B 0 m q is independent of it. The higher order corrections in this case also depend on the LECs for fully local counterterms, H r 1 , H r 2 at order p 4 and C r 91 , . . . , C r 94 at p 6 . When the scalar density is fully defined, measuring these quantities in e.g. lattice QCD and comparing with the ChPT expressions is a well defined procedure.
The functions are defined in massesdecayvev.h and implemented in massesdecayvev.cc.

Masses and decay constants at finite volume
The expressions treated in this section have been derived in [20]. A general remark is that care should be taken to set the precision in the loop integrals sufficiently high. For the one-loop integrals setting it very high is usually no problem. For the sunset integrals the evaluation can become very slow. It is strongly recommended to play around with the settings and compare the outputs for the two ways to evaluate the integral. The theta and Bessel function evaluation approach the correct answer differently. For most cases it is possible to have rsacc set smaller than racc.
For many applications it is useful to calculate the very time consuming parts, those labeled 6RV, once and store them. They only depend nontrivially on the masses and size of the finite volume. The decay constant dependence is very simple and there is dependence on the NLO LECs L r i . The results presented in this section are with periodic boundary conditions and an infinite extension in the time direction. They are also restricted to the case where the particle is at rest, i.e. p = 0.

Masses at finite volume
The finite volume corrections to the masses squared 6 are defined as the difference of the mass squared in finite volume and in infinite volume: These definitions are for a = π, K , η. These functions are available as mpi4Vt, mpi6Vt, mpi6LVt mpi6RVt respectively for the pion. mpi4Vt and mpi6RVt have as arguments a physmass and the length L. The other two have 6 Note that in other papers the corrections to the mass itself are sometimes quoted.
as arguments physmass,Li, double L. The final letter "t" indicates the evaluation using theta functions. With a "b" instead they use evaluation via Bessel functions. The equivalent functions for the kaon (pi to k) and eta (pi to eta) are also available. All these are defined in massdecayvevV.h and implemented in massdeca yvevV.h.

Decay constants at finite volume
The finite volume corrections to the decay constants are defined as the difference of the mass squared in finite volume and in infinite volume: These definitions are for a = π, K , η. Note that the correction is defined to the value of the decay constant, not divided by the lowest order decay constant as in (15). The functions are available as fpi4Vt, fpi6Vt, fpi6LVt fpi6RVt respectively for the pion. fpi4Vt and fpi6RVt have as arguments a physmass and the length L. The other two have as arguments physmass,Li, double L. The final letter "t" indicates the evaluation using theta functions. With a "b" instead they use evaluation via Bessel functions. The equivalent functions for the kaon (pi to k) and eta (pi to eta) are also available. All these are defined in massdecayvevV.h and implemented in massdeca yvevV.h.

Error handling
Error handling has been dealt with in a very simple manner. Most functions print out a message to standard output if something doesn't seem right. In particular, since the subtraction scale is present in several inputs, many functions check if these are the same and print out a message if not.
Errors due to a zero in a denominator are not caught and might lead to a crash.

Warnings
The definition of higher orders in ChPT is not unique. In this program collection, we have consistently chosen to rewrite all results in the physical masses and the physical pion decay constant, but note that even this is not fully unique. While there is usually a standard choice for the lowest order expression, at one-loop order this is often not the case since the Gell-Mann-Okubo relation can be used to rewrite the dependence on the η mass.
The files contain many internal functions as well as some extensions which are not described in this manuscript. These might change in future releases and have in general not been as fully tested as the described ones. Use at your own risk.

Possible extensions
Two-loop results are known for many more ChPT quantities also including isospin violation as well as for two-flavour ChPT and the partially quenched case. In addition at the oneloop level very many extensions exist for inclusion of the internal electro-magnetic interaction, finite volume effects, twisting and various extensions that include finite a effects in lattice gauge theory.
Another extension is how the higher orders are actually defined. In this program collection so far we have consistently chosen to rewrite all results in the physical masses and the physical pion decay constant. Implementations of other choices of higher orders, in particular in terms of lowest-order quantities are planned.
A last extension worth mentioning is the inclusion of the existing two-flavour ChPT results.

Conclusions
This paper describes a library of useful numerical programs in C++ for ChPT at upto two-loop order. Care has been taken to be independent of other libraries. In particular a number of numerical routines has been reimplemented in the numerical algorithm part of the library, libjbnumlib.a. The more ChPT direct functions like the loop integrals, a number of data structures to deal with LECs and the result for the masses, decay constants and decay constants are put in libchiron.a. Finite volume results are included for the masses and decay constants.
A simple Makefile as well as a large number of testing/example programs are included.