text:
# Correcting for lensing bias in the Hubble diagram

## 1 Introduction

During the last decade, gravitational lensing has become one of the most important tools in cosmology. Weak lensing measurements are used to obtain information on the amount and distribution of matter in the universe (see, e.g., van Waerbeke et al., 2001; Hennawi et al., 2001). Strong lensing has long been recognized as a tool to constrain cosmological parameters (see, e.g., Koopmans & Fassnacht, 1999; Browne, 2000, for measurements of the Hubble constant, $`H_0`$). Controversial upper limits on the amount of dark energy have been derived from the statistics of multiply imaged quasars and radio sources (Kochaneck, 1996; Falco et al., 1998). Using a different set of Schechter parameters for E/S0 galaxies Chiba & Yoshii (1999) conclude that the same quasar data set is best fitted with a universe dominated by dark energy, i.e., in agreement with the concordance cosmology, $`(\mathrm{\Omega }_\mathrm{M},\mathrm{\Omega }_X)=(0.3,0.7)`$ whereas Keeton (2002) argue that when calibrated to counts of distant galaxies, the method loses most of its sensitivity to the amount of dark energy.

Further, in Goobar et al. (2002b) it is shown how one could use multiple images of supernovae (SNe) to simultaneously measure $`H_0`$, the fractional energy density components and the equation of state parameter of a possible dark energy component. The classical SN cosmological test however, is to study the distance-redshift relation using Type Ia SNe as standard candles. The effect of gravitational lensing on these measurements is to cause an additional dispersion in the observed magnitudes and thus be a source of systematic error in the cosmological parameter determination (see e.g. Wambsganss et al., 1997; Mörtsell et al., 2001b). In this note we show how the negligence of this effect may lead to a systematic underestimation of $`\mathrm{\Omega }_\mathrm{M}`$. We also investigate how this bias can be diagnosed and corrected for by using the correct probability distribution function (pdf) for the dispersion in SN magnitudes.

In Sects. 2 and 3 we discuss the method and the numerical simulations used in the analysis. In Sect. 4, we investigate what kind of bias we expect if not taking lensing effects into consideration, i.e., assuming a Gaussian magnitude distribution. In Sect. 6, the dispersion in magnitudes due to gravitational lensing is parametrized for different matter distributions and redshifts and in Sect. 7 we correct for lensing effects by performing a maximum likelihood (ml) analysis using the correct pdf’s.

## 2 The distance-redshift cosmological test

The distances and redshifts of Type Ia SNe can be used to fit, e.g., the mass energy density, $`\mathrm{\Omega }_\mathrm{M}`$, and the dark energy density, $`\mathrm{\Omega }_X`$, in a Friedmann-Lemaître universe. The relation between the parameters can be expressed as

$`m(\theta ,,z)`$ $`=`$ $`+5\mathrm{log}_{10}\left[d_L^{}(\theta ,z)\right],`$ (1)
$``$ $`=`$ $`25+M+5\mathrm{log}_{10}(c/H_0),`$ (2)
where $``$ is a nuisance parameter containing the absolute magnitude, $`M`$, of the SNe and $`H_0`$ is the Hubble constant. Here $`\theta `$ represents the parameter vector $`\theta =(\mathrm{\Omega }_\mathrm{M},\mathrm{\Omega }_X,w)`$. Further the luminosity distance $`d_L^{}`$ is given by

$`d_L^{}`$ $`=`$ $`\{\begin{array}{cc}(1+z)\frac{1}{\sqrt{\mathrm{\Omega }_k}}\mathrm{sin}(\sqrt{\mathrm{\Omega }_k}I),\hfill & \mathrm{\Omega }_k<0\hfill \\ (1+z)I,\hfill & \mathrm{\Omega }_k=0\hfill \\ (1+z)\frac{1}{\sqrt{\mathrm{\Omega }_k}}\mathrm{sinh}(\sqrt{\mathrm{\Omega }_k}I),\hfill & \mathrm{\Omega }_k>0\hfill \end{array}`$ (6)
$`\mathrm{\Omega }_k`$ $`=`$ $`1\mathrm{\Omega }_\mathrm{M}\mathrm{\Omega }_X,`$ (7)
$`I`$ $`=`$ $`{\displaystyle _0^z}{\displaystyle \frac{dz^{}}{H^{}(z^{})}},`$ (8)
$`H^{}(z)`$ $`=`$ $`\sqrt{(1+z)^3\mathrm{\Omega }_\mathrm{M}+\mathrm{\Omega }_X(1+z)^{3(w+1)}+(1+z)^2\mathrm{\Omega }_k}.`$
The usual way to proceed is to make a $`\chi ^2`$-fit of $`(\mathrm{\Omega }_\mathrm{M},\mathrm{\Omega }_X,)`$ using Eq. (1) assuming that the dispersion in magnitudes is Gaussian distributed (as the intrinsic dispersion of Type Ia SNe seems to be). The same procedure can be used to fit, e.g., the equation of state parameter, $`w`$, of a more general dark energy component (Goliath et al., 2001).

Gravitational lensing will however induce an asymmetry in the magnitude distribution and in order to avoid a bias in the cosmological parameter estimation, we need to perform a ml-analysis with the correct pdf.

## 3 Simulated data sets

We use the numerical simulation package, SNOC, the SuperNova Observation Calculator (Goobar et al., 2002a), to obtain simulated samples of the intrinsic dispersion and gravitational lensing effects of Type Ia SNe. The intrinsic dispersion and measurement error is represented by a Gaussian distribution with $`\sigma _m=0.16`$ mag. Gravitational lensing effects are calculated by tracing the light between the source and the observer by sending it through a series of spherical cells in which the dark matter distribution can be specified.

We model compact dark matter as point masses and smooth dark matter with the Navarro-Frenk-White \[nfw; Navarro et al. (1997)\] density profile using halo parameters, mass distributions and number densities as outlined in Bergström et al. (2000). Note that the results obtained in this paper are not sensitive to the exact parameterization of the smooth halo profile. Also, the results are independent of the individual masses of the compact objects as well as their clustering properties on galaxy scales (Holz & Wald, 1998; Bergström et al., 2000). Any eventual small subhalo structure in the dark matter halos does not act as a compact component (Mörtsell et al., 2001a).

Two different types of simulated data sets have been used for the analysis presented in this note. The first one, A, assumes that the number of observed Type Ia SNe is constant per comoving volume, which means that the number of events increases rapidly for higher redshifts. The B distribution instead assumes a constant number of SNe per redshift interval, i.e., a uniform $`z`$-distribution. Each data set consists of 2000 SNe in the redshift interval $`0.01<z<2.0`$.

The number 2000 corresponds roughly to one year’s data from the proposed satellite telescope Supernova/Acceleration Probe \[snap; Perlmutter et al. (2000)\]. The exact redshift distribution of the Type Ia SNe to be followed by snap is the subject of several ongoing science and instrumental optimization studies. The current anticipated distribution is roughly a constant rate per comoving volume for $`z<1`$ and a uniform distribution at higher redshifts. In that respect is the A distribution to be regarded as a limiting case of maximal lensing effects due to the large number of high redshift SNe.

To each of these data sets, 300 SNe between $`0.04<z<0.08`$ have been added. These numbers correspond to the predicted results from the Supernova Factory Campaign (Aldering, Supernova Factory Webpage, ) that is scheduled to start in 2003. Since lensing effects get larger at higher redshifts, the magnitudes of theese SNe are unaffected. For example, for a supernova at redshift $`z=0.08`$ the lensing effects are in general of the order of 0.01 magnitudes for a dark matter model where the fraction of compact objects is assumed to be 20 %.

All simulations have been made assuming a universe with $`\mathrm{\Omega }_\mathrm{M}=0.3`$, $`\mathrm{\Omega }_X=0.7`$, $`w=1`$ and $`H_0=65`$ km s$`^1`$Mpc-1.

## 4 The $`\chi ^2`$-fit of cosmological parameters

Making a $`\chi ^2`$-fit of $`(\mathrm{\Omega }_\mathrm{M},\mathrm{\Omega }_X,)`$ using Eq. (1) (assuming $`w=1`$) will for the data set B, described in Sect. 3, result in the solid contour presented in Fig. 1(a) if lensing effects are absent.

If lensing dispersion is added, the fits of the cosmological parameters will be biased in accordance with the figure. All contours show the $`1.51\sigma `$ level, which corresponds to a $`68\%`$ confidence region of including the true value of the two parameters.

The amount of bias depends on the matter distribution responsible for the lensing effects. The dashed contour shows the result if only a nfw halo model is considered, while the dotted (dash-dotted) contour is based on a simulation where the fraction of compact objects in the dark matter model is 20 % (40 %).

Fig. 1(b) shows similar fits for the $`(\mathrm{\Omega }_\mathrm{M},w)`$-plane where a flat universe is assumed $`(\mathrm{\Omega }_X=1\mathrm{\Omega }_\mathrm{M})`$. $``$ is still treated as a nuisance parameter, i.e., no prior knowledge is assumed. Also this figure is based on the data set B, but using the A distribution does not significantly alter the qualitative nature of the results. Gravitational lensing starts to become important for $`z>1.0`$, and in the interval $`1.0<z<2.0`$, the distributions A and B are quite similiar.

In Fig. 1(b) it is evident how important it is to consider lensing effects since a naive analysis of lensed data may in fact misleadingly exclude a cosmological constant as the explanation for the dark energy.

The size of the bias, i.e., the systematic errors due to lensing in the above figures depends on the asymmetry of the magnitude distribution, which is described in Sect. 5. This in turn depends on the redshift for a given dark matter model. The statistical error – the sizes of the ellipses – on the other hand only depend on the sample size so it is interesting to see where these two errors are comparable.

In Fig. 2 the solid curve is the estimated statistical error of $`\mathrm{\Omega }_\mathrm{M}`$ plotted as a function of the sample size. It should be noted here that zero events means zero events from the high-$`z`$ distribution, whereas the 300 events from the Supernova Factory simulations are always present.

It is interesting to note in this figure that the full data set of 2300 SNe is not necessary to get systematic errors comparable to the statistical error, but gravitational lensing has to be considered also for rather small data samples containing high redshift SNe ($`z>1`$) if the fraction of compact objects is non-negligible.

## 5 Determining the lensing dispersion

As we have seen, gravitational lensing may induce sizeable systematic effects when trying to determine cosmological parameters using the distance-redshift relation for standard candle sources. A virtue of gravitational lensing is however that the distribution of luminosities can be used to obtain information on the matter distribution in the Universe, e.g., to determine the fraction of compact objects like primordial black holes or MACHOs.

As shown in Mörtsell et al. (2001a), gravitational lensing effects are quite similar for different smooth dark matter halo distributions, but very sensitive to the fraction of the matter density in compact objects, $`f_\mathrm{p}`$. Assuming that the major part of the matter density in the universe is in either smooth halos or compact objects, we can thus parametrize lensing effects with $`f_\mathrm{p}`$.

Using SNOC, large data sets of synthetic SN observations over a broad redshift range were created with a variable fraction of compact objects ranging from 0 to 40 %, see Fig. 3.

Using the Kolmogorov-Smirnov test, it was shown how the matter density in compact objects can be determined using SN data comparable to set A or B to $`<5`$ % accuracy which we will assume to be the predicted uncertainty in the parameter $`f_\mathrm{p}`$, see also Metcalf & Silk (1999) and Seljak & Holz (1999). We now want to parametrize the magnitude distributions obtained from the Monte-Carlo simulations to obtain approximately correct pdf’s to be used in the ml-analysis.

## 6 Parametrizing the pdf

We parametrize the pdf’s for different fractions of compact objects and redshifts with a Gaussian with a high magnification tail, i.e.,

$`f(m)`$ $`=`$ $`ae^{\frac{(mm_0)^2}{2\sigma ^2}},m>m_\mathrm{c}`$
$`f(m)`$ $`=`$ $`ae^{\frac{(mm_0)^2}{2\sigma ^2}}+b|m|10^{sm},m<m_c.`$ (9)
This distribution is motivated by the fact that the dominating intrinsic Gaussian dispersion will be shifted toward fainter values since a majority of lines-of-sight contain less matter than in a Friedmann-Lemaître universe. The tail represents the effect of, e.g., galaxies lying close to some lines-of-sight causing high magnifications.

We can set $`s=1.5`$ and $`m_c=0`$ for all cases and still have a reasonable $`\chi ^2`$. Since the fitted probability density functions are normalized to give an integrated probability of unity before performing the maximum likelihood analysis, our parameters are $`m_0,\sigma `$ and $`b/a`$.

As an example, see Fig. 4, where Eq. (6) is fitted to a reference sample corresponding to 40 000 sources at $`z=1`$ in a universe with a 15 % mass fraction in compact objects. Here, $`m_0=0.01756,\sigma =0.1653`$ and $`b/a=0.02530`$ with $`\chi ^2/\mathrm{ndf}=86.44/68`$.

For redshifts $`z<0.5`$ (regardless of the value of $`f_\mathrm{p}`$) and $`0.5<z<1`$ with low $`f_\mathrm{p}`$, one may ignore the high magnification tail (using a simple Gaussian), and still have acceptable $`\chi ^2`$. Note however that the mean is still shifted from the zero value. For higher $`z`$ and $`f_\mathrm{p}`$, we need the full pdf. E.g., for $`z=2`$ and 20 % compact objects, we get $`\chi ^2/\mathrm{ndf}10`$ for a Gaussian and $`\chi ^2/\mathrm{ndf}1`$ including the high magnification tail.

## 7 Correcting for the bias

There is apparently a great need to be able to handle gravitational lensing effects in the fits of the cosmological parameters. If the magnitude distribution presented in Sect. 6 is used in a ml-analysis the results presented in Fig. 5 are obtained. The contours represent the $`68\%`$ confidence region for one specific experimental realisation.111A different realisation would give contours of the same size but these may not have the same centres. They should however also include the true value which gives a limit on the how much of a difference there could be between different realisations. The solid contours are as previously the

confidence region when no lensing effects have been considered while the dashed curve is the $`\chi ^2`$-fit of for a dark matter model of 20 % compact objects (cf. Fig. 1). The dash-dotted curves show the ml-fits using the correct pdf, which is reducing the bias significantly as compared with the simple $`\chi ^2`$-fit. The dotted curves show the result of an ml-fit when fractions of 15 % and 25 % compact objects are assumed when the true value is 20 %. These contours represent the expected systematic bias when under- or overestimating the fraction of dark matter in compact objects with 5 %. The results are a big improvement over the $`\chi ^2`$-analysis. Thus, it should be possible to reduce the bias due to lensing in a fairly effective way, even with the simple approach used in this note.

## 8 Summary and discussion

Neglecting gravitational lensing effects may cause systematic errors when determining cosmological parameters using the distance-redshift relation for Type Ia SNe. E.g., if the universe contains a large fraction of compact objects, the matter density $`\mathrm{\Omega }_\mathrm{M}`$ can be severely underestimated or a cosmological constant may be wrongly excluded if one assumes a Gaussian distribution of magnitudes. Using an approximately correct pdf for the SNe magnitudes, this bias can be reduced significantly.

In this note, we show that lensing effects can be parametrized in a simple way by the fraction of compact objects, $`f_\mathrm{p}`$, and the redshift. In general, lensing effects are bigger for larger values of $`f_\mathrm{p}`$ and higher redshifts. If $`f_\mathrm{p}`$ can be determined with $`<5`$ % accuracy, the bias in the cosmological parameter determination can be effectively reduced as compared to the case where the magnitudes are assumed to be Gaussian distributed. For future surveys, e.g., by the proposed snap satellite aiming at doing precision cosmology with Type Ia SNe, this reduction is essential to get systematic errors comparable to the statistical errors.

For the simple analysis performed in this note, we have first determined $`f_\mathrm{p}`$ to be able to choose the correct pdf for the subsequent parameter determination. Another possibility is to make a simultaneous fit of, e.g., $`(f_\mathrm{p},\mathrm{\Omega }_\mathrm{M},\mathrm{\Omega }_X,)`$ which will be the aim of some forthcoming work.

## Acknowledgements

AG is a Royal Swedish Academy Research Fellow supported by a grant from the Knut and Alice Wallenberg Foundation. We would like to thank Greg Aldering for pointing out an assumed error in the redshift range for Supernova Factory.


format: text

text:
# Universal shape ratios for polymers grafted at a flat surface

## 1 Introduction

It has been known since 1934 that polymers are not istantaneously spherical in shape . Since then a lot of work has been devoted to quantify the average deviation from spherical symmetric shape in the limit of large degree of polymerization, since this deviation plays an important role in the hydrodynamic behaviour of dilute polymer solutions . There are several observables which quantify the shape of polymers. One considers the (squared) radius-of-gyration tensor and its eigenvalues $`q_1,\mathrm{},q_d`$, which give the square lengths of the polymer along the principal axes of inertia . The ratios $`q_i/q_j`$ or, alternatively, $`q_i/q_j`$ are universal and may be used to give a universal characterization of the polymer shape. Alternatively, one may use the shape asphericity that allows to distinguish between spherical and rodlike configurations. The equivalent ellipsoid has been experimentally determined in a small-angle x-ray scattering experiment for globular proteins .

The shape of Gaussian polymers has been extensively studied. The shape asphericity, which in its two variants is proportional to the istantaneous variance of the average of eigenvalues of the radius-of-gyration tensor, has been exactly computed . Instead, the ratios of the eigenvalues $`q_i`$ are not exactly known. There are perturbative estimates at second order in $`1/d`$ , which are, however, not very predictive in the physically relevant case $`d=3`$, and several estimates from numerical simulations . The numerical estimates are in good agreement with each other.

In the case of polymers with excluded-volume interactions, asphericity has been computed to second order in $`ϵ=4d`$ in Ref.  and by numerical simulations in Refs. , while the ratios of the eigenvalues of the radius-of-gyration tensor have been computed by Monte Carlo simulations in Refs. . Here we consider the case of polymers with excluded-volume interactions grafted at an impenetrable interface at the ordinary transition . Polymers are modeled as $`N`$-step self-avoiding walks on a simple cubic lattice. We perform an extensive Monte Carlo simulation using a variant of the pivot algorithm , the cut-and-permute algorithm, which has been shown to be optimal for non-adsorbed polymers grafted at a surface . The efficiency of the algorithm allows us to simulate long walks with $`N32000`$ and to study how the presence of the surface changes the shape of the polymer.

We find that the main effect of the presence of the surface is to increase the length of the longest axis of the inertia ellipsoid with respect to the other two, whose relative length remains, instead, of the same order. Nevertheless, being this increment in asphericity only of a few percents, we expect that it can only be detected in high-resolution experiments. The paper is organized as follows. In Sec. 2 we define the observables that will be computed. In Sec. 3 we present our Monte Carlo results, paying due attention to the corrections to scaling. In Sec. 4 we compare our results with those for polymers in the bulk.

## 2 Definitions

We model the polymer as an $`N`$-step self-avoiding walk on a simple cubic lattice, each step representing a single monomer. The polymer has one of the two endpoints grafted at an impenetrable flat surface of equation $`z=0`$, so that the accessible space is $`_+^3=\{𝐫=(𝐫_{},z):𝐫_{}^2,z0\}`$. Let us denote by $`\omega =\{\omega _0,\mathrm{},\omega _N\}`$ a generic polymer configuration, where $`\omega _i=(\omega _i^1,\omega _i^2,\omega _i^3)_+^3`$ is the position vector of the $`i`$-th monomer and $`\omega _0^3=0`$, and let us call $`W_N`$ the set of allowed configurations.

One can introduce several invariant quantities which describe the shape of the polymer. In generic dimension $`d`$ we define the radius-of-gyration or shape tensor

$$Q_{\alpha \beta }=\frac{1}{N}\underset{i>j}{}(\omega _i^\alpha \omega _j^\alpha )(\omega _i^\beta \omega _j^\beta ),$$
(1)
whose trace concides with the squared radius of gyration

$$R_g^2=\mathrm{Tr}Q.$$
(2)
If we denote with $`\{q_1,q_2,\mathrm{},q_d\}`$ the set of eigenvalues of $`Q`$, ordered in such a way that $`q_1q_2\mathrm{}q_d`$, and with $`\overline{q}=\frac{1}{d}\mathrm{Tr}Q`$ their average, we can define the shape asphericity of the polymer as

$$A_d(\omega )=\frac{1}{d(d1)}\underset{\alpha =1}{\overset{d}{}}\frac{(q_\alpha \overline{q})^2}{\overline{q}^2}=\frac{d}{d1}\frac{\mathrm{Tr}(\widehat{Q}_d^2)}{(\mathrm{Tr}Q)^2},$$
(3)
where $`\widehat{Q}=Q\frac{1}{d}(\mathrm{Tr}Q)I_d`$ and $`I_d`$ is the identity matrix. Since the eigenvalues of $`Q`$ correspond to the squared lengths of the equivalent ellipsoid axes and the eigenvectors represent the axes directions, $`A_d(\omega )=0`$ for a spherically symmetric configuration—in this case all eigenvalues take the same value—while $`A_d(\omega )=1`$ if the polymer is completely elongated and all eigenvalues except one vanish. The mean asphericity is obtained by averaging $`A_d(\omega )`$ over the allowed polymer conformations

$$A_d^N=\frac{1}{c_N}\underset{\omega W_N}{}A_d(\omega ),$$
(4)
where $`c_N`$ is the number of allowed configurations and behaves as $`c_N\mu ^NN^{\gamma _11}`$ in the critical limit $`N\mathrm{}`$. Here $`\mu `$ is the critical fugacity, which does not depend on the presence of the surface, and $`\gamma _1`$ is the susceptibility critical surface exponent $`\gamma _1=0.679(2)`$ (Ref. ), $`\gamma _10.680`$ (Ref. ). Other invariant quantities are often considered in order to characterize the polymer shape. For instance, one may consider

$$\widehat{A}_d^N=\frac{d}{d1}\frac{\mathrm{Tr}(\widehat{Q}_d^2)}{(\mathrm{Tr}Q)^2},$$
(5)
or the relative magnitudes of the averages of the eigenvalues, $`q_2/q_1`$, $`q_3/q_1`$. The advantage of studying these quantities is that they do not require computing an average of the ratio of highly fluctuating quantities.

In the limit of large degree of polymerization $`N\mathrm{}`$, these quantities reach universal values with corrections of order $`N^\mathrm{\Delta }`$,

$`A_d^N`$ $`=`$ $`A_d\left(1+{\displaystyle \frac{B}{N^\mathrm{\Delta }}}+\mathrm{}\right),`$ (6)
$`\widehat{A}_d^N`$ $`=`$ $`\widehat{A}_d\left(1+{\displaystyle \frac{\widehat{B}}{N^\mathrm{\Delta }}}+\mathrm{}\right),`$ (7)
$`{\displaystyle \frac{q_2^N}{q_1^N}}`$ $`=`$ $`L_{2,1}^2\left(1+{\displaystyle \frac{l_{2,1}}{N^\mathrm{\Delta }}}+\mathrm{}\right),`$ (8)
$`{\displaystyle \frac{q_3^N}{q_1^N}}`$ $`=`$ $`L_{3,1}^2\left(1+{\displaystyle \frac{l_{3,1}}{N^\mathrm{\Delta }}}+\mathrm{}\right),`$ (9)
where $`\mathrm{\Delta }`$ is the leading correction-to-scaling exponent. On the other hand, the eigenvalues of the gyration tensor scale as

$$q_i^N=b_iN^{2\nu }\left(1+\frac{c_i}{N^\mathrm{\Delta }}+\mathrm{}\right)i=1,2,3.$$
(10)
In the following section we present a Monte Carlo study of the universal behaviour of these quantities.

## 3 Monte Carlo study of the shape asphericity

We computed the eigenvalues of the shape tensor and the asphericities by Monte Carlo simulations. The system was modeled as a self-avoiding walk grafted at the origin of the half cubic lattice $`_+^3`$. In order to sample the phase space $`W_N`$, we used a nonlocal Monte Carlo algorithm, the cut-and-permute algorithm . It is an optimal algorithm for studying polymers grafted at an impenetrable surface at the ordinary transition, since autocorrelation times in CPU units for global observables are simply proportional to the degree of polymerization $`N`$. The algorithm has, therefore, the same optimal behaviour as the pivot algorithm for polymers in the bulk.

The system in the bulk is known to have strong corrections to scaling , since the leading correction-to-scaling exponent is $`\mathrm{\Delta }_{\mathrm{bulk}}0.5`$. It is expected that polymers in the presence of a surface show analogous corrections, and, indeed, our estimates of the leading correction-to-scaling exponent are close to $`0.5`$.

Since the presence of large corrections to scaling introduces systematic errors in the estimates for finite $`N`$, we studied the system for large $`N`$, $`1000N32000`$. At first, we studied the behaviour of the eigenvalues of the shape tensor and of its trace $`R_g^2`$. The results are displayed in Table 1. These quantities are expected to increase as $`N^{2\nu }`$, being $`\nu `$ the same critical exponent as in the bulk. In order to check this prediction, we performed fits of the form

$`q_i_N`$ $`=`$ $`b_iN^{2\nu }i=1,2,3,`$
$`R_g^2`$ $`=`$ $`bN^{2\nu },`$ (11)
where $`b`$, $`b_i`$, and $`\nu `$ are free parameters. In Eq. (11) corrections to scaling are neglected. We expect, therefore, a systematic error in our estimates. In order to take into account the effect of corrections to scaling, we have performed different fits, considering only data with $`NN_{\mathrm{min}}`$. The effective exponents which are obtained decrease with increasing $`N_{\mathrm{min}}`$ for all values of $`N_{\mathrm{min}}`$ considered, so that we can only give an upper bound for the exponent $`\nu `$. The estimates which are less affected by corrections to scaling appear to be those corresponding to the largest eigenvalue. Therefore, taking as an upper bound for $`\nu `$ the estimate from the fit of $`q_3`$ with $`N_{\mathrm{min}}=16000`$, we obtain $`\nu 0.58774(35)`$, which agrees with the best estimates obtained so far for polymers in the bulk, which are $`\nu =0.58758(7)`$ (Ref. ), $`\nu =0.5877(6)`$ (Ref. ).

Let us consider now the ratios of eigenvalues. For such invariant quantities we performed the following three-parameter fits

$$\frac{q_k_N}{q_1_N}=L_{k,1}^2+D\left(\frac{N_0}{N}\right)^\mathrm{\Delta }k=1,2,$$
(12)
where $`L_{k,1}^2`$, $`D`$, and $`\mathrm{\Delta }`$ are free parameters and $`N_0=16000`$ is a fixed parameter introduced in order to reduce the cross-correlations between $`D`$ and $`\mathrm{\Delta }`$ in the correction-to-scaling term. Again, since we are neglecting higher-order corrections to scaling, the estimates of the universal amplitudes are affected by systematic errors. In order to take into account this effect, we performed fits in which only data with $`NN_{\mathrm{min}}`$ were considered and we studied how our estimates approach the asymptotic value as $`N_{\mathrm{min}}`$ increases. The results of the fits are shown in Table 2. In this case the estimates from the three-parameter fits appear to be very stable, so that we can choose as best estimates of the universal ratios the ones for $`N_{\mathrm{min}}=1000`$: $`L_{2,1}^2=2.9272(28)`$, $`L_{3,1}^2=14.6293(66)`$.

Nevertheless, one can notice that, while the correction-to-scaling exponent determined from the ratio $`q_3_N/q_1_N`$ is $`\mathrm{\Delta }=0.450(46)`$, which is in agreement with the estimates for the leading correction-to-scaling exponent for polymers in the bulk , the estimate of $`\mathrm{\Delta }`$ obtained in fits of $`q_2_N/q_1_N`$ is considerably smaller. This can be due to the presence of corrections with opposite sign and could affect our final estimate. For this reason, we also performed fits in which the leading correction-to-scaling exponent was fixed at $`\mathrm{\Delta }=0.5`$, and higher-order corrections were neglected. The estimates obtained are shown in the last three rows of Table 2. Estimates for $`L_{2,1}^2`$ are stable for $`N_{\mathrm{min}}=4000`$, while those for $`L_{3,1}^2`$ are already stable for $`N_{\mathrm{min}}=1000`$. In all cases the $`\chi ^2`$ per degree of freedom is less than 1. Therefore, we can take as final estimates

$`L_{2,1}^2`$ $`=`$ $`2.94931(99),`$ (13)
$`L_{3,1}^2`$ $`=`$ $`14.6449(60).`$ (14)
One can notice that the estimate in Eq. (14) agrees with the one obtained without fixing $`\mathrm{\Delta }`$, while $`L_{2,1}^2`$ in Eq. (13) is not compatible with the previous estimate. Since we know a priori that $`\mathrm{\Delta }`$ should be approximately $`1/2`$, we are inclined to consider (14) a more reliable estimate.

We performed the same kind of analysis for the asphericities, performing fits of the form

$`A_3^N`$ $`=`$ $`A_3+D\left({\displaystyle \frac{N_0}{N}}\right)^\mathrm{\Delta },`$ (15)
$`\widehat{A}_3^N`$ $`=`$ $`\widehat{A}_3+\widehat{D}\left({\displaystyle \frac{N_0}{N}}\right)^\mathrm{\Delta }.`$ (16)
The results are reported in Table 3. Fits are not very stable and the amplitude appearing in the correction to scaling is very small. For $`N_{\mathrm{min}}=1000`$ we estimated $`D(47\pm 10)\times 10^5`$ and $`\widehat{D}(43\pm 7)\times 10^5`$. As before, we also performed a two-parameter fit in which the leading correction-to-scaling exponent was fixed at $`\mathrm{\Delta }=0.5`$. Since all estimates for $`N_{\mathrm{min}}1000`$ are compatible within error bars, our final estimates are

$`A_3`$ $`=`$ $`0.44520(15),`$ (17)
$`\widehat{A}_3`$ $`=`$ $`0.55949(13).`$ (18)
These estimates are in agreement with those obtained letting $`\mathrm{\Delta }`$ be a free parameter.

## 4 Comparison of the results with the case of polymers in the bulk

Let us compare our results with the ones for polymers in the bulk . An early Monte Carlo simulation  of polymers with length $`N220`$, in which the presence of corrections to scaling was taken into account, gave $`A_3=0.431(2)`$. Since then, more computing power and more efficient algorithms became available and much longer polymers could be studied. In Ref. , self-avoiding walks for a single value of $`N`$ ($`N=963`$) were simulated using the pivot algorithm, obtaining $`A_3=0.43333(8)`$ and $`\widehat{A}_3=0.54722(9)`$ for the shape asphericities. From the data in Ref.  for the shape factors, we can compute the eigenvalue ratios $`L_{3,1}^2=14.1102(52)`$ and $`L_{2,1}^2=2.9631(12)`$. Since only one value of $`N`$ was considered, no analysis of the corrections to scaling could be done. Thus, we expect the errors to be grossly underestimated. Indeed, our results show that, for $`N10^3`$, the systematic error due to the scaling corrections is much larger than the quoted statistical error. In Refs.  even longer walks were considered and data were extrapolated to infinite length with leading correction-to-scaling exponents $`\mathrm{\Delta }=0.5`$ and $`\mathrm{\Delta }=0.46`$. The estimates for the asphericities are: $`A_3=0.4302`$ and $`\widehat{A}_3=0.5450`$ for $`\mathrm{\Delta }=0.5`$; $`A_3=0.4303`$ and $`\widehat{A}_3=0.5452`$ for $`\mathrm{\Delta }=0.46`$. From the asymptotic values for shape factors, we obtain for the eigenvalue ratios: $`L_{3,1}^2=13.89`$, $`L_{2,1}^2=2.934`$ for $`\mathrm{\Delta }=0.5`$; and $`L_{3,1}^2=13.92`$ $`L_{2,1}^2=2.939`$ for $`\mathrm{\Delta }=0.46`$. No error bars on the extrapolated values were reported in Refs. . The asphericity $`A_3`$ has also been computed in an extensive Monte Carlo simulation in Ref. , where walks of length $`N40000`$ were considered. If the data reported there are fitted assuming $`\mathrm{\Delta }=0.5`$, one finds $`A_3=0.43074(57)`$, which agrees within error bars with the estimate for $`N=40000`$, $`A_3=0.4309(17)`$, and with the above-reported estimates of Refs. . The same analysis cannot be performed for the eigenvalue ratios, since the estimates at different $`N`$ are not available. Nevertheless, using Table 2 of Ref. , we obtain $`L_{3,1}^214.05`$, $`L_{2,1}^22.963`$. Both values are slightly higher than the corresponding ones in Refs. .

Comparing our estimates (17), (18), with those reported above, we find that, due to the presence of the surface, $`A_3`$ approximately increases by $`3.5\%`$, while $`\widehat{A}_3`$ by $`2.6\%`$. The change in the two asphericity observables is, therefore, of the same order.

Moreover, while $`L_{2,1}^2`$ does not appreciably change, $`L_{3,1}^2`$ increases by about $`5\%`$. This shows that the change in asphericity is mainly due to the fact that the length of the longest axis of the inertia ellipsoid increases with respect to the other two. This is a signal of the fact that the presence of the surface gives the polymer directionality along the normal direction to the surface, which we expect, on average, to coincide with the direction of the larger axis in the inertia ellipsoid. Nevertheless, the difference is only of a few percents. Therefore, we expect that it can be experimentally detected only in high-precision experiments.


text:
# The dominant X-ray wind in massive star binaries

## 1 Introduction

The violent wind-wind collision in massive star binaries creates a region of high temperature shock-heated plasma, which can contribute to the total system emission at radio, infrared, optical, ultraviolet and X-ray wavelengths. Over the last 20 years, theoretical models have focussed mainly on the dynamics of the stellar winds and the wind collision zone (WCZ), and on the resulting X-ray emission (Pittard P2000 (2000) and references therein). As the X-ray emission is dependent on the physical conditions within the WCZ and on the distribution and properties of the unshocked attenuating wind material, observations provide information on basic parameters of the system (see Stevens et al. S1996 (1996); Zhekov & Skinner ZS2000 (2000); Pittard & Corcoran PC2002 (2002)).

To interpret such observations, it is useful to know which wind dominates the X-ray emission. For instance, if there is a velocity difference between the winds, we would like to be able to predict how hot the broad-band emission will be, or if there are abundance differences, how strong or weak certain lines will be (this is of particular relevance now that grating observations can resolve line profiles).

Simple expressions for the X-ray luminosity from each of the shocked winds were presented in a complex, yet elegant paper, which provides some answers to these questions in terms of fundamental parameters of the system (Usov U1992 (1992)). Their applicability has resulted in their common use in the literature and in observing proposals. However, we have recently discovered inconsistencies between results from these expressions and those determined from numerical models, leading to some confusion on the issue of the dominant X-ray emitting wind. In this letter we reinvestigate previous conclusions in the literature.

## 2 Estimates of $`L_1/L_2`$

Let us define $`L_1`$ as the X-ray luminosity from the shocked wind with the greater momentum flux (i.e. the ‘stronger’ wind), and $`L_2`$ the equivalent from the shocked weaker wind. Analytical estimates in the literature of the ratio of $`L_1/L_2`$ exist for two limiting cases: the radiative limit, where the cooling timescale for the hot gas from both winds is assumed small in comparison to the timescale for flow of this gas out of the system (i.e. $`\chi <1`$, see Sec. 2.2); and the adiabatic limit, where the opposite is true. We will first focus on these limiting cases, before discussing the behaviour of $`L_1/L_2`$ between these limits.

### 2.1 $`L_1/L_2`$ in the radiative limit

In this limit, the entire kinetic energy thermalized by the shocks is immediately radiated (normally with the majority at X-ray energies), and the region of shocked gas is thin. Binary systems with a small orbital separation are favoured, and examples near this limit include V444 Cyg (WN5 + O6; P = 4.2 d; Corcoran et al. C1996 (1996)) and DH Cep (O5.5-6V + O6.5-7V; P = 2.11 d; Penny et al. PGB1997 (1997)).

Let us first consider the region of the WCZ which lies directly between the stars. As there is a momentum balance, the relative kinetic energy fluxes are proportional to the ratio of the wind speeds, so emission from this volume should be dominated by gas from the star with the faster wind. If, for simplicity, we initially assume that the pre-shock wind speeds are spatially invariant, we find that the majority of the thermalized energy in the WCZ occurs close to the line of centres (see Fig. 1). Under these conditions, the faster wind dominates the total emission.

If the pre-shock wind speeds are approximately equal, neither star should significantly dominate the total emission. In such circumstances, one might expect secondary considerations, such as the wind momentum ratio, to become important. Luo et al. (LMM1990 (1990)) noted that if $`v_1v_2`$, the weaker wind should be more luminous, based on the premise that it has, on average, a larger velocity component normal to the shock front: referring to V444 Cyg they argue that, “Because of this, the shocked O6 wind dominates the X-ray emission from the shocked region”. However, we show below that even for very large wind momentum ratios, the stronger wind still accounts for almost half of the total emission. Therefore Luo et al. ’s use of the word “dominates” goes too far, although their basic conclusion that the O6 wind is the majority emitter is consistent with our analysis. Complicating factors, such as wind acceleration, are discussed later.

We can obtain analytical estimates for $`L_1/L_2`$ from the relevant equations in Usov (U1992 (1992)). For the radiative limit, the total X-ray luminosity from the external wind shock is (cf. his Eq. 88) $`L_1=2.3\rho _{\mathrm{}}v_1^3z_{\mathrm{OB}}^2`$, where $`\rho _{\mathrm{}}=\dot{M}_1/4\pi D^2v_1`$, $`z_{\mathrm{OB}}`$ is the distance of the contact discontinuity from the centre of the star with the weaker wind ($`=\sqrt{\eta }D/(1+\sqrt{\eta })`$), $`D`$ is the stellar separation, and $`\eta `$ is the momentum ratio of the winds ($`=\dot{M}_2v_2/\dot{M}_1v_1`$). Usov (U1992 (1992)) parameterizes the luminosity from the internal wind shock as $`L_2=0.22\dot{M}_2v_2^2`$.

With the same parameter values as used by Luo et al. (LMM1990 (1990)) for V444 Cyg, these equations yield $`L_1/L_2=0.59`$ (i.e. that the weaker wind is the dominant X-ray emitter), in agreement with the statement in Luo et al. (LMM1990 (1990)). From Usov (U1992 (1992)) we find that $`L_1/L_2\dot{M}_1v_1^2/\dot{M}_2v_2^2\rho _1v_1^3/\rho _2v_2^3`$, which at the stagnation point is proportional to $`v_1/v_2`$, such that as expected the faster wind will normally be the dominant X-ray emitter. While this agreement is satisfying, an indication of potential problems with the equations in Usov (U1992 (1992)) is revealed by considering the predicted value of $`L_1/L_2`$ for identical winds: in this case we obtain $`L_1/L_2=0.2`$, instead of unity!

We have therefore performed a numerical calculation of the kinetic energy flux normal to the wind shocks to precisely determine $`L_1/L_2`$ for various wind parameters. The position of the contact discontinuity (which is also the position of the wind shocks since the wind collision region is thin) was obtained from integration of Eq. 4 in Stevens et al. (SBP1992 (1992)). We find that the value of $`L_1/L_2`$ primarily depends on the ratio of the wind speeds and to a much lesser extent on the wind momentum ratio. Results are shown in Fig. 1 and in Table 1, where we have defined $`\mathrm{\Xi }`$ to be the fractional wind kinetic power normal to the contact discontinuity i.e. $`L_1=0.5\mathrm{\Xi }_1\dot{M}_1v_1^2`$ for the stronger wind.

For $`\eta =(0.01,\mathrm{\hspace{0.33em}0.1},\mathrm{\hspace{0.33em}1.0})`$, we find that $`\mathrm{\Xi }_1=(0.0042,\mathrm{\hspace{0.33em}0.033},\mathrm{\hspace{0.33em}0.167})`$ and $`\mathrm{\Xi }_2=(0.564,\mathrm{\hspace{0.33em}0.403},\mathrm{\hspace{0.33em}0.167})`$ respectively (the latter values being the analytical limit of 1/6), consistent with the solid angle of the wind collision zone as viewed from the star with the stronger (weaker) wind decreasing (increasing) with decreasing $`\eta `$. The ratio of Usov’s equations yield values for $`L_1/L_2`$ which are too low in comparison to the exact numerical calculation by factors of (1.09, 1.74, 4.88) for $`\eta =(0.01,\mathrm{\hspace{0.33em}0.1},\mathrm{\hspace{0.33em}1})`$ respectively, irrespective of the ratio of $`v_1/v_2`$. Thus Usov’s equations for the radiative limit are somewhat in error, being most accurate for low values of $`\eta `$.

Fig. 1 also shows the gradient of $`\mathrm{\Xi }`$ as a function of the off-axis distance, $`r`$. In Table 2 we list the value of $`r`$ at which $`_{r=0}^r(\mathrm{d}\mathrm{\Xi }/\mathrm{d}r)dr`$ is 50% and 90% of the asymptotic value of $`\mathrm{\Xi }`$. For equal strength winds, 50% and 90% of the maximum wind kinetic power thermalized by the shocks (and hence radiated) occurs within $`r/D=0.38`$ and 0.95 (this off-axis distance is much smaller than the corresponding value for adiabatic winds cf. Luo et al. LMM1990 (1990)). Therefore distortion of the WCZ by orbital motion should not significantly alter our previous conclusions. We can also safely relax our assumption of spatially invariant wind speeds as we only need to consider this ratio between the stars.

Finally, for close binaries we should also consider the potential role of radiative acceleration, inhibition (Stevens & Pollock SP1994 (1994)) and braking (Gayley et al. GOC1997 (1997)) of the winds on $`L_1/L_2`$. As already shown, the ratio of the wind speeds is the critical parameter in such systems, and it seems plausible to expect this to continue to be true in the presence of these 3 effects. In such binaries, it is probable that neither wind will have room to reach terminal speed. Since the stronger wind will have more room to accelerate, we would normally expect it to be faster, and to dominate the emission. However, if there is substantial radiative inhibition plus braking, this will cause the strong wind to be by far the slower and will likely shift the dominant energy generation back to the weaker wind.

### 2.2 $`L_1/L_2`$ in the adiabatic limit

In direct contrast to the radiative limit, colliding wind systems are most likely to be near the adiabatic limit if the stellar separation is large. Perhaps the best example of such a system is the Wolf-Rayet binary WR 140 (WC7 + O4V; P = 2900 d; Williams et al. W1997 (1997)). Luo et al. (LMM1990 (1990)) argue that for a WR+O system in the adiabatic limit “the shocked WR stellar wind dominates the X-ray emission.” Myasnikov & Zhekov (MZ1993 (1993)) arrive at a similar conclusion for their “standard system” which is also near this limit: “the whole of the luminosity is due entirely to the emission of the shocked gas of the WR wind”, being more than an order of magnitude greater than the emission from the weaker wind.111Although Myasnikov & Zhekov (MZ1993 (1993)) actually use different abundances for each wind, gas with typical WN abundances has a similar emissivity as gas with solar abundances, such that this difference is unimportant to their findings.

While the appropriate equations in Usov (U1992 (1992)) are known to be lower limits (due to the omission of line emission), the ratio of $`L_1/L_2`$ should not be overly affected and can once again be calculated. Using the same values for the stellar wind parameters as Myasnikov & Zhekov (MZ1993 (1993)), Eqs. 89 and 95 in Usov (U1992 (1992)) yield $`L_1/L_2=1.36`$, indicating that neither wind is particularly dominant. This is clearly in disagreement with the published statements in Luo et al. (LMM1990 (1990)) and Myasnikov & Zhekov (MZ1993 (1993)).

To investigate this issue we have calculated the X-ray luminosity from numerical simulations of the WCZ near the adiabatic limit. Table 3 summarizes our findings for winds with equal, spatially invariant pre-shock speeds. The general trend is for $`L_1/L_2`$ to increase with decreasing $`\eta `$ (which is opposite to the radiative limit), and for this ratio to become very large for small $`\eta `$. We thus find agreement with the work of Luo et al. (LMM1990 (1990)) and Myasnikov & Zhekov (MZ1993 (1993)), and conclude that the appropriate equations in Usov (U1992 (1992)) are again somewhat in error.

The underlying reason for the trend shown in Table 3 concerns the ratio of the cooling timescale to the flow timescale in each of the winds. Stevens et al. (SBP1992 (1992)) noted that, near the local minimum in the cooling curve (which implies that the pre-shock velocity at the stagnation point is in the range $`900\mathrm{km}\mathrm{s}^1<v<3500\mathrm{km}\mathrm{s}^1`$), this ratio can be approximated as $`\chi v_8^4d_{12}/\dot{M}_7`$, where $`v_8`$ is the wind velocity in units of 1000$`\mathrm{km}\mathrm{s}^1`$, $`d_{12}`$ is the distance from the star to the contact discontinuity in units of $`10^{12}\mathrm{cm}`$, and $`\dot{M}_7`$ is the mass-loss rate of the star in units of $`10^7\mathrm{M}_{}\mathrm{yr}^1`$. The shocked wind is radiative when $`\chi <1`$, and approaches the adiabatic limit when $`\chi 1`$. If we ignore the dependence of $`\chi `$ on $`d`$ (i.e. if we were to assume that the relevant distance appropriate to the flow dynamics is the same for each wind), the ratio of this characteristic cooling parameter for the two winds is $`\chi _1/\chi _2\dot{M}_2v_1^4/\dot{M}_1v_2^4`$.

For the results in Table 3 where $`v_1=v_2`$, $`\chi _1/\chi _2\dot{M}_2/\dot{M}_1=\eta `$. Hence as the value of $`\eta `$ decreases in Table 3, the stronger wind becomes more efficient at emitting X-rays relative to the weaker wind. This is consistent with the fact that although the post-shock density at the stagnation point is the same for both winds, its decline with off-axis distance is faster for the weaker wind (as observed also by Myasnikov & Zhekov MZ1993 (1993)). We note also that within a given distance from the stagnation point, the volume occupied by the shocked stronger wind exceeds that occupied by the shocked weaker wind.

One might expect that $`L_1/L_21`$ when $`\chi _1/\chi _2=1`$, irrespective of the value of $`\eta `$. Again ignoring the dependence of $`\chi `$ on $`d`$, we find that $`\eta =(v_2/v_1)^5`$ for $`\chi _1/\chi _2=1`$. Thus to obtain comparable luminosity from each shocked wind when $`\eta =0.1`$, we require $`v_11.6v_2`$. Since, in reality, mass-loss rates from early-type stars can vary by several orders of magnitude, whereas wind velocities lie typically within the range $`10003000\mathrm{km}\mathrm{s}^1`$ (excluding LBV’s), we expect most systems near the adiabatic limit will have a luminosity dominated by the stronger wind.

### 2.3 $`L_1/L_2`$ in-between these limits

Estimates of the luminosity ratio for systems where the shocked winds are in-between the limiting radiative and adiabatic cases must be done on a case-by-case basis. Somewhere in this region of parameter space, the wind which dominates the X-ray emission must switch over from the one with the faster pre-shock speed (radiative limit) to the one with the slower pre-shock speed (adiabatic limit). In contrast, situations where one wind is clearly radiative and the other is closer to being adiabatic will have their X-ray emission dominated by the former (e.g. in $`\gamma ^2`$ Velorum, where $`\chi _11`$ and $`\chi _2>1`$, the X-ray emission is dominated by the shocked WR wind).

## 3 Summary

In this letter we have re-examined the issue of which wind is the dominant X-ray emitter in colliding wind binaries, following the discovery of some confusion in the existing literature. Our work supports the earlier conclusions in Luo et al. (LMM1990 (1990)) and Myasnikov & Zhekov (MZ1993 (1993)), though is sometimes in disagreement with the relevant analytical equations in Usov (U1992 (1992)).

For systems near the radiative limit (typically short period binaries), we find that the primary influence on $`L_1/L_2`$ is the ratio of the wind speeds, $`v_1/v_2`$, since it controls the energy flux ratio. For $`v_1=v_2`$, $`L_1L_2`$ irrespective of the wind momentum ratio, $`\eta `$. When $`v_1v_2`$, the faster wind normally dominates the luminosity, although it is unlikely to do so by more than a factor of 5. The equations in Usov (U1992 (1992)) predict values for $`L_1/L_2`$ which can be in error by up to a factor of 4.

For systems near the adiabatic limit (e.g. long-period, high eccentricity binaries at apastron), we confirm earlier findings that the stronger wind is typically the dominant X-ray emitter, often by an order of magnitude relative to the weaker wind. This is because the dominant driver of $`L_1/L_2`$ is the ratio of $`\chi _1/\chi _2`$, so that the stronger wind can in general more easily radiatively cool. If we force $`\chi _1/\chi _21`$, we find that both shocked winds contribute roughly equally to the X-ray emission irrespective of the value of $`\eta `$. In contrast, the equations in Usov (U1992 (1992)) yield $`L_1L_2`$ irrespective of the assumed wind parameters.

For systems in-between these limits, we anticipate that the dominant luminosity should generally switch from the faster wind (radiative limit) to the slower wind (adiabatic limit). Where one shocked wind is substantially closer to the radiative limit than that of the other, the X-ray emission will naturally be dominated by the former.

This work provides a basic understanding of the dominant factors controlling the luminosity ratio in colliding wind binaries. Detailed calculations are needed to further investigate the effects of wind acceleration/braking and different wind abundances.

###### Acknowledgements.

We would like to thank Ken Gayley, whose comments led to this work and who also as referee helped to clarify some of the important points. JMP would also like to thank PPARC for the funding of a PDRA position.


text:
# Estimating higher order perturbative coefficients using Borel transform

## Abstract

A new method of estimating higher order perturbative coefficients is discussed. It exploits the rapid, asymptotic growth of perturbative coefficients and the information on the singularities in the complex Borel plane. A comparison with other methods is made in several Quantum Chromodynamics (QCD) expansions.

The ordinary perturbative expansions in weak coupling constant in quantum field theories are generally asymptotic expansions, with their perturbative coefficients growing factorially at large orders . In practice, in many cases the asymptotic growth sets in quite early in perturbation, and the rapid growth of the coefficients becomes apparent already at first few orders.

The asymptotic divergence of the perturbative coefficients implies singularities in the Borel plane. There are three kinds of known singularities, all on the real axis: those instanton-induced, ultraviolet renormalons, and infrared renormalons.

In this paper, we show that this rapid growth of the perturbative coefficients and the information on the singularities in the Borel plane can be turned into a useful tool that allows us to estimate the unknown $`(N+1)`$th order coefficient using the known coefficients up to order $`N`$. The method we are going to present can be most easily understood by working out an explicit example, the double-well potential in quantum mechanics. The tunneling between the two potential wells splits the degenerate perturbative ground states into a parity even, and an odd states, and the average energy $`E(\alpha )`$ of the energies of these two states has the perturbative expansion of the form

$$E(\alpha )=\left[\alpha +\underset{n=1}{\overset{\mathrm{}}{}}a_n\alpha ^{n+1}\right],$$
(1)
where $`\alpha `$ denotes the canonical coupling of the model . The Borel transform $`\stackrel{~}{E}(b)`$ of $`E(\alpha )`$, which has the perturbative expansion

$$\stackrel{~}{E}(b)=\left[1+\underset{n=1}{\overset{\mathrm{}}{}}\frac{a_n}{n!}b^n\right],$$
(2)
is known to have multi-instanton–anti-instanton caused singularities at $`b=2nS_0,(n=1,2,3,\mathrm{})`$, where $`S_0=1/6`$ is the one-instanton action (in units of $`1/\alpha `$).

The nature of the singularities can in principle be determined by doing perturbation in the background of multi-instanton–anti-instanton configurations. The closest singularity to the origin at $`b=1/3`$, which determines the leading large order behavior of the expansion (1), can be shown to have the following form

$$\stackrel{~}{E}(b)=\frac{9}{\pi (13b)^2}\left[1\frac{53}{18}(13b)+O[(13b)^2\mathrm{ln}(13b)]\right]+\mathrm{Analytic}\mathrm{part},$$
(3)
which comes from the instanton–anti-instanton contributions to $`E(\alpha )`$. The “Analytic part” denotes terms that are analytic around the singularity.

We now consider the function $`R(b)(13b)^2\stackrel{~}{E}(b)`$ as introduced in , to control the divergence at the singularity. $`R(b)`$ is bounded at $`b=1/3`$, and has a very soft singularity, a suppressed logarithmic cut. Thus the power expansion of $`R(b)`$ around the origin is expected to better behave than that of $`\stackrel{~}{E}(b)`$.

Ignoring the residual logarithmic cut, we would expect the convergence radius of the former is bounded by the second singularity at $`b=2/3`$, thus effectively to become twice that of the latter. To make the expansion of $`R(b)`$ better, we can take a further step of conformally mapping the singularities except the first one as far away as possible from the origin. This way one can hope to have a smoother $`R`$ in the new plane, and a better behavior of the expansion.

The conformal mappings we consider in this paper are such that all the singularities in the Borel plane except the first one are mapped to the unit circle and the first one to a point within the circle . Such a mapping in this case is given by

$$w=\frac{1\sqrt{13b/2}}{1+\sqrt{13b/2}}$$
(4)
which maps the first singularity to $`w_0=(\sqrt{2}1)/(\sqrt{2}+1)=0.171`$ and all others to the unit circle. Without the residual cut-singularity, $`R`$ would be analytic on the unit disc in the $`w`$-plane.

Now the power expansion of $`R(b(w))`$ up to $`O(w^5)`$ is

$`R`$ $`=`$ $`1+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}r_nw^n`$ (5)
$`=`$ $`1+(16+2.667a_1)w+(120+3.556a_2)w^2+(1896.44+3.161a_3)w^3+`$ (7)
$`(22544+2.107a_4)w^4+(254667.7+1.124a_5)w^5+\mathrm{}`$
$`=`$ $`14w+38.22w^2+84w^3+699.1w^4+2538.1w^5+\mathrm{},`$ (8)
where the coefficient at a given order in (7) was calculated using the exact values of $`a_n`$, which are given in Ref. , up to that order less one, and the coefficients in (8) were obtained with the exact values of $`a_n`$. What is interesting about this expansion is that the $`a_n`$-independent constant term in the coefficient at a given order in (7) is much larger than the exact value of the corresponding $`r_n`$ in (8). This then implies that a good approximate estimate of $`a_n`$ can be obtained by simply putting

$$r_n=0$$
(9)
in (7). The estimated values for $`a_n`$ from this prescription are given in Table I, which shows improving accuracy as the order grows, and the accuracy becomes better than 1% at $`n=5`$.

Obviously, the success of this method relies on the exact value of $`r_n`$ being much smaller than the constant term in the corresponding coefficient in expansion (7). We may understand this feature in the following way. We note first the coefficient $`r_n`$ is a linear combination of $`a_i`$’s ($`in`$), which are large numbers. For the sake of argument, let us ignore for the moment the soft cut-singularity of $`R`$ at $`w=w_0`$. Then $`R`$ is analytic on the unit disk in $`w`$-plane, and so we expect the growth of $`r_n`$ is fundamentally limited. To yield a small number out of large numbers, this then suggests a delicate cancellation occur among the large numbers in the expression for $`r_n`$. However, when $`r_n`$’s are written as in expansion (7), the cancellations are incomplete, yielding the large constant terms. This account shows that the essential ingredients for the success of our prescription are: (i) rapid growth of the perturbative coefficients, and (ii) information on the nature and locations of the Borel plane singularities. One may notice the exact values of $`r_n`$’s in (8) still grow quite rapidly. This may be due to the residual logarithmic cut, and the infinite number of singularities on the unit circle.

With this example in mind, we now give the general prescription for estimating uncalculated higher order coefficients. Let us suppose an amplitude $`A(\alpha )`$ and its Borel transform in a given theory have the perturbative expansion of the form

$$A(\alpha )=\alpha +\underset{n=1}{\overset{\mathrm{}}{}}a_n\alpha ^{n+1}$$
(10)
with $`\alpha `$ denoting the canonical coupling of the theory, and

$$\stackrel{~}{A}(b)=1+\underset{n=1}{\overset{\mathrm{}}{}}a_n\frac{b^n}{n!}.$$
(11)
The overall normalization of the expansion (10) was chosen such that the coefficient of the leading term is unit. We assume the Borel transform $`\stackrel{~}{A}(b)`$ has singularities on the real axis on the $`b`$-plane, and the first two singularities on the positive real axis are at $`b=b_0`$, and at $`b=b_1`$, respectively, and the first singularity on the negative real axis is at $`b=b_0`$. We further assume the nature of the first singularity on the positive real axis is known, and is of the form

$$\stackrel{~}{A}(b)=\frac{C}{(1b/b_0)^{1+\nu }}[1+O(1b/b_0)]+\mathrm{Analytic}\mathrm{part},$$
(12)
with $`C`$ a real constant and the ‘Analytic part’ denoting terms analytic about $`b=b_0`$. We now introduce a conformal mapping

$$w=\frac{\sqrt{1+b/b_0}\sqrt{1b/b_1}}{\sqrt{1+b/b_0}+\sqrt{1b/b_1}},$$
(13)
which maps the first singularity on the positive axis to a point within the unit circle and all others to the circle, and the function

$$R(b)=(1b/b_0)^{1+\nu }\stackrel{~}{A}(b).$$
(14)
We then expand $`R(b(w))`$ in power series in $`w`$-plane:

$`R`$ $`=`$ $`1+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}r_nw^n`$ (15)
$`=`$ $`1+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(p_n+q_na_n)w^n.`$ (16)
Note that $`p_n,q_n`$ depend linearly on the leading coefficients to $`a_n`$ only. The estimated $`n`$th order coefficient is then given by

$$a_n=\frac{p_n}{q_n}.$$
(17)
The fractional error of this estimate is

$$\delta _n=\left|\frac{r_n}{p_n}\right|,$$
(18)
and this scheme would work as long as $`r_n`$ is much smaller than $`p_n`$. Of course, it is impossible to know the error a priori, and the reliability of the estimated coefficients should be decided by other circumstantial informations, for instance, as in the following examples, by the pattern of the known terms in expansion (16), or by comparison with the estimated values from using different methods. Since the expansion (16) is in general not convergent on the whole unit disc due to the residual soft singularity, and because of the singularities on the boundary of the disc, it is possible in principle for $`r_n`$ at a certain order to jump to a large value. In such a case the scheme could fail, but that does not seem to occur, at least, at the orders considered in the following QCD expansions.

We now apply our scheme to the various QCD expansions, estimating the uncalculated next higher order coefficients, and compare the obtained values with the estimates from other methods. Of course, the coefficients are renormalization scheme and scale dependent, but for simplicity, in all of the following examples our consideration will be in the $`\overline{\mathrm{MS}}`$ scheme and at the renormalization scale $`\mu =Q`$, where $`Q`$ denotes the energy scale of the problem in consideration. The canonical coupling in these examples is $`a(Q)\alpha _s(Q)/\pi `$, where $`\alpha _s(Q)`$ is the strong coupling constant.

The Bjorken sum rule: Our first example is the QCD correction $`\mathrm{\Delta }`$ for the Bjorken polarized sum rule, which has the perturbative expansion

$$\mathrm{\Delta }(a)=a+\underset{n=1}{\overset{\mathrm{}}{}}d_na^{n+1}.$$
(19)
Incidently, this correction coincides with that of Gross-Lewellyn Smith (GLS) sum rule at next-leading order, and differs only by the small ‘light-by-light’ contribution at next-next-leading order. The first two coefficients $`d_1,d_2`$ are known . The locations of the first singularities of the Borel transform $`\stackrel{~}{\mathrm{\Delta }}(b)`$ and the parameter $`\nu `$ are given as

$$b_0=\frac{1}{\beta _0},b_1=\frac{2}{\beta _0},b_0=\frac{1}{\beta _0},\nu =(\beta _1/\beta _0\gamma _0)/\beta _0,$$
(20)
where $`\beta _0,\beta _1`$ are the first two coefficients of the QCD $`\beta `$-function:

$$\beta _0=(112/3N_f)/4,\beta _1=(10238N_f/3)/16,$$
(21)
and $`\gamma _0`$ is the one-loop anomalous dimension of the twist-four operator appearing in operator product expansion of the sum rule :

$$\gamma _0=(N_c1/N_c)/3=8/9.$$
(22)
Here $`N_f,N_c`$ denote the number of quark flavors and colors, respectively. These are all we need to estimate the unknown next higher order coefficients.

The result is given in Table II, alongside of the estimates from the principle of minimal sensitivity (PMS) method (and the effective charge (ECH) method) by Kataev and Starshenko . The latter is based on the optimization of the perturbative amplitude over the parameter space of the renormalization scheme and scale dependence. It estimates the uncalculated higher order coefficients by reexpanding the optimized amplitude in terms of the coupling in the particular scheme and scale of interest. The recent, exact, partial calculation of the next-next-next-leading order correction  indicates this method, among others , works particularly well .

The estimates from the two methods are in good qualitative agreement, with some exceptions at $`d_3`$. It is remarkable that these two completely different approaches offer such close estimates. As in PMS case, our method works best at $`N_f=3`$, and the estimates become worse as $`N_f`$ increases. Although the two approaches are in an overall agreement, the differences at $`d_3`$ estimates are large enough to be phenomenologically significant. For instance, the new estimate $`d_3156`$ at $`N_f=3`$, which was used in the recent GLS sum rule analysis , results in a less renormalization scale dependent amplitude than with the PMS estimate $`d_3130`$. It is instructive to see this estimate more closely. The expansion of $`R`$, following (16), at $`N_f=3`$, reads

$`R`$ $`=`$ $`1+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}r_nw^n=1+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(p_n+q_nd_n)w^n`$ (23)
$`=`$ $`1+(3.72+1.185d_1)w+(13.488+0.702d_2)w^2+(43.257+0.277d_3)w^3+\mathrm{}`$ (24)
$`=`$ $`1+0.527w+0.709w^2+\mathrm{}`$ (25)
where the last two terms in (25) are exact. Notice that the first two $`p_n`$’s are large compared to the corresponding $`r_n`$’s, satisfying one of the absolutely necessary conditions for our scheme, and the pattern of the first two known coefficients suggests that $`p_3`$ is also likely to be much larger than $`r_3`$. If we believe that the first two $`r_n`$’s have any indication on $`r_3`$, it seems to be safe to assume $`|r_3|2.0`$. This then leads to $`149d_3163`$.

The estimates of the next higher order coefficients $`d_4`$’s are obtained using the $`d_3`$ estimates. Compared to the PMS estimates there is some significant difference at this order. Since the obtained values depend on the estimated $`d_3`$, one should note that there is another uncertainty arising from the error in $`d_3`$ estimate.

The Adler function: The Adler function $`D(a)`$ for the vector current-current correlation function of different quark flavors has the perturbative expansion

$$D(a)=a+\underset{n=1}{\overset{\mathrm{}}{}}d_na^{n+1}$$
(26)
where the first two coefficients are known . Several physical observables can be related to this function through the dispersion relations, and the perturbative coefficients of the former can be obtained once the corresponding coefficients in the Adler function are known.

The locations of the first renormalon singularities and $`\nu `$ in this case are given by

$$b_0=\frac{2}{\beta _0},b_1=\frac{3}{\beta _0},b_0=\frac{1}{\beta _0},\nu =2\beta _1/\beta _0^2,$$
(27)
and the resulting estimates for the first coefficients are given in Table II, along with the PMS, ECH values. Again, our method works best at $`N_f=3`$ for the first two coefficients, and in this case the predicted $`d_3=24.4`$ is slightly smaller than the PMS value $`d_3=27.5`$.

The quark mass: This last example is concerned with the on-shell quark mass. The on-shell quark mass and the $`\overline{\mathrm{MS}}`$ mass are related by

$$\frac{m_{\mathrm{OS}}}{m_{\overline{\mathrm{MS}}}}=1+\frac{4}{3}\left[a+\underset{n=1}{\overset{\mathrm{}}{}}r_na^{n+1}\right],$$
(28)
and the first two coefficients are known . The relevant parameters in this case are

$$b_0=\frac{1}{2\beta _0},b_1=\frac{3}{2\beta _0},b_0=\frac{1}{\beta _0},\nu =\beta _1/2\beta _0^2,$$
(29)
and the estimated first coefficients are in Table III, along with the values from the residue based method of Pineda  (also, for comparison, those of large $`\beta _0`$ approximation ), which is known to work well in this case. The latter method relies on the rapid convergence of the perturbative calculation of the renormalon residue, and the expansion of the Borel transform about the first renormalon . The agreement between these two approaches at $`r_3`$ estimates is remarkable.

To conclude the paper, we presented a new method of estimating higher order unknown coefficients. It is based on the two generic features of the asymptotic expansions, the rapid growth of the coefficients and the Borel plane singularities. It can provide independent estimates for the yet unknown coefficients, which may prove useful in physical analysis.

###### Acknowledgements.

T.L. is thankful to Gorazd Cvetic for useful discussions and comments. This work was supported by the BK21 Core Program. Postscript: After submission of our paper we were informed that in the published version  of hep-ph/9811367 Borel transform with a conformal mapping was tried on a model function to estimate higher order coefficients. Unlike our method, however, it did not remove the first Borel singularity, and because of that it reached to a conclusion substantially different from ours, that the method can work only for large orders (say, $`N7`$) and is not applicable for low orders (N$``$ 3) relevant for present QCD perturbations.


text:
# Vafa-Witten theorem and Lee-Yang singularities

## Abstract

We prove the analyticity of the finite volume QCD partition function for complex values of the $`\theta `$-vacuum parameter. The absence of singularities different from Lee-Yang zeros only permits $``$ cusp singularities in the vacuum energy density and never $``$ cusps. This fact together with the Vafa-Witten diamagnetic inequality implies the vanishing of the density of Lee-Yang zeros at $`\vartheta =0`$ and has an important consequence: the absence of a first order phase transition at $`\theta =0`$. The result provides a key missing link in the Vafa-Witten proof of parity symmetry conservation in vector-like gauge theories and follows from renormalizability, unitarity, positivity and existence of BPS bounds. Generalizations of this theorem to other physical systems are also discussed, with particular interest focused on the non-linear ℂℙN sigma model.

The Vafa-Witten theorem vw  is one of the very few non-perturbative analytic results of QCD. In absence of supersymmetry, non-Abelian continuum gauge theories in 3+1 dimensions do exhibit many very interesting non-perturbative phenomena which have been so far elusive to analytic approaches. However, by a very simple and elegant argument Vafa and Witten vw  explained why we do not observe spontaneous violation of parity symmetry, not only in QCD but in any gauge theory with Dirac fermions in 3+1 space-time dimensions. The relevance of the physical implications of this theorem require to further clarify some questions raised about the Vafa-Witten argument (see ch  for a review and references therein).

Vafa and Witten claimed that their argument should hold for any parity violating order parameter, but this might not be true in general. In fact, order parameters including higher order composite operators might be non-renormalizable and it is unclear how the argument proceeds in such a case ch ; wu . On the other hand, non-local order parameters like a Wilson loop around a very small compactified space direction acquire non-trivial vev’s. Thus, one should at least restrict the type of order parameter to only local operators which are renormalizable, parity odd and preserve all other symmetries of the system (including gauge symmetry and Lorentz invariance). One of the order parameters matching all these requirements in 3+1 dimensions is the topological density $`ϵ^{\mu \nu \sigma \rho }F_{\mu \nu }F_{\sigma \rho }`$. However, even in that case the argument of Vafa and Witten can be questioned ag ; Ji . The aim of this note is to rigorously prove the Vafa-Witten theorem for this particular order parameter.

In this case there is no doubt about the existence of a well defined free energy for any real value of the $`\theta `$ parameter. This is guaranteed by unitarity of the theory, which is translated into Osterwalder-Schrader (OS) positivity of the Euclidean functional measure of the theory with real $`\theta `$-terms es ; am . We recall that OS positivity is not just the standard positivity used in the Vafa-Witten argument. This special property excludes the possibility that pathological scenarios where the free energy is ill-defined ag  might occur in these theories. The Vafa-Witten argument can be correctly applied to show that the free energy density $`_\theta =E_\theta /V`$ at $`\theta =0`$ is lower or equal than at $`\theta 0`$, i.e. $`_0_\theta `$.

However, this fact does not mean that the theorem is already proven. One needs to exclude the existence of a cusp in the energy density at $`\theta =0`$. The appearance of the cusp would signal the existence of a first order phase transition and spontaneous breaking of parity at $`\theta =0`$ without violation of the Vafa-Witten inequality for free energies. The only extra property one needs to prove is smoothness of $`_\theta `$ at $`\theta =0`$, which would make the existence of such a cusp not possible. In other terms, $`\theta =0`$ is always an absolute minimum of the free energy density but only if $`_\theta `$ is smooth we will have $`_0^{}=0`$. Only in that case we can make sure that the vacuum expectation value of the topological density vanishes and the Vafa-Witten theorem holds for this particular order parameter.

In fact, due to Bragg symmetry $`_{\pi +\theta }=_{\pi \theta }`$, which follows from CP symmetry $`_\theta =_\theta `$ and $`\theta `$-periodicity $`_\theta =_{\theta +2\pi }`$, the two classically CP symmetric points $`\theta =0,\pi `$ are extremals of the vacuum energy density. However, $`_\theta `$ is not always smooth at $`\theta =\pi `$: there are cases where CP is spontaneously broken at that point. In principle, the same pathology could appear also at $`\theta =0`$.

In order to prove smoothness, the role of Lee-Yang singularities becomes very relevant. Because of Bragg symmetry and unitarity, non-analyticities are harder to trace in the pure real $`\theta `$ sector. Non-analyticity in $`\theta `$ can be better inferred from the lack of an analytic continuation into the complex $`\theta `$ plane. Indeed, in that case OS positivity (i.e. unitarity) is not preserved and there is even no guarantee that the free energy is well defined, in which case it will not make sense to talk about smoothness. Essentially, there are two possible ways in which the $`\theta `$ theory might give rise to pathological non-analyticities in the complex sector. One possibility is that the analytic continuation to complex values of $`\theta `$ is not defined at all, i.e. the partition function itself becomes divergent. We will see that this possibility does not occur in gauge theories. The other possibility is that the partition function exhibits a sequence of zeros converging to $`\theta =0`$ in the infinite volume limit. This is the scenario advocated by Lee and Yang ly  which signals in many cases the existence of a phase transition like in the Ising model in presence of an external (imaginary) magnetic field.

The partition function splits into a sum of functional integrals over the different topological sectors

$`𝒵_\theta (g)={\displaystyle \underset{q=\mathrm{}}{\overset{\mathrm{}}{}}}\mathrm{e}^{q\mathrm{Im}\theta +\dot{ı}q\mathrm{Re}\theta }{\displaystyle {}_{c_2(A)=q}{}^{}\delta A\mathrm{e}^{S_{\mathrm{YM}}(q)/2g^2}det(\text{/}D+m)}`$ (1)
In a strict sense, the functional integral is UV divergent but it can be regularized in such a way that its positivity properties are preserved aff1 . We will consider such a geometrical regularization whenever it is required. On the other hand, in the infinite volume limit $`𝒵_\theta (g)`$ vanishes unless the vacuum energy is renormalized to 0. But it is just the dependence on $`\theta `$ of the vacuum energy what we want to analyze. Therefore, we will consider throughout the paper a compact space-time with large but finite volume $`VT<\mathrm{}`$.

Positivity of the Yang-Mills Euclidean measure amm ; aff1  and fermionic determinants vw  implies the existence of a real effective action $`S_{\mathrm{eff}}(g,q)=S_{\mathrm{YM}}(q)/2g^2\mathrm{log}det(\text{/}D+m)`$. On the other hand the Yang-Mills action $`S_{\mathrm{YM}}(q)`$ is bounded below by the BPS bound $`S_{\mathrm{YM}}(q)8\pi ^2|q|`$ in each $`q`$-topological sector. We shall see that the total effective action $`S_{\mathrm{eff}}(g,q)`$ is also bounded from below by a similar linear bound involving the absolute value of the topological charge $`|q|`$. In the regularized theory, the coefficient of this linear bound increases with the UV regulating scale. This growth of the effective action is enhanced in the chiral limit of light quark masses by the contribution of fermionic determinants. In general, from positivity of fermionic determinants and Yang-Mills measure and the BPS bounds we have the inequality

$`\left|𝒵_\theta (g)\right|{\displaystyle \underset{q=\mathrm{}}{\overset{\mathrm{}}{}}}\mathrm{e}^{|q\mathrm{Im}\theta |}{\displaystyle {}_{c_2(A)=q}{}^{}\delta A\mathrm{e}^{S_{\mathrm{eff}}(g,q)}}{\displaystyle \underset{q=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle {}_{c_2(A)=q}{}^{}\delta A\mathrm{e}^{S_{\mathrm{eff}}(\stackrel{~}{g},q)}}=𝒵_0(\stackrel{~}{g})`$ (2)
where $`1/\stackrel{~}{g}^2=1/g^2(1/8\pi ^2)|\mathrm{Im}\theta |`$. This shift of the gauge coupling constant can be considered as a change of the renormalization scale. This shows that the partition function of the original theory with a complex $`\theta `$-term is bounded by the partition function of a similar theory with $`\theta `$$`=0`$ but with a different coupling constant. Since the theory at $`\theta =0`$ is unitary and renormalizable, its partition function $`𝒵_0`$ is finite and from (2) it follows that $`\left|𝒵_\theta \right|<\mathrm{}`$ is also finite for small values of $`\mathrm{Im}\theta `$. If in addition the theory is asymptotically free, then the renormalization of $`g^2`$ in the UV fixed point can absorb any value of $`\mathrm{Im}\theta `$, which implies that the radius of convergence of the sum (1) in the complex $`\theta `$ plane, in fact, becomes infinity. To make the argument more precise one should consider the regularized theory. The analyticity of $`𝒵_\theta `$ also appears naturally in lattice regularizations, even with non-topological regularizations of the $`\theta `$–term. However, to prove of absence of poles the topological charge bound (2) is essential. The bare coupling constant $`g`$ has to be fine tuned according to the renormalization group to yield the appropriate continuum limit. But because of asymptotic freedom this coupling goes to zero as the UV regulator is removed and the shift from $`g`$ to $`\stackrel{~}{g}`$ induced by $`\mathrm{Im}\theta `$ simply implies a change in the effective scale of the continuum theory. However, the proof of the Vafa-Witten theorem only requires a finite radius of convergence of $`𝒵_\theta `$ around $`\theta =0`$. This feature does not require asymptotic freedom, just renormalizability.

We remark that one cannot use the same kind of arguments for other order parameters. For instance, the partition function associated to higher odd powers of the topological density might diverge for any complex value of the perturbation parameter, as they are non-renormalizable. A similar remark applies to the perturbation by a chemical potential term in theories at finite density. In that case the missing property is the lack of a BPS bound for the perturbation.

The only remaining possible source of non-analyticity is the presence of Lee-Yang singularities, i.e. zeros in the partition function $`𝒵_\theta `$ which could prevent the existence of a unique limit of the free energy $`\mathrm{log}𝒵_\theta `$ or its derivatives at $`\theta =0`$.

Because of Bragg symmetry and $`\theta `$–periodicity $`𝒵_\theta `$ is a function of $`\mathrm{sin}^2(\vartheta /2)`$ and its zeros come in 4-plets $`\{\vartheta _n,\vartheta _n,\vartheta _n^{},\vartheta _n^{}\}`$ in the modular strip $`\vartheta (\pi ,\pi ]\times `$) of the complex plane. They can be parametrized by their representatives in the positive quadrant $`\vartheta _n`$ , which cannot lie on the real or imaginary axes for any finite volume.

Only zeros converging to the origin in the thermodynamic limit can give rise to a cusp. Moreover, only an infinity of such zeros can produce the cusp.

The Weierstrass factorization theorem implies that the vacuum energy density

$$_V(\vartheta )=\frac{1}{V}\mathrm{ln}Z_V(\vartheta ),$$
(3)
can be rewritten in the functional form ww ; nw

$$_V(\vartheta )=\frac{2}{V}\mathrm{Re}\underset{n=1}{\overset{\mathrm{}}{}}\rho _n\left\{\mathrm{ln}\left(1\frac{\mathrm{sin}^2(\vartheta /2)}{\mathrm{sin}^2(\vartheta _n/2)}\right)+h_n(\mathrm{sin}^2(\vartheta /2),V)\right\}+g__V(\mathrm{sin}^2(\vartheta /2)),$$
(4)
where $`\vartheta _n`$ denotes the order $`n`$ Lee-Yang zero in the positive quadrant, $`\rho _n`$ its degeneracy and $`h_n`$ and $`g`$ are, respectively, polynomial and analytic functions of $`\mathrm{sin}^2(\vartheta /2)`$. The polynomial functions $`h_n`$ are necessary to compensate the divergent effect of zeros with slow approach to the infinite modulus for a given volume.

In the thermodynamic limit the degeneracy factor $`\rho _n`$ of zeros converging to the origin defines a density of zeros $`\rho (z)`$ as a function of the continuous variable $`z=\vartheta _n/V`$. Then, in the thermodynamic limit

$`(\vartheta )`$ $`=`$ $`\mathrm{\hspace{0.17em}2}\mathrm{Re}{\displaystyle _0^\pi }{\displaystyle _0^{\mathrm{}}}dz_1dz_2\rho (\mathrm{sin}^2(z/2))\left\{\mathrm{ln}\left(1{\displaystyle \frac{\mathrm{sin}^2(\vartheta /2)}{\mathrm{sin}^2(z/2)}}\right)+h(\mathrm{sin}^2(\vartheta /2),z)\right\}`$ (5)
$`+g(\mathrm{sin}^2(\vartheta /2)).`$
Notice that only if $`\rho (0)0`$ the contribution from the logarithm at $`z=0`$ could give rise to a cusp at $`\vartheta =0`$. The rest of the contributions from the integral have continuous derivatives for real $`\vartheta `$ at $`\vartheta =0`$. In fact, the logarithmic term in (5) corresponds to a Coulomb potential, in two dimensions, generated by the charge density $`2\rho (z)`$. The cusp at $`\theta =0`$ is only possible if there is a non-vanishing linear charge density at $`\theta =0`$ evans evans1 , i.e. $`\rho (z)\rho _0\delta (z_1mz_2^\beta )`$ for $`|z|1`$ with $`\rho _00`$. A cusp corresponds to a gap in the potential derivative, which is proportional to the charge density $`2\rho _0`$. Since this charge density is non-positive, the corresponding potential is always repulsive, which means that any cusp is always a $``$ cusp instead of a $``$ cusp. Since the vacuum energy density $`(\vartheta )`$ has a local minimum at $`\theta =0`$ because of the Vafa-Witten inequality such a cusp cannot exist. The theorem follows from two very simple facts, the Vafa-Witten inequality and the repulsive character of the negative charge density defined by $`2\rho (z)`$, which follows from the fact that the partition function $`Z_V(\vartheta )`$ is analytic and has zeros but no poles. The only way of having a a $``$ cusp at $`\vartheta =0`$ is with an attractive charge density which can only only appear from poles in $`Z_V(\vartheta )`$. In fact, this is what happens in 1+1 massless scalar field theories with pseudo-periodic boundary conditions adm  where the partition function is non-analytic in the complex plane. However, in QCD the analyticity of the partition function is based in very simple fundamental properties of the theory: unitarity, renormalizability, positivity and existence of BPS bounds. This makes impossible the existence of a cusp in the vacuum energy density at $`\vartheta =0`$. In summary, the theory cannot undergo a first order phase transition at $`\theta =0`$ with parity symmetry breaking and we have a complete proof of the Vafa-Witten theorem.

From these arguments we cannot, however, exclude the existence of first order phase transitions for other values of $`\theta 0`$. In particular, we cannot discard CP symmetry breaking for $`\theta =\pi `$.

We remark that many other theories share the same properties and therefore the Vafa-Witten theorem can also be extended for all of them. One particularly interesting case is that of ℂℙN non-linear sigma models. Indeed, from a similar analysis one can conclude that there is no spontaneous CP symmetry breaking in ℂℙN models at $`\theta =0`$. This result can be checked for the ℂℙ1 model where the exact solution is known for $`\theta =0`$ zamm ; pw .

The solution of the ℂℙ1 model at $`\theta =0`$ describes the interaction of a massive scalar particle in the adjoint representation of $`SU(2)`$, and CP symmetry is preserved confirming the result of the Vafa-Witten theorem.

Another case where an exact solution is also known is the ℂℙN model in the large $`N`$ limit adl . In this case the system also describes a weakly interacting and parity preserving massive scalar particle in the adjoint representation of $`SU(N)`$. Let us analyse in some detail how the theorem works in this particular case. The effective action in a topological sector $`q0`$ with periodic boundary conditions, is given at leading order in $`N`$ by

$`S_{\mathrm{eff}}(m,q)={\displaystyle \frac{Nm^2V}{4\pi }}\mathrm{log}{\displaystyle \frac{4\pi |q|}{\mu ^2V}}qN\mathrm{log}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{1}{2}+\frac{m^2V}{4\pi |q|}\right)}{\sqrt{2\pi }}}{\displaystyle \frac{Nm^2V}{2g^2}}`$ (6)
where $`2m`$ is the mass of the physical particles, $`\mu =m\mathrm{exp}(2\pi /g^2)`$ is an energy scale fixing the value of the coupling constant $`g`$, $`V=LT`$ is the volume of the Euclidean space-time 𝕋2 and $`\mathrm{\Gamma }`$ is the Euler gamma function. The effective action has a simple linear asymptotic regime for large values of the topological charge $`q>m^2V`$,

$`S_{\mathrm{eff}}(m,q)={\displaystyle \frac{Nm^2V}{4\pi }}\mathrm{log}{\displaystyle \frac{4\pi |q|}{\mu ^2V}}+{\displaystyle \frac{N|q|}{2}}\mathrm{log}2+𝒪(1)`$ (7)
The radius of convergence of the partition function $`𝒵_\theta `$ in the pure imaginary axis $`\theta =\dot{ı}\vartheta `$ is $`\vartheta _c=N/2\mathrm{log}2`$. The critical value $`\theta _c`$ increases with $`N`$, and in the large $`N`$ limit the whole pure imaginary axis is free of Lee-Yang singularities in agreement with the general theorem. In the complex $`\theta `$-plane the radius of convergence around the origin is larger than $`\pi /4`$ even in the infinite volume limit, which also confirms the Vafa-Witten theorem in the large $`N`$ limit. Now, the theorem holds not only in that limit, but for any finite value of $`N`$ as well.

Finally, we remark that the asymptotic behavior of the effective action for small values of the topological charge $`q<m^2V`$ is quadratic instead of linear in $`q`$,

$`S_{\mathrm{eff}}(m,q)={\displaystyle \frac{Nm^2V}{4\pi }}(1\mathrm{log}{\displaystyle \frac{m^2}{\mu ^2}})+{\displaystyle \frac{N\pi q^2}{6m^2V}}+𝒪(q^4).`$ (8)
Since $`V=LT`$ the change of asymptotic behavior from (7) to (8) can be associated with a finite temperature $`\beta =1/T`$ crossover from the low temperature regime $`\beta >m^2L/|q|`$ to the high temperature regime $`\beta <m^2L/|q|`$ and cannot be related to any phase transition aff . One might expect a similar behavior in QCD, although in that case there is a finite temperature phase transition kpt .

The arguments used above cannot exclude the existence of higher-order phase transitions which do not require spontaneous parity symmetry breaking. Indeed, if $`\rho _0=0`$ but $`\rho (z)\rho _1z_2^\alpha \delta (z_1mz_2^\beta )`$ with $`0<\alpha 1`$ for $`|z|1`$, the vacuum energy density presents a second order singularity and the topological susceptibility diverges at $`\vartheta =0`$. This behavior is compatible with the Vafa-Witten inequality and might explain the behavior of the topological susceptibility of the ℂℙ1 model in numerical simulations lh .

We thank J. Barbón, I. Cirac, V. Laliena, E. Seiler and E. Witten for illuminating discussions.