text:
# A Lossy Method for Compressing Raw CCD Images

## 0.1 Introduction

Optical and infrared instruments now routinely produce huge amounts of image data. The largest current common-user CCD mosaic has $`12\mathrm{k}\times 8\mathrm{k}`$ pixels (Veillet 1998); a single image from such a mosaic is 192 MB in size. The largest current common-user infrared detector mosaic has $`2\mathrm{k}\times 2\mathrm{k}`$ pixels (Beckett et al. 1998), but makes up for its smaller size by being read more frequently. More than a few nights of data from such instruments can easily overwhelm workstation-class computers. Compression is a solution to some of the problems generated by these large quantities of data. Similarly, compression can improve the effective bandwidth to remote observatories (in particular space observatories), remote data archives, and even local storage devices.

This paper describes in detail a lossy method for compressing raw images and presents a quantitative comparison to other lossy methods. It is organized as follows: § 2 reviews the limitation of lossless compression and the motivation for lossy compression; § 3 briefly summarizes the most relevant previous work on lossy compression; § 4 describes the new method; § 5 presents results on the distribution of differences between the original and compressed images; § 6 investigates the performance of the method with particular reference to hcomp; § 7 compares the performance of the method to other similar methods; § 8 discusses the suitability of the method for compressing raw data and investigates some consequences of such use; § 9 discusses how the method might be improved; and § 10 presents a brief summary.

## 0.2 Lossless Compression

### 0.2.1 The Shannon Limit

Compression methods can be classified as lossless or lossy depending on whether compression changes the values of the pixels. The difficulties of lossless compression of CCD images have been discussed by White (1992), Press (1992), and Véran & Wright (1994). The noise in images from an ideal CCD consists of Poisson noise from the signal and uncorrelated Gaussian read noise. If the signal is measured in electrons, the variance of the Poisson noise is equal to the signal and the variance of the read noise is typically 10–100. However, CCD images are often quantized to sample the read noise, with typical analog-to-digital converter gains being 1–5 electrons. This means that the larger noises associated with larger signal levels can be hugely over-sampled. If this is the case, the low-order bits in an image are essentially uniform white noise, contain no information, and cannot be compressed.

To illustrate this, series of $`512\times 512\times \text{16-bit}`$ FITS images of pure Gaussian noise were created. The images were written with with a BSCALE of $`q\sigma `$, so the standard deviation is sampled by a factor of $`1/q`$. Shannon’s first theorem (Shannon 1948ab, 1949) allows us to calculate the optimal compression ratio for these images. The theorem states that if a stream of bits is divided into fixed-length “input code units”, then the minimum number of bits required to encode each input code unit is just the Shannon entropy,

$$H\underset{p_i0}{}p_i\mathrm{log}_2p_i,$$
(1)
where $`p_i`$ is the normalized frequency of the input code unit $`i`$. From this we can derive that the optimal input code unit is the longest one for which the bits are still correlated; since in this case, each pixel is independent, the optimal input code unit is simply a 16-bit pixel. The resulting optimal compression ratio (the ratio of the compressed size to original size) will be $`H/16`$.

Figure 1 shows as a function of $`q`$ the optimal compression ratio determined by calculating the Shannon entropy explicitly. As can be seen, the optimal compression ratio is a constant minus $`(\mathrm{log}_2q)/16`$. Intuitively, this is expected: the number of incompressible low-order bits grows as $`\mathrm{log}_2q`$ and the fraction of the compressed image occupied by these bits grows as $`(\mathrm{log}_2q)/16`$. This result is not new; Romeo et al. (1999) investigated the Shannon limit for compression of a quantized Gaussian distribution, and their equation (3.3), in my notation, is

$$H\mathrm{log}_2\sqrt{2\pi e}\mathrm{log}_2q,$$
(2)
This equation was derived in the limit of small $`q`$, but Gaztañaga et al. (2001) have shown that the corrections for finite $`q`$ are small for $`q`$ at least as large as 1.5. Not surprisingly, the Shannon entropies and compression ratios calculated using this equation and calculated explicitly for the Gaussian noise images are in almost perfect agreement.

Real astronomical images contain information in addition to noise, but nevertheless this exercise demonstrates the futility of attempting to achieve good compression ratios for images that contain over-sampled noise. For example, if the noise is over-sampled by a factor of 60 (which corresponds to $`q=60`$), it is *impossible* to compress the image losslessly to better than half its size; too much of the image is occupied by white noise in the low-order bits. This exercise also explains why lossless methods are often more successful on images with low signals, such as biases or short exposures, and less successful on images with high signals, such as flats, deep exposures, and infrared images.

We are forced to the conclusion that if we wish to significantly and consistently compress raw images, including those with high background levels, there is no alternative but to use lossy compression methods.111The results of White & Becker (1998) appear to contradict this statement. These authors achieve compression ratios of 17%–32% on four WFPC2 images using lossless Rice compression. However, WFPC2 images have low backgrounds, due to the small pixels of WFPC2 and low sky brightness in orbit, and the gain normally under-samples the read noise. It is likely that none of their images has a large white noise component, and as such their sample of images is not representative of the full range of raw astronomical images.

### 0.2.2 Real Lossless Compression Methods

The Shannon limit is a theoretical limit and is independent of any particular real compression method. Nevertheless, it is useful to consider how well real lossless compression methods perform on Gaussian noise both to introduce these methods and to investigate how well they perform in this idealized case. The methods considered are Huffman coding, Arithmetic coding, hcomp, bzip2, gzip, and lzop.

Huffman coding (Huffman 1952; Press et al. 1992, § 20.4) and Arithmetic coding (Whitten, Neal, & Cleary 1987; Press et al. 1992, § 20.5) make use of the frequency distribution of fixed input code units. As such, if 16-bit pixel values are taken as the input code unit, the compression ratio should be close to the Shannon limit (provided the image is large compared to the frequency table). Indeed, Gaztañaga et al. (2001) show that Huffman coding of quantized Gaussian noise is very nearly optimal for $`q1`$ (see their Figure 1). For this reason, we need not investigate actual implementations of these methods but instead can approximate their compression ratios by the Shannon limit.

Hcomp (White 1992) is specialized to compressing 16-bit astronomical images and uses a complex algorithm based on a wavelet transformation of the pixel values, optional quantization of the coefficients, and quad-tree compression of the coefficients. Hcomp can be used losslessly by omitting the quantization of the coefficients.

Bzip2 (Seward 1998), gzip (Gaily 1993), and lzop (Oberhumer 1998) are methods for compressing general byte streams and use dictionary-based algorithms with variable-length input code units. They are portable, efficient, and free from patents. Sources and executables are freely available. They represent different trade-offs between compression ratio and speed, with bzip2 typically being slowest but achieving the best compression ratios, lzop typically being the fastest but achieving the worst compression ratios, and gzip being intermediate in both speed and compression ratio.

Figure 1 shows the compression ratio achieved by hcomp (used losslessly), bzip2, gzip, and lzop. The compression ratios have the same trend as the Shannon limit – they also cannot compress white noise in the low order bits – but the constant is larger. (Lzop is a slight exception in that it is initially unable to compress the image but then follows the expected trend.) Both bzip2 and hcomp are quite close to the Shannon limit. The surprisingly good performance of bzip2 suggests that it might be useful in compressing images; a large portion of the rest of this paper will be devoted to evaluating this suggestion.

## 0.3 Prior Work on Lossy Compression

Louys et al. (1999) have recently scrutinized lossy compression methods, including hcomp (White 1992). Their investigation concentrated on relatively high compression of reduced images. Compressing raw images is somewhat different, in that we wish to maintain all of the information in the image and avoid introducing spurious features. Even weak features can be important, as subsequent co-addition, spatial filtering, or other processing may increase their significance. Therefore, for raw images we require methods that suppress only those parts of an image that contains no information (the white noise in the low-order bits) and preserve perfectly those parts that do contain information (the mid-order and high-order bits). None of the lossy methods considered by Louys et al. (1999) can give this guarantee. For raw images, we must consider other methods.

More promising is the work of White and Greenfield (1999), Nieto-Santisteban et al. (1999), and Gaztañaga et al. (2001). These authors have all presented lossy compression methods that work by discarding white noise and then using a lossless compression method. White and Greenfield (1999) discard white noise from generic floating-point images by quantizing each row by a fraction of the estimated noise and converting to integers, Nieto-Santisteban et al. (1999) discard white noise from simulated NGST data by scaling and taking the square root, and Gaztañaga et al. (2001) discard white noise from simulated data from cosmic microwave background experiments by quantizing by a fraction of the standard deviation. White and Greenfield (1999) and Nieto-Santisteban et al. (1999) use the Rice compression algorithm (Rice, Yeh, & Miller 1993) and Gaztañaga et al. (2001) use Huffmann coding (Huffman 1952; see Press et al. 1992, § 20.4).

Comparisons to the work of White and Greenfield (1999) and Nieto-Santisteban et al. (1999) are made in § 0.7; the method suggested by Gaztañaga et al. (2001) is very similar to the one suggested here, but has a completely different context.

## 0.4 Method

The quantization compression method is so simple that it is somewhat embarrassing to describe it as a “method”. It consists of lossy quantization (or resampling) in brightness of the original image to produce a quantized image, followed by lossless compression of the quantized image. Decompression consists of reversing the lossless compression to recover the quantized image. The resampling quantum will normally be a specified fraction of the expected noise in the background. The aim of quantization is to reduce the over-sampling of the noise, reduce the number of bits of white noise, and permit the lossless compression method to perform better. The quantization method is intended to be used on raw data from CCD-like devices.

### 0.4.1 Quantization

The first stage of the quantization method is to produce a quantized image. For definiteness, the quantization is described in terms of integer FITS images (Wells, Greisen, & Harten 1981). This is the most common format for archived raw data. The extension to other image formats is obvious.

Each integer $`n`$ in the original FITS image is related to a data value $`b`$ by

$$b=b_{\mathrm{zero}}+nb_{\mathrm{scale}},$$
(3)
where $`b_{\mathrm{zero}}`$ and $`b_{\mathrm{scale}}`$ are taken from the BZERO and BSCALE records in the FITS header and have defaults of 0 and 1. The data quantum $`Q_b`$ is chose to be an integral multiple $`Q_n`$ of $`b_{\mathrm{scale}}`$. The quantized brightness $`\widehat{b}`$ is thus

$$\widehat{b}=b_{\mathrm{zero}}+Q_b[b/Q_b]b_{\mathrm{scale}}.$$
(4)
Here, $`[x]`$ is the nearest integer to $`x`$ with half integers rounded to the nearest even integer to help avoid biases. The quantized image is written with the same BZERO and BSCALE values as the original image. Thus, the quantized integer $`\widehat{n}`$ in the quantized FITS file is related to‘ the quantized data value $`\widehat{b}`$ by

$$\widehat{n}=(\widehat{b}b_{\mathrm{zero}})/b_{\mathrm{scale}}.$$
(5)
It can be seen that this quantization can be carried out entirely in integer arithmetic as

$$\widehat{n}=Q_n[n/Q_n].$$
(6)
Certain values in the original image are special and are treated accordingly. Thus, “blank” pixels (those which have the value specified by the BLANK record in the FITS header) are not changed and any pixel whose quantized value would be the blank value is not changed. The values corresponding to the maximum and minimum integer values are treated similarly, as these values are often used to flag saturated data or implicitly blank pixels. (Some images have ancillary information, such as the exposure time and engineering information, encoded in “data” pixels; this information should be transferred to the header prior to compression.)

The quantum $`Q_b`$ should be determined according to the noise model for the detector. For CCD-like devices, the estimated standard deviation $`\sigma _b`$ of the noise is given by

$$\sigma _b=\{\begin{array}{cc}((bb_{\mathrm{bias}})g+r^2)^{1/2}/g\hfill & b>b_{\mathrm{bias}}\hfill \\ r/g\hfill & bb_{\mathrm{bias}}\hfill \end{array}.$$
(7)
Here $`g`$ is the gain in electrons, $`r`$ is the read noise in electrons, and $`b_{\mathrm{bias}}`$ is the bias level. We first determine the bias level $`b_{\mathrm{bias}}`$ and data background level $`b_{\mathrm{data}}`$ in the image, perhaps from the mean, median, or mode in certain regions, and hence the noise by applying Equation (7). We then determine the quantum as

$$Q_b=\{\begin{array}{cc}b_{\mathrm{scale}}q\sigma _b/b_{\mathrm{scale}}\hfill & q\sigma _b>b_{\mathrm{scale}}\hfill \\ b_{\mathrm{scale}}\hfill & q\sigma _bb_{\mathrm{scale}}\hfill \end{array},$$
(8)
or

$$Q_n=\{\begin{array}{cc}q\sigma _b/b_{\mathrm{scale}}\hfill & q\sigma _b>b_{\mathrm{scale}}\hfill \\ 1\hfill & q\sigma _bb_{\mathrm{scale}}\hfill \end{array}.$$
(9)
Here $`x`$ is the largest integer no greater than $`x`$. Thus, the parameter $`q`$ controls the coarseness of the quantization by setting the size of the quantum $`Q_b`$ to be the smallest integer multiple of $`b_{\mathrm{scale}}`$ no larger than $`q`$ times the noise. Suitable values for $`q`$ are 0.5 to 2 (see § 0.5 and § 0.6).

There is no requirement that the quanta be equal over the entire image. The data and overscan regions of an image should probably be quantized with different quanta as they likely have different characteristic noises. Furthermore, if the background and hence the noise vary significantly within the data region, it makes sense to subdivide the data region and use different quanta within each subregion. However, it is important the quanta are constant within regions and are not determined on a pixel-to-pixel basis, as this can introduce small but nevertheless real biases (see § 0.7).

### 0.4.2 Compression

The quantized FITS file is no smaller than the original FITS file. Therefore, the second stage of the quantization method is to compress the quantized FITS file using a lossless compression method. The methods investigated here are bzip2 (Seward 1998), gzip (Gaily 1993), lzop (Oberhumer 1998), and hcomp (White 1992). Hypothetical optimal 16-bit and 32-bit encoders (which will approximate Huffman or Arithmetic encoders) are also considered.

### 0.4.3 Decompression

The quantized image is recovered as a standard FITS file by reversing the lossless compression using bunzip2, gunzip, lzop, or hdecomp.

## 0.5 Distribution of Differences

One advantage of the simplicity of the quantization method is that we can derive useful information on the distribution of the differences between the quantized and original images. These give the confidence required to apply the quantization method to raw data. Note that the distribution of differences is determined only by the quantization stage; the compression stage is lossless.

We define the normalized difference $`\mathrm{\Delta }`$ between the quantized and original data values to be

$$\mathrm{\Delta }=\frac{\widehat{b}b}{\sigma _b}$$
(10)
A firm upper bound on the magnitude of $`\mathrm{\Delta }`$ is

$`|\mathrm{\Delta }|{\displaystyle \frac{Q_b}{2\sigma _b}}0.5q.`$ (11)
Thus, no difference is larger in magnitude than $`q/2`$ standard deviations in the background. Provided $`q`$ is not larger than $`2`$, we expect $`\mathrm{\Delta }`$ to be roughly uniformly distributed about zero (on discrete points corresponding to integer multiples of $`b_{\mathrm{scale}}/\sigma _b`$) and so we expect the mean and standard deviation of $`\mathrm{\Delta }`$ (the mean and RMS normalized difference) to be approximately

$`\mu _\mathrm{\Delta }0`$ (12)
and

$`\sigma _\mathrm{\Delta }{\displaystyle \frac{q}{\sqrt{12}}}0.29q.`$ (13)
The quantity $`\sigma _\mathrm{\Delta }`$ corresponds to the square-root of the “digital distortion” $`𝒟`$ discussed by Gaztañaga et al. (2001). They show that a quantization of a distribution of standard deviation $`\sigma _b`$ will result in a distribution with a standard deviation $`\sigma _{\widehat{b}}`$ given by

$$\frac{\sigma _{\widehat{b}}}{\sigma _b}\sqrt{1+\frac{q^2}{12}}1+0.144q,$$

which comes from adding the variances in quadrature and assuming that the difference $`\widehat{b}b`$ is uncorrelated with $`b`$ (which is a good approximation provided $`q`$ is small).

## 0.6 Performance

This section presents tests and comparisons of the performance of the quantization compression method proposed here, qcomp, in terms of the distribution of differences, compression ratio, and speed of compression and decompression. Hcomp and optimal 16-bit and 32-bit encoders serve as points of reference.

### 0.6.1 Trial Images

Five trial images were selected to sample different noise regimes and scenes. The images were obtained with a CCD at the Observatorio Astronómico Nacional at Sierra San Pedro Mártir and saved as $`512\times 512\times 16\text{-bit}`$ FITS files with a BSCALE of 1. The images were (a) a bias exposure, (b) a shallow exposure of a sparse field (a standard star), (c) a deep exposure of a relatively sparse field (a small galaxy), (d) a deep exposure of a complex field (a nebula), and (e) a flat-field exposure. To make the comparison with hcomp fair, the overscan sections were removed, as hcomp does not distinguish between the data and overscan sections.

The bias and background levels in each image were estimated from the mode in the exposed and overscan regions, and the noise $`\sigma _b`$ estimated by applying Equation 7. The standard deviation $`\widehat{\sigma }_b`$ of the pixel values was also determined (with $`5\sigma `$ clipping). Figure 2 shows the images and the histograms of the pixel values, along with the expected contribution contribution to the histogram width from the noise. The histograms of the images (a), (b), and (c) are dominated by noise. However, the histogram of image (d) is much wider than the noise because of structure in the scene and the histogram of image (e) is much wider because of structure in the flat field (large-scale structure of around 1.5% and pixel-to-pixel variations of around 0.5%).

Any choice of trial images is likely to be somewhat arbitrary. However, one advantage of these images is that they cover a good range of noise regimes and scenes. This implies that if a compression method can produce good results for all of the trial images, it is likely to be produce good results for other sets of images, regardless of the relative frequency of each type of image. Thus, we see that a successful compression method must achieve not only a good mean compression ratio for the whole set of trial images but also consistently good compression ratios for each individual image.

### 0.6.2 Compression

The images were compressed using the following methods:

1. Lossless compression by hcomp with the “-s 0” option.
2. Lossless compression using hypothetical optimal 16-bit and 32-bit methods. The compression ratios for these methods were determined by calculating the Shannon entropy for each image, given by Equation 1. Although these methods are hypothetical, they will approximate Huffman and Arithmetic coding.
3. Lossless compression by direct application of bzip2, gzip, and lzop with options varying from “-1” (fastest/worst) to “-9” (slowest/best).
4. Lossy compression using hcomp, with the degree of loss being controlled by value of the “-s” option. The value of this option was $`2q\sigma _b`$, which gives roughly the same RMS difference as the quantization method with the same value of $`q`$.
5. Lossy compression using the quantization method described in § 0.4 implemented by the program qcomp (available from the author), with hypothetical optimal 16-bit and 32-bit compression, bzip2, gzip, lzop, and lossless hcomp used for the actual compression stage.

The distributions of differences, compression ratios, and speeds are shown in Tables 2 and 3 and Figures 3 and 4. The lossy methods are shown for $`q`$ values of 0.5, 1, 2, 4, 8, and 16. Table 2 and Figure 3 show the RMS, maximum absolute, and mean normalized difference $`\mathrm{\Delta }`$ between the compressed and original images. (The normalization is to $`\sigma _b`$, the noise in the image; see Equation 10.) Obviously, the lossless methods do not appear as their differences are uniformly zero. Table 3 shows the compression ratio achieved for each image (the ratio of the compressed size to original size), the mean compression ratio for all five images, and the speed at which the five images were compressed and decompressed. These speeds were measured on a computer with a single 1.9 GHz Pentium IV processor. Speeds are not given for the optimal 16-bit and 32-bit methods, as these methods are hypothetical.

The compression speeds quoted in Table 3 do not include the time taken to determine the background level and calculate the noise. Determining these using the mean, median, or mode proceeds at about 11 Mpixel/s and reduces the compression speeds to 0.79 Mpixel/s for qcomp/bzip2 (9% slower), 2.84 Mpixel/s for qcomp/gzip (26% slower), 4.03 Mpixel/s for qcomp/lzop (37% slower), 1.47 Mpixel/s for qcomp/hcomp (13% slower), and 1.94 Mpixel/s for hcomp (18% slower) for $`q=1`$.

Figure 4 shows $`q`$, the RMS difference, and maximum absolute differences as functions of the compression ratio for the more effective lossy compression methods.

### 0.6.3 Lossless Methods

The upper section of Table 3 gives the compression ratios for the lossless methods. The best methods are hcomp, bzip2, and the optimal 16-bit and 32-bit methods, which achieve mean compression ratios of about 0.4.

In contrast to the compression of pure Gaussian noise discussed in § 0.2, hcomp and bzip2 achieve better mean compression ratios than the optimal 16-bit method. How can this be? The explanation is that in the deep and flat-field images, (c), (d), and (e), there are correlations between adjacent pixels that both hcomp and bzip2 can exploit – hcomp because its wavelets can extend over more than one pixel and bzip2 because it can increase the size of its input code unit to cover more than one pixel. This also explains why the optimal 32-bit method, which considers pixel in pairs, performs better than the optimal 16-bit method, which considers pixels individually. However, in the bias and shallow images, (a) and (b), the noise component is dominant, and the optimal 16-bit method performs slightly better than hcomp and bzip2, just as for pure Gaussian noise. The optimal 32-bit method performs slightly better than the optimal 16-bit method even for the bias image, which suggests that there are small pixel-to-pixel correlations even in this image, presumably due to the bias structure. (Another way of looking at this is that optimal 16-bit compression uses only the histogram of pixel values whereas the other methods use both the histogram and, to varying degrees, the image; the image adds information on correlations between pixels, and allows for more efficient compression.)

Nevertheless, the mean compression ratio of about 0.4 masks variations of factors of 2–3 in the compression ratio for individual images, with the bias and shallow images, (a) and (b), being compressed well to about 0.25 and the deep and flat-field images, (c), (d), and (e), being compressed poorly to 0.4–0.7. This is as expected; the bias and shallow images compress well as their noises are only moderately over-sampled, but the deep and flat-field images compress badly as their noises are significantly over-sampled. In short, the lossless methods produce inconsistent results, and while they will compress some sets of data very well (for example, those dominated by observations of standard stars or images in low-background narrow-band filters), they will compress others very poorly (for example, those dominated by deep images with high backgrounds). In all, these results confirm what we already knew: if we want to achieve consistently good compression ratios, we have no choice but to adopt lossy methods.

(Table 3 also shows the effect of the options to bzip, gzip, and lzop. For bzip2 and gzip there is not much improvement in the compression ratio between the “-1” and “-9” options, but the compression speed drops noticeably. For lzop there is a significant improvement, but the compression speed drops precipitously. It seems clear that the “-1” option is preferable, and for this reason it was used in the lossless compression stage by qcomp.)

### 0.6.4 Lossy Methods

Table 2 and Figure 3 shows that the distributions of differences for qcomp are behaving much as predicted in § 0.5. For $`q2`$, the RMS difference is roughly $`0.3q\sigma _b`$, the maximum absolute difference is $`0.5q\sigma _b`$, and the mean difference is roughly zero. For larger values of $`q`$, the mean difference ceases to be close to zero, as might be expected when the quantization is too coarse to adequately sample both sides of a sharply-peaked distribution. This suggests that qcomp should not be used beyond $`q2`$.

Comparison of the distributions of differences for qcomp and hcomp in Figure 3 reveals two trends. First, for $`q2`$, the RMS differences are very similar. This is not surprising as the hcomp option was set to give this result. However, qcomp has lower values of the maximum absolute difference and mean difference. For larger values of $`q`$, the RMS differences and mean differences for hcomp are much better than those for qcomp.

Moving from the performance of the quantization step to the performance of the compression step, we first investigate what appears to be a worrying anomaly: when the bias image, (a), is quantized with $`q=1`$, Table 3 shows that optimal 16-bit compression gives a compression ratio of 0.163, whereas Equation 2 indicates that optimal 16-bit compression of quantized Gaussian noise with $`q=1`$ should give a compression ratio of about 0.128. The bias image has very little structure and so should to a good approximation be simply Gaussian noise, so why are the compression ratios so different? Only a small part of the difference can be explained by corrections to Equation 2 for finite $`q`$; numeric experiments give a true compression ratio for Gaussian noise with $`q=1`$ of 0.134. The explanation for the remaining difference is that the noise $`\sigma _b`$ in the bias image is 2.5 and is not an integer multiple of $`b_{\mathrm{scale}}`$, and so the quantum $`Q_b=b_{\mathrm{scale}}q\sigma _b/b_{\mathrm{scale}}`$ is rounded down to 2. Thus, the proper comparison to compressing the bias image with $`q=1`$ is compressing Gaussian noise with $`q=0.8`$. By Equation 2, this would have a compression ratio of 0.148, and numerical experiments give a value of 0.153, which is much closer to the 0.163 of the bias image. The remaining difference is presumably due to structure in the bias and non-Gaussian contributions to the noise (for example, the slightly flat top to the histogram in Figure 2a). Similarly, the shallow image (b) should be properly compared with Gaussian noise quantized with $`q=0.7`$, which has a compression ratio of 0.162. The other images have large values of $`\sigma _b`$, so the truncation of $`Q_b`$ does not effect them significantly.

Figure 4a shows the mean compression ratios achieved by qcomp/optimal-16, qcomp/optimal-32, qcomp/bzip2, and qcomp/gzip as functions of $`q`$. (Neither qcomp/lzop nor qcomp/hcomp are shown, as they are significantly worse.) The best mean compression ratios are achieved by qcomp/bzip2; again, bzip2 appears to be finding correlations between pixels that are not exploited by either the optimal 16-bit or 32-bit methods, especially in the cases of the deep image of the nebula (d) and the flat-field image (e). For $`q=1`$ and $`q=2`$ the mean compression ratios are 0.19 and 0.13, that is, factors of two and three better than lossless compression with bzip2. Furthermore, Table 3 shows that qcomp/bzip2 is achieving consistently good compression ratios which vary by only 20% from 0.176 to 0.214 for $`q=1`$.

Appropriate figures of merit for comparing hcomp and qcomp/bzip2, are the RMS and maximum absolute difference at a given mean compression ratio, and these are shown in Figure 4b and c. For compression ratios between 0.13 and 0.25 (which correspond to $`q`$ between 0.5 and 2), qcomp/bzip2 is better than hcomp both in RMS difference and in maximum absolute difference. For compression ratios smaller than 0.1 (which correspond to $`q>2`$), qcomp is significantly worse than hcomp, first in the RMS difference but eventually also in the maximum absolute difference. This is not surprising, as we saw above that the quantization becomes too coarse for $`q>2`$. Thus, in terms of the distribution of differences at a given compression ratio, of the methods we have considered, qcomp/bzip2 is best for compression ratios of 0.13 and above (which corresponds to $`q2`$) and hcomp is best for compression ratios smaller than 0.1 (which correspond to $`q4`$). The one advantage that hcomp has over qcomp/bzip2 for $`q2`$ is that it compresses about three times faster, although both decompress at about the same speed.

At the same value of $`q`$ (i.e., the same distribution of differences), qcomp/gzip2 produces significantly worse compression ratios then qcomp/bzip2. For example, at $`q=1`$ the mean compression ratio of qcomp/gzip is 0.31 compared to 0.19 for qcomp/bzip2. This appears to be not much better than the 0.41 achieved by lossless compression by bzip2, but qcomp/gzip has two advantages: it is roughly four times faster then either lossless bzip2 or lossy qcomp/bzip2, both at compression and decompression, and it achieves more consistent compression ratios than lossless compression.

## 0.7 Comparison to Other Quantization Methods

This section compares the quantization method proposed here, qcomp, to the similar quantization methods proposed by White and Greenfield (1999) and Nieto-Santisteban et al. (1999).

### 0.7.1 Distribution of Differences

White and Greenfield (1999) calculate the quantum from the empirical standard deviation in the image rather than using a noise model. This can cause the actual noise in raw data to be significantly over-estimated. Consider the flat field image (e), which has standard deviation of more than 500 (after 3 sigma clipping), but a noise of only about 75. Even if we eliminate the contribution from large-scale structure, the pixel-to-pixel standard deviation is about 160. Obviously, as the empirical standard deviations is larger than the noise, quantizing on the basis of the empirical standard deviation would remove real information on the variations in the flat field. Similar problems can occur in any images with high backgrounds or more generally, in which the empirical standard deviation is not dominated by noise. This comparison is somewhat unfair, in that White and Greenfield (1999) do not claim that their method is suitable for raw data, but it does illustrate the danger of assuming that the noise in an image is related to the empirical standard deviation.

Nieto-Santisteban et al. (1999) quantize by taking the square root of scaled data. This is equivalent to calculating the quantum on a pixel-by-pixel basis and assuming that the read noise is zero. To investigate this, a set of images were created with mean values from 1 to 200 electrons, a gain of 1 electron, no read noise, and with appropriate Gaussian noise (although strictly speaking the noise should be Poissonian), and written with BZERO of 0 and BSCALE of 1. The images were quantized two different ways: first using a single quantum for each image calculated by applying Equations 7 and 8 with $`q=1`$ to the median pixel value (the quantization method proposed here) and second using individual quanta for each pixel calculated by applying Equations 7 and 8 to each pixel value (equivalent to the quantization method suggested by Nieto-Santisteban et al.). The mean differences between the quantized and compressed images are shown in Figure 5. Quantizing on a pixel-to-pixel basis results in much larger mean differences than using a single quantum, and the differences have an oscillatory structure with a change of gradient every whole integer squared. This probably results from the quantum changing with pixel value, so that the two sides of the distribution of pixel values are quantized slightly differently. The lower envelope appears to be roughly $`0.15/\sigma _b`$. The mean difference is small even in the pixel-to-pixel case (of order 1% of the noise for typical signals), but nevertheless quantizing with the same quantum gives a much smaller mean difference (of order 0.05% of the noise).

Thus, it appears that the quantization mehtod proposed here has advantages over the other quantization methods when we consider the fidelity of the compressed image to the original image.

### 0.7.2 Compression Ratio and Speed

Rice, Yeh, & Miller (1993) describe a special-purpose compression algorithm designed to compress low-entropy 16-bit integer data. This method is widely referred to as “Rice compression” and has been used by White and Greenfield (1999), Nieto-Santisteban et al. (1999), and, in a modified form, White & Becker (1998).

The compression ratio achieved by Rice compression appears to be only very slightly better than that achieved by gzip; White & Becker (1998) report 0.255 for Rice and 0.272 for gzip while Nieto-Santisteban et al. (1999) report 0.193 for Rice and 0.202 for gzip (for 2 bits of noise).

White & Becker (1998) report that gzip is “far slower” than their implementation of Rice compression and Nieto-Santisteban et al. (1999) report that gzip is 20 times slower than White’s implementation of Rice compression. Both of these articles focus on compressing data on a satellite with only limited computing resources, and so understandably do not consider the speed of decompression. Neither mention the flag given to gzip, but the default is equivalent to “-6” and is about half as fast at compression as “-1”. If we assume that the default was used, then Rice compression appears to be about 10 times faster than gzip with the “-1” flag.

It would appear then that the methods proposed by White and Greenfield (1999) and Nieto-Santisteban et al. (1999) would produce similar compression ratios to qcomp/gzip2 but would would compress much more quickly; neither method is likely to produce compression ratios that are as good as those produced by the much slower qcomp/bzip2. Of course, it would be possibly to adapt qcomp to use Rice compression or the other methods to use bzip2, in which case they should produce very similar compression ratios and speeds.

## 0.8 Discussion

### 0.8.1 Suitability for Compressing Raw Data

The quantization compression method described here, qcomp, was designed to compress raw data. With the results of the tests and comparisons in § 0.6 and § 7, we can now critically evaluate its suitability for this task. We will see that that qcomp is better suited for compressing raw data than either hcomp (White 1992) or the quantization methods proposed by White and Greenfield (1999) and Nieto-Santisteban et al. (1999).

The fundamental advantage of qcomp over hcomp is that its straightforward nature allows firm limits to be placed on the distribution of differences (see § 0.5). So, for example, an image can be compressed with the knowledge that no pixel will change by more than a specified amount. Furthermore, the distribution of the errors is predictable. Hcomp does give such firm guarantees on the distribution of differences.

However, this advantage would be moot if qcomp were inferior in compression ratio and speed. The most interesting range of $`q`$ for compressing raw images is 0.5–2, that is, quantizations that change the values in the image by at most 0.25–1 standard deviations. In this range, qcomp/bzip2 gives mean compression ratios of 0.27–0.13 and gives roughly the same RMS difference as hcomp for similar compression ratios but with smaller maximum absolute differences and mean differences (see Figure 4). Furthermore, its decompression speed is very similar to hcomp. Its only disadvantage is that its compression speed is only one third that of hcomp, although it is still adequately fast; its compression rate of roughly 800 kpixel/s on a 1.9 GHz Pentium IV is much faster than the typical read rate of 40 kpixel/s for single-output CCDs and 160 kpixel/s for quad-output CCDs.

It is also worth considering qcomp/gzip, which achieves lower compression ratios of 0.38–0.23 for $`q`$ in 0.5–2, while compressing slightly faster and decompressing very much faster than hcomp. It also has the great advantage of portability; gzip is installed almost universally on workstations.

The current implementation of qcomp has a further small advantage over the current implementation of hcomp in it distinguishes between the data and overscan sections of an image and quantizes each appropriately; the current implementation of hcomp treats both identically, and so the two have to be separated and compressed individually to avoid either drastic over-compression of the overscan section or under-compression of the data section.

The advantages of the qcomp over the other quantization methods are that the manner in which White and Greenfield (1999) determine the quantization makes their method unsuitable for use with raw data (and indeed such a claim was never made) and the manner in which Nieto-Santisteban et al. (1999) quantize introduces a small bias that is not present with qcomp. When using the same lossless compression method, all of the quantization methods should be similar in speed. The one disadvantage of qcomp is that is requires that the noise in the background be roughly constant within fixed regions, whereas the method of Nieto-Santisteban et al. (1999) estimates the noise on a pixel-by-pixel basis and makes no such requirement. However, this may not be too restrictive as the noise varies only as the square root of the background (and more slowly if read noise is significant), so changes in the background produce much smaller changes in the noise.

Obviously, confidence in the use of qcomp to compress raw data would be significantly improved if end-to-end tests on compressed data gave results that were statistically indistinguishable from those derived from the original data. Such tests should include astrometry of stars, aperture and PSF-fitting photometry of stars, and photometry of low surface brightness sources. These tests are especially important for applications which have stringent requirements on the absence of biases, such as far infrared and sub-millimeter imaging and ultra low surface brightness studies at other wavelengths. Such tests will be presented in Paper II (Watson, in preparation).

### 0.8.2 Distribution

An important advantage of compressing files with qcomp/gzip is that they can be distributed with the expectation that no additional software will be needed to decompress them; gzip is installed almost universally on workstations. As bzip2 becomes more widespread (it already forms part of several Linux distributions), it too will gain this advantage. This is a major advantage over software that uses special-purpose compression methods.

### 0.8.3 Media Capacities

If we take a compression ratio of 0.2 as typical for qcomp/bzip2 with $`q=1`$, then media have effective capacities that are 5 times larger than their raw capacities. Thus, a 650 MB CD-ROM has an effective capacity of 3 GB, a 12 GB DAT DDS-3 tape has an effective capacity of 60 GB, and an 80 GB disk has a capacity of 400 GB. This is roughly equivalent to a generation of technology; in terms of capacity, a CD-ROM is almost a DVD-ROM, a DAT tape is effectively a DLT tape, and a single disk is effectively a large RAID array. Thus, compression can allow one to work with large data sets without acquiring the devices often considered essential.

### 0.8.4 Bandwidth

If images can be compressed and decompressed sufficiently quickly, compression increases the effective bandwidth of a device or transport mechanism by a factor equal to the inverse of the compression ratio. Thus, it may be faster to store images in their compressed form, even though there is an overhead in compressing or decompressing them.

Table 4 shows for several devices or transport mechanisms the raw bandwidths and the current and maximum effective bandwidths for compressed data, assuming a compression ratio of 0.3 appropriate for compression with qcomp/gzip with $`q=1`$. The maximum bandwidths assume an arbitrarily fast processor, so the overheads of compression and decompression drop to zero. The current read and write bandwidths were measured or estimated for qcomp/gzip on an early-2002 computer with single 1.9 GHz Pentium IV processor, which can compress at about 3 Mpixel/s and decompress at about 12 Mpixel/s. (A similar table could be constructed for qcomp/bzip2; the maximum effective bandwidths would be higher and the current effective bandwidths would be lower, except for the 56 kbps link which would have a current effective bandwidth of 0.015 Mpixel/s.)

Table 4 shows that even now compression results in an improvement in effective bandwidth for CD-ROMs, tapes, single disks, and all but the very fastest network connections. As processors become faster, the actual bandwidths will approach the maximum bandwidths, and, for fast enough processors, all devices will benefit. For example, achieving 83 Mpixel/s read bandwidth from a fast RAID array will require processors only 7 times faster than a 1.9 GHz Pentium IV; if processors continue to double in speed every 18 months, according to Moore’s law, such processors should be available in about 2006. Achieving 83 Mpixel/s write speed will require processors 28 times faster which should be available in about 2009. Alternatively, if we allow 8-way parallelism (see § 0.9.2), such speeds should be possible now for reading and in about 2005 for writing. (Again assuming 8-way parallelism, qcomp/bzip2 should give 125 Mpixel/s read bandwidth around 2006 and 125 Mpixel/s write bandwidth around 2008.)

## 0.9 Possible Improvements

### 0.9.1 Compression Ratio

The compression ratio could be improved by using a compression method that performed better than bzip2. However, the experiments above indicate that bzip2 is better than both optimal 16-bit and optimal 32-bit compression, so this may be difficult.

Furthermore, it is easy to show that bzip2 is performing very close to optimally. The mean compression ratio achieved by qcomp/bzip2 with $`q=1`$ is 0.193. The mean compression ratios achieved by bzip2 and optimal 16-bit compression of Gaussian noise with the same $`q`$ values (which, as mentioned above, are really $`q=0.8`$ and $`q=0.7`$ for images (a) and (b)) are 0.172 and 0.143. This suggests that a hypothetical compression method that compressed the information as well as bzip2 and compressed the noise as well as optimal 16-bit compression would achieve a compression ratio of about 0.164, an improvement of only 18%. Even a hypothetical compression method that compressed the information to nothing and compressed the noise as well as optimal 16-bit compression would achieve a compression ratio of about 0.143, an improvement of only 35%. It seems that bzip2 is close to optimal, and any further improvements will be small; there are no more factors of two to be gained.

One means to improve the compression ratio in high-background images would be to roughly flatten and/or subtract the sky prior to compression. This should narrow the histogram of pixel values and thereby improve subsequent compression. However, again, there’s not much to be gained. The difference between the compression ratios for the bias image (a) and the flat-field image (e) is only 10% for qcomp/bzip2 with $`q=1`$; bzip2 is doing a very good job of encoding both the pixel-to-pixel and large-scale structure in the flat-field image, and so would presumably perform similarly well on other high-background images.

### 0.9.2 Speed

An obvious improvement to the writing speed would be to parallelize the quantization and compression steps. This has already been implemented in qcomp, as the quantized image is piped to the lossless compressor as it is created. The improvements are up to 50% for qcomp/gzip and 10% for qcomp/bzip2.

To obtain further improvements, compression and decompression must be made faster. For both gzip and bunzip2, compression can be trivially parallelized by splitting the input between processes and concatenating the output from each process. Unfortunately, the file formats for gzip and bzip2 do not permit decompression to be parallelized. (The only way to determine the length of a compressed block is decompress it.) However, it would be a simple matter to modify the format of the compressed files to permit decompression to be parallelized (by prefixing each compressed block with its length); this would not be so conveniently portable, but might be worthwhile for private use.

Rice compression might be attractive if more speed is required and compression ratios similar to those of gzip are adequate. That said, even if the cost of lossless compression could be reduced to zero, determining the noise level and quantizing would still require effort and would limit the compression speed to only about 3 times that of qcomp/gzip. (It is possible that the noise level could be determined more rapidly by determining the noise from only a subsample of pixels, in which case the compression speed might improve by a somewhat larger factor.)

### 0.9.3 Robustness against Varying Backgrounds

One problem with the current implementation of qcomp is that it assumes that the whole data section of the image can be adequately characterized by a single background and hence a single background noise. (It does, however, allow the data and bias sections to have different backgrounds.) This will present a problem for images that have very different background levels, for example, spectra in the thermal infrared in which the thermal background increases rapidly with wavelength.

This situation could be improved by subdividing the data section into regions such that each region is adequately characterized by a single background. This could be achieved manually or perhaps adaptively. One adaptive scheme would be to subdivide regions recursively until the backgrounds in the subregions are not significantly different from the background in the region.

Another option would be to estimate the noise and hence the quantum on a pixel-by-pixel basis, following Nieto-Santisteban et al. (1999), but as has been shown this introduces a small bias in the mean difference between the quantized and original images.

## 0.10 Summary

A lossy method for compressing raw images has been presented. The method is very simple; it consists of lossy quantization (resampling in brightness) followed by lossless compression. The degree of quantization is chosen to eliminate the low-order bits that over-sample the noise, contain no information, and are difficult or impossible to compress.

The method is lossy but gives certain guarantees about the distribution of differences between the compressed and original images. In particular, it gives guarantees on the maximum absolute value, the expected mean value, and the expected RMS value of the difference. These guarantees make it suitable for use on raw data.

The method consistently reduces images to 1/5 of their original size while changing no value by more than 1/2 of a standard deviation in the background. This is a dramatic improvement on the compression ratios achieved by lossless compression. The method is adequately fast on current computers and would be relatively simple to parallelize.

A key feature of the method is that data can be uncompressed using tools that are widely available on modern workstations, which means that one can distribute compressed data and expect that it can be used without the need to install specialized software. This is achieved by writing the quantized image as a normal FITS file and compressing it with gzip and bzip2, which widely available general-purpose compression tools. It appears that bzip2 is compressing the data within a few tens of percent of optimally.

The next step in the development of this method is real-world testing with compressed raw data to ensure that the method does not degrade the results of astronomical analyses. Such end-to-end tests will be presented in Paper II (Watson, in preparation).

###### Acknowledgements.

I thank Enrique Gaztañaga his referee’s report, possibly the most useful I have ever received, and in particular for his suggestion to compare my results to the Shannon entropy limit. I thank Rick White for useful comments on hcomp and for blazing the trail. I thank Julian Seward, Jean-Loup Gailly, Mark Adler, and Markus Oberhumer for doing the hard work of writing bzip2, gzip, and lzop. I thank the staff of the Observatorio Astronómico Nacional de México at San Pedro Mártir for their warm hospitality in February 1999, when the ideas presented here were conceived.


format: text

text:
# Low Energy Pion–Nucleon Scattering in the Heavy Baryon and Infrared Schemes

## I Introduction

Pion-nucleon scattering is a fundamental process that one would like to describe using the low-energy realization of quantum chromodynamics, namely chiral perturbation theory . This is an attractive approach because it not only embodies chiral symmetry, which is fundamental to low-energy physics, but also offers a systematic expansion in powers of the momentum. Further it ensures unitarity order by order. Gasser and Leutwyler  have shown that chiral perturbation theory works nicely for mesons and high accuracy can now be achieved for the $`\pi \pi `$ scattering lengths .

However, the power counting fails when baryons are introduced . The power counting can be restored in heavy baryon chiral perturbation theory (henceforth referred to as HB)  where the heavy components of the baryon fields are integrated out. An alternative – the infrared regularization scheme (henceforth referred to as IR) – has been proposed , based on the ideas in Ref. . This preserves the chiral power counting and has the advantage that it is manifestly Lorentz invariant and avoids the voluminous effective Lagrangian of the HB approach. Generally favorable results have been obtained with the IR scheme in a number of applications . By suitable approximation of the IR expressions the HB formulae can be regained. A recent review of the HB and IR formalisms has been given by Meißner .

The first fit to the pion-nucleon scattering phase shift data using chiral perturbation theory was carried out in Ref. , henceforth referred to as I, using the HB scheme. In addition to the nucleon and pion fields, the $`\mathrm{\Delta }`$ resonance field was included explicitly since the intent was to fit out to energies in the $`\mathrm{\Delta }`$ resonance region. The calculation was carried to $`𝒪(Q^3)`$, where $`Q`$ is a generic small momentum scale, and it was found possible to obtain a reasonable fit up to energies slightly below the $`\mathrm{\Delta }`$ resonance. Subsequently Fettes and Meißner carried out several HB studies both with  and without  explicit inclusion of the $`\mathrm{\Delta }`$ field. The calculation with the $`\mathrm{\Delta }`$ field to $`𝒪(Q^3)`$ appeared to be a little better than the $`𝒪(Q^4)`$ calculation without it; for further discussion see Ref. . Fettes and Meißner  have also studied isospin violation in the $`\pi N`$ system, although here we shall focus on the isospin symmetric case. All of these calculations were carried out with the HB approach and it is natural to examine the IR method in this context. Becher and Leutwyler  have used the IR approach in the sub-threshold region, where the HB scheme is inappropriate, and found that the IR representation of the scattering amplitude was not sufficiently accurate to allow the extrapolation of the experimental data to this region. The purpose of the present work is to study the IR scheme in the physical region at $`𝒪(Q^3)`$ and compare it with the HB approach in order to see whether an improved fit to the phase shift data can be obtained.

The organization of this paper is as follows. In Sec. II we review our notation for the effective Lagrangian and discuss the $`𝒪(Q^3)`$ calculation of the scattering amplitudes in the IR scheme and their reduction to HB form. Formulae for the $`\sigma `$ term and effective vertex couplings are also given. In Sec. III our fit to the phase shift data is described, first for the HB scheme and then for the IR approach. Our conclusions are presented in Sec. IV. Expressions for the IR integrals mentioned in the text are given in the Appendix, together with their HB reductions.

## II Formalism

### A Effective Lagrangian

In the IR or HB schemes the Feynman diagrams for $`\pi N`$ scattering follow Weinberg’s power counting rule . Define

$$\nu =1+2L+\underset{i}{}V_i\left(d_i+\frac{1}{2}n_i2\right),$$
(1)
where $`L`$ is the number of loops, $`V_i`$ is the number of vertices of type $`i`$ characterized by $`n_i`$ baryon fields and $`d_i`$ pion derivatives or pion mass factors. Then a given diagram is of leading order $`Q^\nu `$, where $`Q`$ is a small or “soft” momentum scale, for example, the pion mass $`m`$, the pion momentum or the mass splitting between the $`\mathrm{\Delta }`$-isobar and the nucleon, $`\delta =M_\mathrm{\Delta }M`$.

We briefly recapitulate the effective chiral Lagrangian given in I. Characterising it by $`d_i+\frac{1}{2}n_i`$, the order $`Q^2`$ part is

$`_2`$ $`=`$ $`\overline{N}(i\text{/}𝒟+g_\mathrm{A}\gamma ^\mu \gamma _5a_\mu M)N+\frac{1}{4}f_\pi ^2\mathrm{tr}(_\mu U^{}^\mu U)+\frac{1}{4}m^2f_\pi ^2\mathrm{tr}(U+U^{}2)`$ (3)
$`+\overline{\mathrm{\Delta }}_\mu ^a\mathrm{\Lambda }_{ab}^{\mu \nu }\mathrm{\Delta }_\nu ^b+h_\mathrm{A}\left(\overline{𝚫}_\mu \mathbf{}𝒂^\mu N+\overline{N}𝒂^\mu \mathbf{}𝚫_\mu \right)+\stackrel{~}{h}_\mathrm{A}\overline{\mathrm{\Delta }}_\mu ^a\gamma ^\nu \gamma _5a_\nu \mathrm{\Delta }_a^\mu ,`$
where the trace is taken over the isospin matrices. The isotriplet pion fields enter in the $`SU(2)`$ matrix

$$U(x)\xi ^2=\mathrm{exp}(2i\pi (x)/f_\pi ),$$
(4)
where $`f_\pi `$ is the pion decay constant and $`\pi (x)𝝅(x)\mathbf{}\frac{1}{2}𝝉`$, with $`𝝉`$ denoting the Pauli matrices. The axial vector field $`a_\mu (x)`$ and vector field $`v_\mu (x)`$ are defined by

$`a_\mu `$ $``$ $`{\displaystyle \frac{i}{2}}(\xi ^{}_\mu \xi \xi _\mu \xi ^{})=a_\mu ^{}=\frac{1}{2}𝒂_\mu \mathbf{}𝝉={\displaystyle \frac{1}{f_\pi }}_\mu \pi {\displaystyle \frac{1}{3f_\pi ^3}}\pi [\pi ,_\mu \pi ]+\mathrm{},`$ (5)
$`v_\mu `$ $``$ $`{\displaystyle \frac{i}{2}}(\xi ^{}_\mu \xi +\xi _\mu \xi ^{})=v_\mu ^{}=\frac{1}{2}𝒗_\mu \mathbf{}𝝉={\displaystyle \frac{i}{2f_\pi ^2}}\left(1{\displaystyle \frac{\pi ^2}{3f_\pi ^2}}\right)[\pi ,_\mu \pi ]+\mathrm{},`$ (6)
both of which contain one derivative. The covariant derivative on the nucleon field $`𝒟_\mu N=_\mu N+iv_\mu N`$. For the $`\mathrm{\Delta }`$ resonance we have introduced an isovector field $`𝚫_\mu =𝑻\mathrm{\Delta }_\mu `$ in terms of the standard $`2\times 4`$ isospin $`\frac{3}{2}`$ to $`\frac{1}{2}`$ transition matrix and the labels $`a`$ and $`b`$ in Eq. (3) are isospin indices. The kernel tensor in the $`\mathrm{\Delta }`$ kinetic energy term is

$$\mathrm{\Lambda }^{\mu \nu }=(i\text{/}𝒟M_\mathrm{\Delta })g^{\mu \nu }+i(\gamma ^\mu 𝒟^\nu +\gamma ^\nu 𝒟^\mu )\gamma ^\mu (i\text{/}𝒟+M_\mathrm{\Delta })\gamma ^\nu ,$$
(7)
suppressing isospin indices, and the covariant derivative is defined by

$$𝒟_\mu 𝚫_\nu =_\mu 𝚫_\nu +iv_\mu 𝚫_\nu 𝒗_\mu \times 𝚫_\nu .$$
(8)
We have chosen the simplest form for the $`\pi N\mathrm{\Delta }`$ interaction in Eq. (3) since modifications of the type suggested by Pascalutsa  or modifications of the standard off-shell $`Z`$ parameter can be absorbed in the other parameters of the Lagrangian .

The order $`Q^3`$ and $`Q^4`$ parts of the Lagrangian are:

$`_3`$ $`=`$ $`{\displaystyle \frac{\beta _\pi }{M}}\overline{N}N\mathrm{tr}(_\mu U^{}^\mu U){\displaystyle \frac{\kappa _\pi }{M}}\overline{N}v_{\mu \nu }\sigma ^{\mu \nu }N`$ (10)
$`+{\displaystyle \frac{\kappa _1}{2M^2}}i\overline{N}\gamma _\mu \underset{\nu }{\overset{}{𝒟}}N\mathrm{tr}(a^\mu a^\nu )+{\displaystyle \frac{\kappa _2}{M}}m^2\overline{N}N\mathrm{tr}(U+U^{}2)+\mathrm{},`$
$`_4`$ $`=`$ $`{\displaystyle \frac{\lambda _1}{M}}m^2\overline{N}\gamma _5(UU^{})N+{\displaystyle \frac{\lambda _2}{M^2}}\overline{N}\gamma ^\mu D^\nu v_{\mu \nu }N`$ (14)
$`+{\displaystyle \frac{\lambda _3}{M^2}}m^2\overline{N}\gamma _\mu [a^\mu ,UU^{}]N+{\displaystyle \frac{\lambda _4}{2M^3}}i\overline{N}\sigma _{\rho \mu }\underset{\nu }{\overset{}{𝒟}}N\mathrm{tr}(a^\rho D^\mu a^\nu )`$
$`+{\displaystyle \frac{\lambda _5}{16M^4}}i\overline{N}\gamma _\rho \{\underset{\mu }{\overset{}{𝒟}},\underset{\nu }{\overset{}{𝒟}}\}\tau ^aN\mathrm{tr}(\tau ^a[D^\rho a^\mu ,a^\nu ])`$
$`+{\displaystyle \frac{\lambda _6}{2M^2}}m^2\left\{\overline{𝚫}_\mu \mathbf{}\mathrm{tr}[i^\mu (UU^{})𝝉]N+\overline{N}\mathrm{tr}[i^\mu (UU^{})𝝉]\mathbf{}𝚫_\mu \right\}+\mathrm{}.`$
Here we have used the definitions

$`\underset{\mu }{\overset{}{𝒟}}`$ $`=`$ $`𝒟_\mu (\underset{\mu }{\overset{}{}}iv_\mu );v_{\mu \nu }=_\mu v_\nu _\nu v_\mu +i[v_\mu ,v_\nu ]=i[a_\mu ,a_\nu ],`$ (15)
$`D_\mu a_\nu `$ $`=`$ $`_\mu a_\nu +i[v_\mu ,a_\nu ];D_\sigma v_{\mu \nu }=_\sigma v_{\mu \nu }+i[v_\sigma ,v_{\mu \nu }],`$ (16)
In Eqs. (10) and (14) the ellipsis represents terms that do not contribute to the $`\pi N`$ scattering amplitude.

### B Scattering Amplitudes

Following the standard notation of Höhler  for $`\pi N`$ scattering we write the $`T`$ matrix as

$$T_{ba}\pi _b|T|\pi _a=T^+\delta _{ab}+\frac{1}{2}[\tau _b,\tau _a]T^{},$$
(17)
where the isospin symmetric and antisymmetric amplitudes are

$$T^\pm =\overline{u}(p^{})\left[A^\pm +\frac{1}{2}(\text{/}q+\text{/}q^{})B^\pm \right]u(p),$$
(18)
and $`u`$ denotes a nucleon spinor. Here, as shown in Fig. 1, $`q`$ and $`q^{}`$ are the c.m. momenta of the incoming and outgoing pions with isospin labels $`a`$ and $`b`$ respectively. The c.m. momenta of the incoming and outgoing nucleons are labelled $`p`$ and $`p^{}`$ respectively. The amplitudes $`A^\pm `$ and $`B^\pm `$ are functions of the Mandelstam invariants $`s=(p+q)^2`$, $`t=(qq^{})^2`$, and $`u=(pq^{})^2`$.

In Fig. 1 we show the tree level Feynman diagrams arising from the contact terms and from one nucleon exchange. The vertex in Fig. 1(a) arises from any of the interactions in $`_3`$ and $`_4`$ (except for the $`\lambda _1`$ and $`\lambda _6`$ terms), as well as the Weinberg term $`\overline{N}\gamma ^\mu v_\mu N`$ of $`_2`$. The amplitudes for these tree diagrams were given in I.

#### 1 $`\mathrm{\Delta }`$ Exchange

When the $`\mathrm{\Delta }`$ appears as an intermediate state for $`\pi N`$ scattering, the tree-level $`T`$-matrix diverges at $`s=M_\mathrm{\Delta }^2`$. In I we argued that the power counting should work only for irreducible diagrams. Therefore we summed diagrams containing the one-particle irreducible self-energy insertions shown in Fig. 2 to all orders so as to replace the free propagator by the dressed propagator which is finite. To $`𝒪(Q^3)`$ this gave for the real part of the $`\mathrm{\Delta }`$ propagator in $`d`$ dimensions

$`\mathrm{}G_{\mu \nu }(k)`$ $`=`$ $`\mathrm{\Delta }_R(k^2)(\text{/}k+M_\mathrm{\Delta })(P^{3/2})_{\mu \nu }{\displaystyle \frac{1}{\sqrt{d1}M_\mathrm{\Delta }}}(P_{12}^{1/2}+P_{21}^{1/2})_{\mu \nu }`$ (20)
$`+{\displaystyle \frac{(d2)}{(d1)M_\mathrm{\Delta }^2}}(\text{/}k+M_\mathrm{\Delta })(P_{22}^{1/2})_{\mu \nu },`$
where the spin projection operators , denoted by $`\left(P_{ij}^I\right)_{\mu \nu }`$, are generalized by replacing factors of 2 and 3 by $`(d2)`$ and $`(d1)`$, respectively. In Eq. (20) we have defined

$$\mathrm{\Delta }_R(k^2)=\frac{k^2M_\mathrm{\Delta }^2\mathrm{\Pi }_\mathrm{\Delta }(k^2)}{[k^2M_\mathrm{\Delta }^2\mathrm{\Pi }_\mathrm{\Delta }(k^2)]^2+M_\mathrm{\Delta }^2\mathrm{\Gamma }_\mathrm{\Delta }^2(k^2)}.$$
(21)
(Obviously the bare propagator corresponds to setting $`\mathrm{\Pi }_\mathrm{\Delta }`$ and $`\mathrm{\Gamma }_\mathrm{\Delta }`$ to zero.)

The diagrams which contribute to the self energy are shown in Fig. 2, where the open box denotes a free $`\mathrm{\Delta }`$ propagator. Only the $`g_{\mu \nu }`$ part is needed in leading order and the renormalized self energy is

$$\mathrm{\Sigma }_\mathrm{\Delta }^{\mathrm{ren}}(k)=\mathrm{\Sigma }_\mathrm{\Delta }(k)\mathrm{}\mathrm{\Sigma }_\mathrm{\Delta }(k)|_{/k=M_\mathrm{\Delta }}\frac{}{\text{/}k}\mathrm{}\mathrm{\Sigma }_\mathrm{\Delta }(k)|_{/k=M_\mathrm{\Delta }}(\text{/}kM_\mathrm{\Delta }).$$
(22)
The integrals which arise are evaluated in the IR scheme, which is briefly discussed in the Appendix. For the $`\gamma `$-matrix algebra we evaluate terms whose leading contribution is at $`𝒪(Q^3)`$ and we discard terms for which the leading contribution is of higher order. We include polynomial terms obtained from the product of a $`1/(d4)`$ singularity with $`(d4)`$ factors which were were dropped in I. Following this procedure the real part of the $`\mathrm{\Delta }`$ polarization in the $`s`$\- or $`u`$-channel is given by

$`\mathrm{\Pi }_\mathrm{\Delta }(s)`$ $`=`$ $`{\displaystyle \frac{h_A^2}{(4\pi f_\pi )^2}}\{(s+2MM_\mathrm{\Delta }+M_\mathrm{\Delta }^2)[\overline{}_{11}^{(2)}(s,M)\overline{}_{11}^{(2)}(M_\mathrm{\Delta }^2,M)]`$ (27)
$`2(sM_\mathrm{\Delta }^2)M_\mathrm{\Delta }(M_\mathrm{\Delta }+M){\displaystyle \frac{}{M_\mathrm{\Delta }^2}}\overline{}_{11}^{(2)}(M_\mathrm{\Delta }^2,M)\}`$
$`{\displaystyle \frac{5\stackrel{~}{h}_A^2}{36(4\pi f_\pi )^2}}\{2(s+9M_\mathrm{\Delta }^2)[\overline{}_{11}^{(2)}(s,M_\mathrm{\Delta })\overline{}_{11}^{(2)}(M_\mathrm{\Delta }^2,M_\mathrm{\Delta })]`$
$`20(sM_\mathrm{\Delta }^2)M_\mathrm{\Delta }^2{\displaystyle \frac{}{k^2}}\overline{}_{11}^{(2)}(k^2,M_\mathrm{\Delta })|_{k^2=M_\mathrm{\Delta }^2}`$
$`+{\displaystyle \frac{13(sM_\mathrm{\Delta }^2)}{36s^2}}((s+3M_\mathrm{\Delta }^2)[6sm^2(sM_\mathrm{\Delta }^2)^2]24s^2m^2)\},`$
where the integral $`\overline{}_{11}^{(2)}`$ is specified in the Appendix. By suitably approximating the above expression the HB form given in I is obtained (apart from polynomial terms). The width is zero in the $`u`$-channel and in the $`s`$-channel it is

$$\mathrm{\Gamma }_\mathrm{\Delta }(s)=\frac{\pi h_A^2}{12M_\mathrm{\Delta }s^2(4\pi f_\pi )^2}\left[(sM^2+m^2)^24sm^2\right]^{\frac{3}{2}}(s+M^2+2MM_\mathrm{\Delta }m^2),$$
(28)
which gives the exact result  on shell at $`s=M_\mathrm{\Delta }^2`$.

The $`\mathrm{\Delta }`$ exchange tree-level diagram for $`\pi N`$ scattering is pictured in Fig. 3, where the solid box denotes the dressed propagator discussed above (a similar notation is used in Fig. 5 below). The corresponding contribution to the $`T`$-matrix can be obtained from the expression given in I.

Including also the loop diagrams, discussed below, the real part of the total $`T`$-matrix yields the real part of the elastic scattering amplitude, $`f_\alpha `$, by means of the standard partial wave expansion. Here the isospin-spin partial wave channels are labelled by $`\alpha (l,2I,2J)`$ with $`l`$ the orbital angular momentum, $`I`$ the total isospin, and $`J=l\pm \frac{1}{2}`$ the total angular momentum. The phase shifts $`\delta _\alpha `$ are then given by

$$\mathrm{}f_\alpha =\frac{1}{|𝒒|}\mathrm{}e^{i\delta _\alpha }\mathrm{sin}\delta _\alpha ,$$
(29)
where $`|𝒒|`$ is the magnitude of the c.m. three-momentum. We shall refer to this as the “$`S`$-matrix” approach.

An alternative, which has been espoused by Fettes and Meißner , is the $`K`$-matrix approach. This is introduced by setting

$$f_\alpha =\frac{K_\alpha }{1i|𝒒|K_\alpha }\mathrm{where}K_\alpha =\frac{1}{|𝒒|}\mathrm{tan}\delta _\alpha .$$
(30)
The calculated amplitude $`\mathrm{}f_\alpha `$ is then assumed to actually be $`K_\alpha `$ so that an infinite $`f_\alpha `$ on resonance corresponds to $`\delta =\frac{1}{2}\pi `$ (for further discussion see ). This allows the free $`\mathrm{\Delta }`$ propagator to be used everywhere. Thus it is not necessary to sum the self energy insertions of Fig. 2 and we can expand $`\mathrm{\Delta }_R`$ to first order in $`\mathrm{\Pi }_\mathrm{\Delta }`$. The zeroth order term gives the tree diagram of Fig. 3, but with the free propagator, and the first order term is treated in the same way as the other third order loop diagrams. We will compare this method with the $`S`$-matrix approach.

#### 2 One-Loop Diagrams

A set of one-loop diagrams that contribute to the $`\pi N`$ $`T`$-matrix at $`𝒪(Q^3)`$ is shown in Fig. 4. Here we can use the free $`\mathrm{\Delta }`$ propagator, denoted by an open box, since no singularities are generated in the $`T`$-matrix. We illustrate our procedure by discussing the evaluation of Fig. 4(g), which has the same value as Fig. 4(h). Using the standard Feynman rules, we find

$`T_{4g}`$ $`=`$ $`\overline{u}(p^{}){\displaystyle \frac{g_A^2h_A^2}{4f_\pi ^4}}i\mu ^{4d}{\displaystyle _I}{\displaystyle \frac{d^d\mathrm{}}{(2\pi )^d}}\mathrm{}^\nu T^cG_{\nu \mu }(p^{}+\mathrm{})T^bq^\mu `$ (32)
$`\times {\displaystyle \frac{1}{\text{/}p+\text{/}q+\text{/}\mathrm{}M+iϵ}}\text{/}q\gamma _5\tau ^a{\displaystyle \frac{1}{\text{/}p+\text{/}\mathrm{}M+iϵ}}\text{/}\mathrm{}\gamma _5\tau ^c{\displaystyle \frac{1}{\mathrm{}^2m^2+iϵ}}u(p),`$
where $`\mu `$ is the renormalization scale, $`T^i`$ denotes a component of the isospin $`\frac{3}{2}`$ to $`\frac{1}{2}`$ transition matrix and the subscript $`I`$ denotes that the integral is to be evaluated using infrared regularization. Rather than using the projection operators in (20) it is more convenient to write the free $`\mathrm{\Delta }`$ propagator in the form

$$G_{\nu \mu }(k)=\frac{1}{\text{/}kM_\mathrm{\Delta }+iϵ}\left[g_{\nu \mu }+\frac{1}{d1}\gamma _\nu \gamma _\mu +\frac{1}{(d1)M_\mathrm{\Delta }}(\gamma _\nu k_\mu k_\nu \gamma _\mu )+\frac{(d2)}{(d1)M_\mathrm{\Delta }^2}k_\nu k_\mu \right].$$
(33)
To $`𝒪(Q^3)`$ the $`\text{/}\mathrm{}`$ terms in the baryon denominators can be dropped. Then the integral involves $`\mathrm{}^\nu \mathrm{}^\rho `$ in the numerator and to this order only the $`g_{\nu \rho }`$ contribution is needed. Carrying out the $`\gamma `$-matrix and isospin algebra, we obtain

$`T_{4g}`$ $`=`$ $`\overline{u}(p^{})\frac{1}{2}[\tau ^b,\tau ^a]{\displaystyle \frac{4g_A^2h_A^2}{9f_\pi ^2(4\pi f_\pi )^2}}M^2(M+M_\mathrm{\Delta })(sM^2m^2+\frac{1}{2}t2M\text{/}q)`$ (35)
$`\times \left[\overline{}_{13}^{(3)}(s,t,M_\mathrm{\Delta },M,M){\displaystyle \frac{1}{24M^4}}(sM_\mathrm{\Delta }^2m^2)\left(\frac{4}{3}+32\pi ^2L+\mathrm{ln}{\displaystyle \frac{m^2}{\mu ^2}}\right)\right]u(p).`$
The integral $`\overline{}_{13}^{(3)}`$ is defined in the Appendix, as is the divergent quantity $`L`$. In the usual way the terms involving $`\left(2L+\frac{1}{(4\pi )^2}\mathrm{ln}\frac{m^2}{\mu ^2}\right)`$ can be absorbed in the low energy constants so that the result is independant of the renormalization scale. This expression can be reduced to HB form. Here, as in I, we define the average baryon mass, $`\overline{M}=\frac{1}{2}(M+M_\mathrm{\Delta })`$ and $`\eta =(sM^2)/(2\overline{M})`$. Then evaluating the expressions to leading order and using the reduction of $`\overline{}_{13}^{(3)}`$ given in the Appendix, we obtain

$$T_{4g}\overline{u}(p^{})\frac{1}{2}[\tau ^b,\tau ^a]\frac{4g_A^2h_A^2}{27f_\pi ^2(4\pi f_\pi )^2}\left[M\text{/}q\overline{M}\eta +\frac{1}{4}(2m^2t)\right]\left[F_2(\eta ,\delta )\frac{2}{3}(\eta \delta )\right]u(p),$$
(36)
recalling that $`\delta `$ denotes the delta-nucleon mass difference. This agrees with the result given in I, apart from the polynomial term which we include here, but which was subsumed in the low energy constants in I. To these results in the $`s`$-channel should be added the contribution in the $`u`$-channel with the replacement $`su`$ and the interchanges $`ab`$ and $`qq^{}`$.

Proceeding in similar fashion the real parts of the $`T`$-matrix for the diagrams of Fig. 4 may be calculated. Note that the diagrams in Fig. 4(m) and (n) involve the nucleon one-loop self-energy which is real for the energies of interest here. As with the $`\mathrm{\Delta }`$, we make on-shell mass and wavefunction counterterm subtractions to obtain the renormalized self-energy, namely

$$\mathrm{\Sigma }_N^{\mathrm{ren}}(k)=\mathrm{\Sigma }_N(k)\mathrm{\Sigma }_N(k)|_{/k=M}\frac{}{\text{/}k}\mathrm{\Sigma }_N(k)|_{/k=M}(\text{/}kM).$$
(37)
Since these diagrams do not give singular contributions to the $`T`$-matrix we do not sum these self-energy insertions.

We also need to evaluate in similar fashion the diagrams of Fig. 5, where the solid boxes denote the dressed propagator discussed in Subsec. II.B.1. The diagrams of Fig. 5(a)–(f) modify the $`\pi NN`$ tree vertex in $`_2`$, while diagrams (g)–(l) similarly modify the $`\pi N\mathrm{\Delta }`$ tree vertex.

In a few cases the IR results can be checked against Ref. . The reduction of the IR expressions to HB form also provides a check since the results of I should be reproduced (these, in turn, were checked against Mojžiš  for the cases without a $`\mathrm{\Delta }`$). In the course of checking we found a phase error in one of the HB results in I. This means that the expression given in I for diagram 4(i), which is equal to that for diagram 4(j), should be multiplied by a factor of 5. This has little impact on the fits presented in I, but the parameters are modified somewhat. Since the IR results will not prove to be satisfactory, we shall not list the rather lengthy expressions for all the diagrams; they are available on request.

### C $`\sigma `$ Term and Effective Couplings

We may obtain the nucleon $`\sigma `$ term from the Feynman-Hellman theorem,

$$\sigma (0)=m^2\frac{M}{m^2}.$$
(38)
The nucleon mass receives contributions from the $`\kappa _2`$ term in $`_3`$ and diagrams (m) and (n) of Fig. 4. In the IR scheme these yield

$`\sigma (0)`$ $`=`$ $`{\displaystyle \frac{4\kappa _2m^2}{M}}{\displaystyle \frac{9g_A^2Mm^2}{4(4\pi f_\pi )^2}}\overline{}_{11}(M^2,M)`$ (40)
$`{\displaystyle \frac{2h_A^2(M+M_\mathrm{\Delta })m^2}{(4\pi f_\pi )^2}}\left[\overline{}_{11}(M^2,M_\mathrm{\Delta })+{\displaystyle \frac{(M^2M_\mathrm{\Delta }^2)(M^2+2M_\mathrm{\Delta }^2)}{6M^2M_\mathrm{\Delta }^2}}\right].`$
In the HB approximation this reduces to

$$\sigma (0)\frac{4\kappa _2m^2}{M}\frac{9\pi g_A^2m^3}{4(4\pi f_\pi )^2}\frac{4h_A^2m^2}{(4\pi f_\pi )^2}[J(\delta )\delta ],$$
(41)
where the integral $`J`$ is defined in I, see also the Appendix. This agrees with the result of Fettes and Meißner  and, modulo a polynomial contribution in $`m^2\delta `$, with I.

The $`\pi NN`$ vertex up to one-loop order consists of the tree vertex generated from the axial $`a_\mu `$ term in $`_2`$ and the one-loop diagrams shown in the upper part of Fig. 5 – diagrams (a) to (f). As in I, we can calculate the one-loop vertex function $`\mathrm{\Gamma }^a(k,k^{},q)`$, where $`k`$($`k^{}`$) is the incoming (outgoing) momentum of the nucleon and $`q=k^{}k`$ is the momentum transfer. The $`\pi NN`$ coupling for on-shell nucleons is then obtained from

$$\overline{u}(k^{})\mathrm{\Gamma }^a(k,k^{},q)u(k)=g_{\pi NN}(q^2)\overline{u}(k^{})\gamma _5\tau ^au(k).$$
(42)
At zero momentum transfer we obtain, to $`𝒪(Q^2)`$,

$`\zeta _{\pi NN}`$ $``$ $`{\displaystyle \frac{g_{\pi NN}(0)f_\pi }{Mg_A}}=1{\displaystyle \frac{2m^2\lambda _1}{g_AM^2}}+{\displaystyle \frac{g_A^2m^2}{2(4\pi f_\pi )^2}}+{\displaystyle \frac{64h_A^2m^2M}{27(M_\mathrm{\Delta }M)(4\pi f_\pi )^2}}`$ (47)
$`\times \{\overline{}_{11}(M^2,M)+{\displaystyle \frac{[(M^2M_\mathrm{\Delta }^2)^22m^2(M^2+M_\mathrm{\Delta }^2)]}{4m^2M^2}}\overline{}_{11}(M^2,M_\mathrm{\Delta })`$
$`+{\displaystyle \frac{(M_\mathrm{\Delta }^2M^2)}{24m^2M^4}}[(M^2M_\mathrm{\Delta }^2)^26m^2M^2]\}`$
$`{\displaystyle \frac{10h_A^2\stackrel{~}{h}_A(M+M_\mathrm{\Delta })}{81g_AM_\mathrm{\Delta }^2(4\pi f_\pi )^2}}\{[7MM_\mathrm{\Delta }^2+8M_\mathrm{\Delta }^3M^2(2M+3M_\mathrm{\Delta })]\overline{}_{11}^{(1)}(M^2,M_\mathrm{\Delta })`$
$`+{\displaystyle \frac{1}{12M^4}}[35MM_\mathrm{\Delta }^2+37M_\mathrm{\Delta }^32M^2(2M+3M_\mathrm{\Delta })][2m^2M^2(M^2M_\mathrm{\Delta }^2)^2]\}.`$
With the parameters obtained from fits to the $`\pi N`$ phase shifts this allows a test of the Goldberger-Treiman relation.

Similarly, we can calculate the $`\pi N\mathrm{\Delta }`$ vertex from the tree-level $`h_A`$ term and the one-loop diagrams shown in the lower part of Fig. 5 – diagrams (g) to (l). Here, in the vertex function $`\mathrm{\Gamma }^{\mu a}(k,k^{},q)`$ the label $`k`$ now refers to the incoming $`\mathrm{\Delta }`$ momentum. The $`\pi N\mathrm{\Delta }`$ coupling is obtained from

$$\overline{u}(k^{})\mathrm{\Gamma }^{\mu a}(k,k^{},q)u_\mu (k)=g_{\pi N\mathrm{\Delta }}(q^2)\overline{u}(k^{})q^\mu T^au_\mu (k),$$
(48)
where $`u_\mu (k)`$ is the $`\mathrm{\Delta }`$ spinor. At zero momentum transfer we obtain, to $`𝒪(Q^2)`$,

$`{\displaystyle \frac{g_{\pi N\mathrm{\Delta }}(0)f_\pi }{Mh_A}}`$ $`=`$ $`1{\displaystyle \frac{2m^2\lambda _6}{h_AM^2}}{\displaystyle \frac{4g_A^2M^2}{(4\pi f_\pi )^2}}\overline{}_{12}^{(3)}(M_\mathrm{\Delta }^2,M,M)|_{m=0}{\displaystyle \frac{65g_A\stackrel{~}{h}_A(M+M_\mathrm{\Delta })m^2}{108M(4\pi f_\pi )^2}}`$ (52)
$`{\displaystyle \frac{2h_A^2M(M+M_\mathrm{\Delta })}{9(4\pi f_\pi )^2}}\{\overline{}_{12}^{(3)}(M_\mathrm{\Delta }^2,M,M_\mathrm{\Delta })|_{m=0}+{\displaystyle \frac{\left[(M_\mathrm{\Delta }^2M^2)^26m^2M^2\right]}{18M^4}}\}`$
$`{\displaystyle \frac{50\stackrel{~}{h}_A^2M_\mathrm{\Delta }(M+M_\mathrm{\Delta })}{81(4\pi f_\pi )^2}}\{\overline{}_{12}^{(3)}(M_\mathrm{\Delta }^2,M_\mathrm{\Delta },M_\mathrm{\Delta })|_{m=0}`$
$`+{\displaystyle \frac{(15M+M_\mathrm{\Delta })[(M_\mathrm{\Delta }^2M^2)^26m^2M^2]}{360M^4M_\mathrm{\Delta }}}\}.`$
Note that for the first two integrals $`\overline{}_{12}^{(3)}`$ both real and imaginary parts need to be considered. The value of $`g_{\pi N\mathrm{\Delta }}`$ is complex because the intermediate pion and nucleon states for Fig. 5(i) and (j) can go on shell. We therefore define $`\zeta _{\pi N\mathrm{\Delta }}|g_{\pi N\mathrm{\Delta }}(0)f_\pi /(Mh_A)|`$. The HB reduction of Eqs. (47) and (52) agrees with the expressions in I up to polynomial terms.

## III Results

As fixed input parameters, we use the standard baryon and pion masses: $`M=939`$MeV, $`M_\mathrm{\Delta }=1232`$MeV, and $`m=139`$MeV. We also take $`f_\pi =92.4`$MeV from charged pion decay, $`g_\mathrm{A}=1.26`$ from neutron $`\beta `$ decay, and $`h_\mathrm{A}=1.46`$ from Eq. (28) for $`\mathrm{\Gamma }_\mathrm{\Delta }(M_\mathrm{\Delta }^2)`$ using the central value of the $`\mathrm{\Delta }`$ width, $`\mathrm{\Gamma }_\mathrm{\Delta }=120\pm 5`$MeV. We have ten low energy constants: $`\beta _\pi `$, $`\kappa _\pi `$, $`\kappa _1`$, $`\kappa _2`$, and $`\lambda _1`$ to $`\lambda _6`$, plus the $`\pi \mathrm{\Delta }\mathrm{\Delta }`$ coupling, $`\stackrel{~}{h}_A`$. These are obtained by optimizing the fit of our calculated $`\pi N`$ $`S`$\- and $`P`$-wave phase shifts to the 2001 data of the VPI/GW group . Since errors are not given we compared results obtained by assigning all the data points the same relative weight to those obtained using weightings suggested in the literature . The differences were small so we used the same relative weightings, minimizing

$$\chi ^2=\frac{1}{N}\underset{i=1}{\overset{N}{}}\left(\frac{\xi _i^{\mathrm{theory}}\xi _i^{\mathrm{experiment}}}{0.08\xi _i^{\mathrm{experiment}}}\right)^2,$$
(53)
where $`\xi _i=\mathrm{sin}2\delta _i`$ for the $`S`$-matrix calculation and $`\xi _i=2\mathrm{tan}\delta _i`$ for the $`K`$-matrix calculation.

### A HB Results

Our fits to the $`S`$\- and $`P`$-wave phase shifts in the HB calculation are shown in Figs. 6 and 7 using the $`S`$-matrix and $`K`$-matrix methods, respectively. Here the solid curves correspond to an unconstrained fit. As in I and Ref.  it was possible to fit the phase shift data out to a pion c.m. kinetic energy, $`ϵ`$, of 100 MeV, slightly below the $`\mathrm{\Delta }`$ resonance at $`ϵ=127`$ MeV. More precisely the data were fitted from $`\sqrt{s}=1090`$ to 1200 MeV in 3 MeV steps. The short horizontal section in Fig. 6 indicates a small region where unitarity, which is only enforced perturbatively, is slightly violated and a phase shift cannot be determined; this, of course, does not occur in the $`K`$-matrix approach, by definition. The corresponding parameters are collected in Table I where, in common with previous work, it is observed that the deduced values of the $`\sigma `$ term are much larger than the oft-quoted figure of $`45\pm 8`$ due to Gasser et al. . Therefore we performed two additional fits where $`\kappa _2`$ was constrained to produce values of the $`\sigma `$ term of 75 MeV and 45 MeV. These are indicated, respectively, by the dotted and dashed curves in Figs. 6 and 7. For the $`K`$-matrix case the dotted curve is barely distinuishable from the solid curve since the change in the $`\sigma `$ term is small here.

As expected the $`S`$-matrix fits in Fig. 6 are quite similar to those given in I. The fitted parameters, however, differ since here we include the polynomial terms and we subsume the renormalization scale dependance into the low energy constants. In many of the partial waves the $`K`$-matrix fits (Fig. 7) are a little better, particularly for $`\delta _{P33}`$, however $`\delta _{P31}`$ begins to deviate significantly from the data at $`ϵ80`$ MeV. Something of the same trend is visible in the work of Fettes and Meißner . The $`𝒪(Q^2)`$ low energy constants given there are related to our constants by $`c_1=\kappa _2/M`$, $`c_2=\kappa _1/(2M)`$, $`c_3=2\beta _\pi /M`$ and $`c_4=\kappa _\pi /M`$. Their fit $`1^{}`$ to Matsinos’ data  produces values which are fairly similar to ours when the different value of the $`\pi N\mathrm{\Delta }`$ coupling, $`h_A`$, and the constraint $`\sigma (0)=59`$ MeV are taken into account. (We should note, however, that substantial differences are seen for their fit to the Karlsruhe data .) The $`𝒪(Q^2)`$ constants in Table I are reasonably similar in the $`S`$\- and $`K`$-matrix calculations, however sizeable differences appear at $`𝒪(Q^3)`$. Further they are not of natural size (unity). These features may be due to the strong cancellations which occur among the various terms and the fact that the parameters are not uniquely determined since it is possible to obtain similar values of $`\chi ^2`$ when $`\sigma (0)`$ is constrained. A significantly larger $`\chi ^2`$ is only obtained for the $`S`$-matrix case with $`\sigma (0)=45`$ MeV.

The ratios of the effective $`\pi NN`$ coupling constant to the bare value, $`\zeta _{\pi NN}`$, in Table I show a Goldberger-Treiman discrepancy of a few percent. Schröder et al.  have recently determined a precise result, namely $`\zeta _{\pi NN}=1.027_{0.008}^{+0.012}`$. A similar value is obtained from the latest $`\pi NN`$ coupling constant obtained by the George Washington University/TRIUMF group . These values favor a $`\sigma `$ term of between 75 and 45 MeV which is in accord with several recent analyses. In Ref.  64 MeV is quoted with an error of about 10%. A sum rule determination by Olsson  gave 55 MeV with a 16% error. Finally Schröder et al.  have indicated that the value of 45 MeV extracted  from the Karlsruhe data should be increased by 13 MeV giving 58 MeV. Lastly we note that since $`\zeta _{\pi N\mathrm{\Delta }}`$ remains fairly close to unity the coupling is not changed too much from the bare value, which seems intuitively reasonable.

It is also interesting to examine the threshold results. We give in Table II the $`S`$-wave isoscalar and isovector scattering lengths, $`b_0=(a_1+2a_3)/3`$ and $`b_1=(a_3a_1)/3`$ (with the notation $`a_{2I}`$) and the $`P`$-wave scattering volumes ($`a_{2I\mathrm{\hspace{0.17em}2}J}`$). The experimental values come from various sources. The $`S`$-wave isoscalar scattering lengths were obtained  from the pionic atom data  with an improved treatment of the $`\pi `$-$`d`$ scattering length. The $`I=\frac{1}{2}`$ $`P`$-wave scattering volumes are the old Karlsruhe values , while the $`I=\frac{3}{2}`$ volumes are from a recent analysis by Fettes and Matsinos . Our calculated results for $`b_0`$ at the extreme values of $`\sigma (0)`$ are clearly not favored by the data which prefer a value in the middle. Notice that, at this level of accuracy, there is a noticeable difference between the $`S`$\- and $`K`$-matrix results. The value of $`a_{33}`$ appears to be somewhat low for the cases where $`\sigma (0)=45`$ MeV. Apart from this, the remaining $`P`$-wave results and the values of $`b_1`$ show little sensitivity to the calculation employed and are in reasonable agreement with the data. Of course this is not surprising since the low energy phase shifts are included in the fit.

### B IR Results

Turning to the IR scheme the results are disappointing. It was only found possible to fit the phase shift data in the low energy regime out to $`ϵ=40`$ MeV. Specifically we fitted for $`\sqrt{s}`$ from 1090 to 1123 MeV in steps of 3 MeV. The IR $`S`$-matrix results are shown by the dashed curves in Fig. 8; here the restriction $`\stackrel{~}{h}_A0.5`$ was imposed, see the discussion below. The $`K`$-matrix results are denoted by dotted curves in the figure. For contrast the solid line shows the HB results for a fit over the same energy range; the HB results from the $`S`$\- and $`K`$-matrix methods are very similar and we choose to display the $`K`$-matrix results. While the HB phase shifts follow the data approximately for $`ϵ>40`$ MeV, in a number of cases the IR results rapidly diverge from the experimental values. The parameters corresponding to these results are collected in Table III. While the HB parameters differ from those given in Table I, particularly those of $`𝒪(Q^3)`$, we caution that they are not well determined by a fit over this limited energy range. Notice that, for the $`K`$-matrix cases, $`\chi ^2`$ is substantially larger in the IR scheme than the HB one. In the IR case $`\zeta _{\pi NN}`$ differs substantially from unity, implying a huge Goldberger-Treiman discrepancy. This is unphysical and is a signal of the contortions that the parameters are going through in order to fit over even this limited range of energies. The quantity $`\zeta _{\pi N\mathrm{\Delta }}`$ is also much less than unity implying an effective $`\pi N\mathrm{\Delta }`$ coupling very different from the bare value.

We have, of course, carried out a number of additional IR calculations with a view to improving the results. For example constraining the value of the $`\sigma `$ term to be 75 or 45 MeV increases the value of $`\chi ^2`$, but the results are qualitatively very similar to those discussed above. In order to see whether the problem was due to the $`\mathrm{\Delta }`$ resonance contributions we removed them entirely by setting the coupling $`h_A`$ to zero. As expected the parameters changed substantially, but the fit to the phase shifts remained qualitatively the same and in fact $`\chi ^2`$ increased somewhat. Thus we have not been able to qualitatively improve the IR results discussed in detail above.

In order to further investigate the reason for the problem with the IR approach we compare the IR polarization $`\mathrm{\Pi }_\mathrm{\Delta }`$ of Eq. (27) and a few IR integrals (defined in the Appendix) with their HB counterparts. In Table IV we make the comparison in the $`s`$-channel. For the most part the IR and HB results are quite comparable. The picture changes for the $`u`$-channel in Table V, where for a given $`s`$ we have chosen the minimum value, $`u_{\mathrm{min}}=2(M^2+m^2)s`$, so as to enhance the contrast. Here the IR integrals are much larger than the corresponding HB results, by a factor of 6 on average. Further $`\mathrm{\Pi }_\mathrm{\Delta }(u)`$ changes sign. Since $`\mathrm{\Delta }_R(u)=1/[uM_\mathrm{\Delta }^2\mathrm{\Pi }_\mathrm{\Delta }(u)]`$ and $`uM_\mathrm{\Delta }^2`$ is negative this change of sign can produce a pole; this was the reason for our restriction of $`\stackrel{~}{h}_A0.5`$ in the $`S`$-matrix IR fit. Of course the problem does not arise if $`\mathrm{\Pi }_\mathrm{\Delta }(u)`$ is expanded out of the denominator as in the $`K`$-matrix approach, although the mathematical justification for doing so is slight. The problem does not arise in the HB case where $`\mathrm{\Pi }_\mathrm{\Delta }(u)`$ is always positive.

It is also useful to make the point graphically by studying the basic integral involving one baryon and one meson propagator. Figure 9 shows $`\overline{}_{11}(s,M_\mathrm{\Delta })`$ as a function of $`\sqrt{s}`$. It is observed that the IR integral (solid line) and the HB approximation thereto (dashed curve) agree quite closely. Figure 10 shows the corresponding result in the $`u`$-channel; $`\overline{}_{11}(u_{\mathrm{min}},M_\mathrm{\Delta })`$ as a function of $`\sqrt{s}`$. It is seen that the magnitude of the IR integral (solid line) is much larger than that of its HB counterpart (dashed curve); note the ordinate scale here in comparison to Fig. 9. The IR integral is governed by the value of $`x`$ (see Appendix) which in the $`s`$-channel is $`(sM_\mathrm{\Delta }^2+m^2)/(2m\sqrt{s})`$. In the HB approximation this is replaced by $`(sM_\mathrm{\Delta }^2)/(2m\overline{M})=(\eta \delta )/m`$, in the notation of I. At $`\sqrt{s}=1200`$ MeV, $`x`$ is $`0.18`$ ($`0.26`$) in the IR (HB) case. Correspondingly in the $`u`$-channel, at $`\sqrt{u_{\mathrm{min}}}=601.7`$ MeV, $`x`$ is $`6.79`$ ($`3.83`$) in the IR (HB) case and this large difference is the heart of the problem with the IR approach. However, it is disquieting that the rather reasonable HB results are achieved by the replacement of $`\sqrt{u}`$ in the denominator by the square root of average baryon mass, $`\sqrt{\overline{M}}`$, since this can be quite a rough approximation.

## IV Conclusions

We have performed a least-squares fit to the $`S`$\- and $`P`$-wave phase shift data using chiral perturbation theory to $`𝒪(Q^3)`$ with explicit $`\pi `$, $`N`$ and $`\mathrm{\Delta }`$ fields. The heavy baryon formulation used here represents an improvement compared to I in that the polynomial terms were included explicitly and the renormalization scale dependance was subsumed into the low energy constants. We contrasted an “$`S`$-matrix” method (as in I) with a $`K`$-matrix analysis (as in Ref. ). The former seems to be more physical since it employs dressed $`\mathrm{\Delta }`$ propagators which are finite, while the latter uses bare propagators which give a divergent scattering amplitude at the $`\mathrm{\Delta }`$ resonance energy. Broadly speaking the two methods give similar results and allow a fit to the phase shift up to a c.m. pion kinetic energy, $`ϵ`$, of 100 MeV. In detail, however, there are some differences. For instance, the value of the $`\sigma `$ term (if unconstrained) and the Goldberger-Treiman discrepancy are smaller with the $`K`$-matrix.

Our main interest was to see whether the infrared regularization scheme, which is manifestly Lorentz invariant, could provide an equally satisfactory account of the phase shift data. It could not. We were only able to fit the data up to $`ϵ=40`$ MeV and the resulting Goldberger-Treiman discrepancy was ridiculously large. The difficulty was not due to the explicit inclusion of the $`\mathrm{\Delta }`$ resonance, but rather stemmed from the $`u`$ channel where the magnitudes of the integrals were much larger than the corresponding heavy baryon ones. This leads to uneasiness regarding the latter since they are obtained by approximating the former. The heavy baryon scheme thus pushes significant contributions into higher orders where they are combined with the other terms of that order. The only information we have as to whether the net effect is small comes from the comparison of third and fourth orders carried out by Fettes and Meißner . In some cases there are sizeable differences, although these authors suggest that in most cases the chiral series appears convergent. Overall, in our opinion, the description of pion-nucleon scattering via chiral perturbation theory is not yet in a satisfactory state.

## Acknowledgements

We acknowledge support from the Department of Energy under grant No. DE-FG02-87ER40328.

##

Integrals involving only meson propagators need no discussion since the standard results are used and integrals involving only baryon propagators give no contribution in the IR or HB schemes. Infrared regularization is needed for integrals involving both baryon and meson propagators. Here we give only those integrals that are referred to in the main text, together with their HB reduction. The integrals needed can be derived from the basic infrared integral which has been discussed by Becher and Leutwyler  (note that our conventions differ from theirs), see also I. Explicitly this is

$`i\mu ^{4d}{\displaystyle _I}{\displaystyle \frac{d^d\mathrm{}}{(2\pi )^d}}{\displaystyle \frac{1}{[(\mathrm{}+k)^2^2+iϵ][\mathrm{}^2m^2+iϵ]}}`$ (54)
$`=i\mu ^{4d}{\displaystyle \frac{d^d\mathrm{}}{(2\pi )^d}\frac{1}{[2k\mathrm{}+k^2^2+m^2+iϵ][\mathrm{}^2m^2+iϵ]}}`$ (55)
$`=_{11}^L(s,)+(4\pi )^2\overline{}_{11}(s,),`$ (56)
where $`s=k^2`$. The total integral, without the subscript $`I`$, will have, in addition to the infrared part, a regular part but this can be discarded since it can be expanded with the various terms absorbed into the low energy constants. The divergent IR contribution is

$`_{11}^L(s,)={\displaystyle \frac{(s^2+m^2)}{2s}}\left[2L+{\displaystyle \frac{1}{(4\pi )^2}}\mathrm{ln}{\displaystyle \frac{m^2}{\mu ^2}}\right],`$ (57)
$`\mathrm{where}L={\displaystyle \frac{1}{32\pi ^2}}\left({\displaystyle \frac{2}{d4}}+\gamma 1\mathrm{ln}4\pi \right),`$ (58)
with $`\gamma `$ denoting Euler’s constant. The real part of the finite piece is

$`\mathrm{}\overline{}_{11}(s,)`$ $`=`$ $`\{\begin{array}{cc}\frac{m}{\sqrt{s}}\left[x+2\sqrt{1x^2}\mathrm{cos}^1(x)\right]& (|x|1)\\ \frac{m}{\sqrt{s}}\left[x+\sqrt{x^21}\mathrm{ln}\frac{x+\sqrt{x^21}}{x\sqrt{x^21}}\right]& (|x|1)\end{array}`$ (59)
$``$ $`{\displaystyle \frac{1}{\overline{M}}}J\left({\displaystyle \frac{s^2}{2\overline{M}}}\right),`$ (60)
where $`x=(s^2+m^2)/(2m\sqrt{s})`$. The last line here makes contact with the HB integral $`J`$ in I for which $`\sqrt{s}`$ in the denominator is approximated in leading order by the average baryon mass $`\overline{M}=\frac{1}{2}(M+M_\mathrm{\Delta })`$. The imaginary part of the IR integral is

$$\mathrm{}\overline{}_{11}(s,)=\frac{\pi }{s}\sqrt{(s^2+m^2)^24sm^2}\theta (s[+m]^2),$$
(61)
where $`\theta `$ denotes the Heaviside step function.

In the following we shall implictly give only the real part of the IR integral and the $`iϵ`$ in the integrand denominators will be suppressed. For the integral

$$i\mu ^{4d}_I\frac{d^d\mathrm{}}{(2\pi )^d}\frac{\mathrm{}_\mu }{[(\mathrm{}+k)^2^2][\mathrm{}^2m^2]}=_{11}^{(1)}(s,)k_\mu .$$
(62)
The finite part is $`(4\pi )^2\overline{}_{11}^{(1)}`$, with

$$\overline{}_{11}^{(1)}(s,)=\left[\frac{s^2+m^2}{2s}\right]\overline{}_{11}(s,)\left[\frac{s^2}{2\overline{M}^3}\right]J\left(\frac{s^2}{2\overline{M}}\right),$$
(63)
with the last expression giving the HB reduction.

For the integral

$$i\mu ^{4d}_I\frac{d^d\mathrm{}}{(2\pi )^d}\frac{\mathrm{}_\mu \mathrm{}_\nu }{[(\mathrm{}+k)^2^2][\mathrm{}^2m^2]}=_{11}^{(2)}(s,)g_{\mu \nu }+_{11}^{(3)}(s,)k_\mu k_\nu .$$
(64)
The finite part of the first term is $`(4\pi )^2\overline{}_{11}^{(2)}`$, with

$`\overline{}_{11}^{(2)}(s,)`$ $`=`$ $`\left[{\displaystyle \frac{4m^2s(s^2+m^2)^2}{12s}}\right]\overline{}_{11}(s,){\displaystyle \frac{(s^2)}{36s^2}}\left[6m^2s(s^2)^2\right]`$ (65)
$``$ $`{\displaystyle \frac{1}{3\overline{M}}}[m^2{\displaystyle \frac{(s^2)^2}{4\overline{M}^2}}]J\left({\displaystyle \frac{s^2}{2\overline{M}}}\right){\displaystyle \frac{(s^2)}{18\overline{M}^2}}\left[3m^2{\displaystyle \frac{(s^2)^2}{2\overline{M}^2}}\right],`$ (66)
with the last line giving the HB reduction.

One of the integrals involving one pion and two baryon denominators is

$`i\mu ^{4d}{\displaystyle _I}{\displaystyle \frac{d^d\mathrm{}}{(2\pi )^d}}{\displaystyle \frac{\mathrm{}_\mu \mathrm{}_\nu }{[(\mathrm{}+p^{}+q^{})^2_A^2][(\mathrm{}+p^{})^2_B^2][\mathrm{}^2m^2]}}`$ (67)
$`=_{12}^{(3)}(s,_A,_B)g_{\mu \nu }+\mathrm{}.`$ (68)
The finite part is $`(4\pi )^2\overline{}_{12}^{(3)}`$, with

$`\overline{}_{12}^{(3)}(s,_A,_B)=\left({\displaystyle \frac{}{_A^2}}+{\displaystyle \frac{}{_B^2}}\right){\displaystyle \underset{0}{\overset{1}{}}}𝑑z\overline{}_{11}^{(2)}(\overline{k}^2,\overline{}),`$ (69)
$`\mathrm{where}\overline{k}^2=m^2z^2+(sM^2m^2)z+M^2,`$ (70)
$`\mathrm{and}\overline{}^2=m^2z^2+(_A^2_B^2m^2)z+_B^2.`$ (71)
This was evaluated by performing the differentiation analytically and carrying out the integration numerically using Gaussian quadrature. It is reduced to HB form by dropping terms involving $`m^2`$ from the expression $`\overline{k}^2\overline{}^2+m^2`$ and setting $`\overline{k}^2M^2`$ in denominators. This allows the integral to be performed, giving

$`\overline{}_{12}^{(3)}(s,M_\mathrm{\Delta },M){\displaystyle \frac{1}{6\overline{M}^2}}\left[F_3(\eta \delta ,0)m^2+\frac{2}{3}(\eta \delta )^2\right],`$ (72)
$`\mathrm{where}F_3(\alpha ,\beta )={\displaystyle \frac{1}{\alpha \beta }}\left[(\alpha ^2m^2)J(\alpha )(\beta ^2m^2)J(\beta )\right],`$ (73)
and we recall the definitions of I: $`2\overline{M}\eta =sM^2`$ and $`\delta =M_\mathrm{\Delta }M`$.

The integral involving one pion and three baryon denominators is

$`i\mu ^{4d}{\displaystyle _I}{\displaystyle \frac{d^d\mathrm{}}{(2\pi )^d}}{\displaystyle \frac{\mathrm{}_\mu \mathrm{}_\nu }{[(\mathrm{}+p^{})^2_A^2][(\mathrm{}+p^{}+q^{})^2_B^2][(\mathrm{}+p)^2_C^2][\mathrm{}^2m^2]}}`$ (74)
$`=_{13}^{(3)}(s,t,_A,_B,_C)g_{\mu \nu }+\mathrm{}.`$ (75)
The finite part is $`(4\pi )^2\overline{}_{13}^{(3)}`$, with

$`\overline{}_{13}^{(3)}(s,t,_A,_B,_C)=\left({\displaystyle \frac{}{_A^2}}+{\displaystyle \frac{}{_B^2}}+{\displaystyle \frac{}{_C^2}}\right)^2{\displaystyle \underset{0}{\overset{1}{}}}𝑑xx{\displaystyle \underset{0}{\overset{1}{}}}𝑑y\overline{}_{11}^{(2)}(\overline{k}^2,\overline{})`$ (76)
$`\mathrm{where}\overline{k}^2=\left[sM^2m^2+m^2x(1y)\right]x(1y)txy(1x)+M^2,`$ (77)
$`\mathrm{and}\overline{}^2=m^2x^2(1y)^2+\left[_A^2_B^2+m^2t(1x)\right]xy`$ (78)
$`+(_B^2_C^2m^2)x+_C^2.`$ (79)
This was evaluated with a combination of numerical differentiation and Gaussian numerical integration. The HB reduction is carried out by making similar approximations to those mentioned above and this allows both integrals to be performed. For the cases mentioned in the text this gives

$`\overline{}_{13}^{(3)}(s,t,M_\mathrm{\Delta },M,M_\mathrm{\Delta }){\displaystyle \frac{1}{12\overline{M}^3}}\left[F_1(\eta ,\delta )+\frac{2}{3}(\eta 2\delta )\right],`$ (80)
$`\mathrm{where}F_1(\alpha ,\beta )={\displaystyle \frac{(\alpha ^2m^2)J(\alpha )+(2\beta ^23\alpha \beta +m^2)J(\beta )+(\alpha \beta )(2\beta ^2m^2)}{(\alpha \beta )^2}},`$ (81)
and

$`\overline{}_{13}^{(3)}(s,t,M_\mathrm{\Delta },M,M){\displaystyle \frac{1}{12\overline{M}^3}}\left[F_2(\eta ,\delta )+\frac{2}{3}(\eta \delta )\right],`$ (82)
$`\mathrm{where}F_2(\alpha ,\beta )={\displaystyle \frac{(\alpha ^2m^2)\beta J(\alpha )(\beta ^2m^2)\alpha J(\beta )+\pi m^3(\alpha \beta )}{\alpha \beta (\alpha \beta )}}.`$ (83)