text:
# New Physics Beyond the Standard Model at 𝛾𝛾 Colliders

## 1 Introduction

As is well known there are two ways in which new physics beyond the SM may manifest itself. In the direct scenario, a new particle (or set of particles) may be singly or pair produced at a collider. In the indirect scenario, precision measurements would determine at a high level of confidence that deviations from the SM are found for a particular observable or set of observables. This set of deviations may or may not be sufficient in their magnitude or direction to point to the correct new physics scenario from which they arose.

In the case of direct production it is important to note that the cross sections for pairs of scalars, fermions or vectors particles are all significantly larger at $`\gamma \gamma `$ colliders than they are at $`e^+e^{}`$ colliders (even after angular acceptance cuts) with comparable luminosities anticipated at both machines. However, by comparison, the $`\gamma \gamma `$ search reach is reduced due to the high energy cutoff in the photon fluxes. While polarized $`\gamma \gamma `$ colliders allow one to isolate the couplings of new physics to photons and probe spin configurations not accessible in $`e^+e^{}`$ collisions, they cannot produce neutral particles except via loops, though such cross sections may in some cases be large. For indirect searches both polarized $`e^+e^{}`$ and $`\gamma \gamma `$ colliders lead to a similar number of final states which have somewhat comparable statistical power as far as looking for deviations are concerned. However, in some cases the $`\gamma \gamma `$ collider has the edge due to the larger cross sections. Based on these arguments alone, and not knowing the source of new physics a priori, we would expect a general complementarity for $`e^+e^{}`$ and $`\gamma \gamma `$ for new physics searches. Unfortunately, analyses of new physics scenarios are not yet as evolved in $`\gamma \gamma `$ collisions as they are for $`e^+e^{}`$.

## 2 Some Examples

There are a number of examples of new physics scenarios which clearly display this complementarity. Perhaps the most well known amongst these possibilities are the anomalous gauge couplings of the $`W`$ and searches for leptoquarks. While there have been extensive analysesrev  of anomalous triple gauge boson couplings at $`e^+e^{}`$ colliders, the corresponding processes $`\gamma eW\nu `$ and $`\gamma \gamma WW`$ have generally gotten less attention in the literatureold . Note that while $`\gamma eW\nu `$ and $`\gamma \gamma WW`$ isolate the anomalous photon couplings to the $`W`$, $`e^+e^{}WW`$ also potentially involves anomalous $`Z`$ couplings which may also be present. This process is also dominated by the large $`t`$-channel neutrino exchange diagram which can be removed using beam polarization. Perhaps the best example of this complementarity is displayed in the study of Choi and Shrempp that compares the sensitivities of these three processes to the anomalous couplings $`\mathrm{\Delta }\kappa _\gamma =1\kappa _\gamma `$ and $`\lambda _\gamma `$ which is shown in Fig.1. Though now outdated, one sees from this analysis that the sensitivities of the three experiments are qualitatively similar and their overlapping region is substantially smaller than that obtained from any single measurement. It would be very interesting, and perhaps quite important, to repeat this analysis using modern luminosities which are more than an order of magnitude larger than used in this study and to include the increased set of observables that can be used to constrain these anomalous couplings. It is easy to imagine that the resulting area of overlap would now be smaller by more than two orders of magnitude. In the case of anomalous quartic gauge couplings, they can be probed at $`\gamma \gamma `$ colliders through processes which occur at order $`\alpha ^2`$ whereas they can only be accessed in $`e^+e^{}`$ collisions at order $`\alpha ^3`$ giving a great advantage to $`\gamma \gamma `$ colliders.

Leptoquarkslq  have been well studied in $`e^+e^{}`$, $`\gamma e`$ and $`\gamma \gamma `$ colliders. Buchmüller, Rückl and Wyler have long ago shown that under reasonably general assumptions leptoquarks can arise in any of 14 different varieties. Leptoquarks may have fermion number $`F=0,\pm 2`$, come in singlets, doublets or triplets of weak isospin and are either scalars or vectors. (There may even be some evidence for the existence of leptoquarks in the mass range near 400 GeVzarn .) Since the LHC will run before any high energy $`e^+e^{}`$ or $`\gamma \gamma `$ collider and discover(or not) any leptoquarks in the mass range accessible to these machines, their true role is not leptoquark discovery but leptoquark identification. Assuming these particles are elementary, a given leptoquark species is easily identified in $`e^+e^{}`$ collisions through the energy dependence of it’s total cross section and angular distribution. However, if form factors are present and/or leptoquarks have gauge strength Yukawa couplings it is possible that confusion can arise and a polarized $`\gamma \gamma `$ collider will be needed to disentangle the various ID possibilities. Since leptoquark charges, $`Q`$, vary in magnitude between $`1/3`$ and $`5/3`$ and their production cross sections at $`\gamma \gamma `$ are proportional to $`Q^4`$, a wide range of rates can be expected. This is shown explicitly in Fig.2 from bel  for a $`\sqrt{s}`$=1 TeV collider.

Note that similar complementarity between $`e^+e^{}`$ and $`\gamma \gamma `$ colliders also arises in the case of excited fermion productionstu .

## 3 $`\gamma \gamma WW`$ in Theories with Large Extra Dimensions

In theories with large extra dimensionsadd  the exchange of Kaluza-Klein(KK) graviton towers result in large set of new dimension-8 operators that can lead to new contributions to the matrix elements for a number of different processespheno . These exchanges can cause substantial deviations in cross sections, angular distributions and asymmetries. These operators in principle have an arbitrary sign, $`\lambda =\pm 1`$, and an associated cutoff mass scale, $`M_H`$, which is of order the higher dimensional Planck scale. Of the many processes examined so far, $`\gamma \gamma WW`$ provides the greatest reach for $`M_H`$ in comparison to the overall collider center of mass energy. (It is difficult to compare the $`\sqrt{s}`$ dependence of this reach for hadron colliders and $`e^+e^{}`$/$`\gamma \gamma `$ colliders.) The main reasons for this are that the $`WW`$ final state offers many observables which are particularly sensitive to the initial electron and laser polarizations as well as the very high statistics available (due to the 80 pb cross section) with which to probe graviton KK contributions. The differential cross sections are shown in Fig.3 for the SM as well as when the KK tower is included for $`\lambda =1`$.

Note that within the SM there is no dramatically strong sensitivity to the initial state lepton and laser polarizations and all of the curves have roughly the same shape. When the graviton tower contributions are included there are several effects. First, all of differential distributions become somewhat more shallow at large scattering angles but there is little change in the forward and backward directions due to the dominance of the SM poles. Second, there is now a clear and distinct sensitivity to the initial state polarization selections. In some cases, particularly for the $`(++)`$ and $`(++)`$ helicity choices, the differential cross section increases significantly for angles near $`90^o`$ taking on an m-like shape. This shape is, in fact, symptomatic of the spin-2 nature of the K-K graviton tower exchange since a spin-0 exchange leads only to a flattened distribution. Given the very large statistics available with a typical integrated luminosity of 100-300 $`fb^1`$, it is clear that the $`\gamma \gamma W^+W^{}`$ differential cross section is quite sensitive to $`M_H`$ especially for the two initial state helicities specified above.

In the SM, the final state $`W`$’s are dominantly transversely polarized. Due to the nature of the spin-2 graviton exchange, the K-K tower leads to a final state where both $`W`$’s are instead completely longitudinally polarized. This is shown explicitly in Fig.4 where we see the fraction of longitudinally polarized $`W`$’s falling rapidly as the scale $`M_H`$ is increased in comparison to $`\sqrt{s}`$. Thus we might expect that a measurement of the $`W`$ polarization will be a sensitive probe $`M_H`$. Other observables can be constructed out of the polarization dependent cross sections themselves. For the six possible initial state polarizations one can construct three independent polarization asymmetries of the form

$$A_{pol}(z)=\frac{d\sigma (P_{e1},P_{l1},P_{e2},P_{l2})d\sigma (P_{e1},P_{l1},P_{e2},P_{l2})}{d\sigma (P_{e1},P_{l1},P_{e2},P_{l2})+d\sigma (P_{e1},P_{l1},P_{e2},P_{l2})},$$
(1)
where $`z=\mathrm{cos}\theta `$, which are shown in Fig.4. These asymmetries are not only functions of $`z`$ but are also dependent on the $`W`$ pair invariant mass $`M_{WW}`$ in a way that is also sensitive to graviton exchange as shown in the left panel of Fig.5.

By performing a combined fit to the total cross sections and angular distributions, the $`LL`$ and $`LT+TL`$ helicity fractions for various initial state polarization choices and the polarization asymmetries we are able to discern the discovery reaches for $`M_H`$ as a function of the total $`\gamma \gamma `$ integrated luminosity; this is shown in Fig.6. Here we see a search reach in the range of $`M_H1113\sqrt{s}`$, which is the larger than that obtained from all other processes examined so far. By comparison, a combined analysis of the processes $`e^+e^{}f\overline{f}`$ with the same integrated luminosity leads to a search reach of $`67\sqrt{s}`$.

We note in passing that other $`\gamma \gamma `$ final states are also sensitive to graviton exchange, one example being $`ZZ`$ final state, which yield smaller search reaches. The right panel of Fig.5 shows the differential cross section due to graviton exchange for this process for different polarization states. These values are significantly larger than those arising from the SM.

## 4 Conclusion

As can be seen from the discussion above, $`e^+e^{}`$ and $`\gamma \gamma `$ colliders are quite complementary when it comes to discovering and exploring new physics scenarios. In some cases, such as graviton exchange, the $`\gamma \gamma `$ reach is superior to that of other colliders.


format: text

text:
# QCD Factorization For B Decays To Two Light Pseudoscalars Including Chirally Enhanced Corrections *footnote **footnote *Supported in part by National Natural Science Foundation of China and State Commission of Science and Technology of China

## Acknowledgements

This work is supported in part by National Natural Science Foundation of China and State Commission of Science and Technology of China.


text:
# All-Sky Convolution for Polarimetry Experiments

## I Introduction

Over the past decade a number of increasingly sophisticated experiments have reported detections of the temperature anisotropy in the cosmic microwave background (CMB). Following the detection of degree-scale anisotropy by the COBE satellite , ground-based and balloon-borne experiments have pushed back the limits on resolution and sensitivity to provide estimates of the anisotropy power spectrum $`C_l^\text{T}`$ up to Legendre multipoles of $`l700`$ (see e.g. Ref.  for a review of the situation pre-BOOMERANG  and MAXIMA-1 ). The anisotropy power spectrum encodes a wealth of cosmological information (e.g. Ref. ) in a highly compressed form, making it a very convenient data product from which to determine cosmological parameters (see Ref.  and references therein).

The combination of Thomson scattering and the non-zero temperature quadrupole, as the radiation begins to free stream through recombination, leads to the robust prediction that the CMB should be linearly polarized, with an r.m.s. level of a few percent of the temperature anisotropies . Detection of polarization would provide complementary information to that obtained from temperature measurements, e.g. the unique polarization signature of gravitational waves (and vector modes) provides the best hope of detecting their presence at last scattering . Currently, only upper limits exist on the degree of linear polarization (e.g. Refs. ; see also Ref.  and references therein for a recent review) but detections should be made with the MAP  and Planck  satellites, several experiments from the ground , and the flights of the MAXIPOL  and enhanced BOOMERANG  balloon experiments.

In an ideal linear polarimetry experiment, a given detector is configured to respond only to a single component of the electric field of the incident radiation, along incident directions contained within a small solid angle $`\mathrm{\Delta }\mathrm{\Omega }`$ (the beam size). For such an ideal experiment, the detector measures the flux $`(I+Q)\mathrm{\Delta }\mathrm{\Omega }/2`$ where $`I`$ (total intensity) and $`Q`$ are the Stokes parameters of the incident radiation along the beam direction, and the polarization basis vectors have been chosen with the $`x`$ direction aligned with the polarimeter. In practice, such ideals are never achieved, and for precision polarimetry experiments, e.g. Planck, it is essential to take full account of several beam effects on the measured signal. The polarization on the sky must be convolved with the response pattern of the detector, which will not be perfectly axisymmetric. Furthermore, the system will generally have some cross-polar contamination (non-zero response to more than one polarization component). Often the beam response pattern has non-negligible far side-lobes which are highly polarized due to reflection and diffraction effects in the instrument, making a full-sky convolution necessary if the effects of stray light from bright regions (such as the sun, moon and galaxy in CMB experiments) are to be properly accounted for. The implications of a subset of these effects for the analysis pipeline of total power (unpolarized) experiments have been studied recently in Ref. .

In this paper we present a formalism that allows all such non-ideal effects to be taken account of exactly and efficiently in multipole space. The form of the convolution for polarized data in multipole space is very similar to that for unpolarized radiation. This fact allows the recently suggested algorithm of Wandelt and Górski  to be used to compute rapidly the detector output for an arbitrary pointing direction and orientation of the detector. The formalism presented here should prove useful in simulation and modelling of precision polarimetry experiments, as well as in the actual analysis of experimental data. Within the multipole formalism, it is simple to invoke approximations, such as axisymmetry of the beam, where appropriate, to reproduce the approximate results used in previous analyses of beam effects in CMB polarimetry experiments .

The paper is arranged as follows. In Sec. II we discuss the derivation of the detector response from its far-field radiation pattern, and introduce the beam response tensor. The multipole expansion of the beam in spherical scalar and tensor harmonics is presented in Sec. III. The scanning of the detector on the sky is described by a (time-dependent) rotation of the beam from some standard configuration. This rotation is most conveniently handled in multipole space, and is described in Sec. IV, where an efficient algorithm for performing the rotation and convolution is also described. Section V deals with the case of an axisymmetric beam. It is shown that for a certain geometry of the polarized beam response, it is possible to describe the convolution of the beam with the incident linear polarization in terms of spin-weight 2 window functions, first introduced in Ref. . We give the window function for Gaussian beams of arbitrary angular size, and show that in the small-angle limit we reproduce the results of Ng and Liu . We summarise our discussion in Sec. VI. An appendix provides some details on the rotation properties of the tensor harmonics which are needed for Sec. IV.

## II Far-field radiation pattern and the beam response

It is convenient to characterise the response of the detector and feed system by the far-field radiation pattern it would emit if used as a transmitter, rather than a receiver . Assuming a (quasi-)monochromatic system, the electric field in the far field is of the form

$$𝑬\frac{1}{r}\mathrm{}\{\stackrel{~}{𝓔}\mathrm{exp}[i(kr\omega t)]\},$$
(1)
where $`\stackrel{~}{𝓔}`$ is a complex, transverse vector function on the sphere. (We refer to geometric objects having no components outside the surface of the sphere as being transverse.) Here $`r`$ is radial distance and $`\omega =ck`$ is the (mean) angular frequency of the radiation, where $`k1/r`$ is the wavenumber. The most general detector and feed system will produce a partially polarized signal in transmission, so that $`\stackrel{~}{𝓔}`$ will generally be a slowly varying function of time (compared to $`\omega `$). However, only the stationary statistical properties of $`\stackrel{~}{𝓔}`$ are important for determining the response of the system to incident radiation. With the system in some specified orientation, the power received $`\text{d}W_{\text{tot}}`$ when illuminated by the sky along some direction $`𝒆`$ is proportional to the intensity in that polarization mode which is the time reverse of Eq. (1). It follows that

$$\text{d}W_{\text{tot}}|𝓔\stackrel{~}{𝓔}|^2\text{d}\mathrm{\Omega },$$
(2)
where $`𝓔`$ is the complex representative (or analytic signal) of the incident electric field propagating along $`𝒆`$, $`\text{d}\mathrm{\Omega }`$ is the element of solid angle, and angle brackets denote time averaging over the slow variations in $`𝓔`$ and $`\stackrel{~}{𝓔}`$. Writing the components of the fields $`𝓔`$ and $`\stackrel{~}{𝓔}`$ on the orthonormal basis vectors $`𝝈_\theta `$ and $`𝝈_\varphi `$ of a spherical polar coordinate system as e.g. $`_\theta `$ and $`_\varphi `$, the contribution to $`W_{\text{tot}}`$ can be written as

$$\frac{\text{d}W_{\text{tot}}}{\text{d}\mathrm{\Omega }}\frac{1}{2}(I\stackrel{~}{I}+Q\stackrel{~}{Q}+U\stackrel{~}{U}V\stackrel{~}{V}),$$
(3)
where $`\{I,Q,U,V\}`$ are the Stokes parameters of the incoming radiation on the $`\{𝝈_\theta ,𝝈_\varphi \}`$ basis:

$`I`$ $`=`$ $`|_\theta |^2+|_\varphi |^2,`$ (4)
$`Q`$ $`=`$ $`|_\theta |^2|_\varphi |^2,`$ (5)
$`U`$ $`=`$ $`2\mathrm{}_\theta _\varphi ^{},`$ (6)
$`V`$ $`=`$ $`2\mathrm{}i_\theta _\varphi ^{},`$ (7)
and $`\{\stackrel{~}{I},\stackrel{~}{Q},\stackrel{~}{U},\stackrel{~}{V}\}`$ are effective Stokes parameters for the beam:

$`\stackrel{~}{I}`$ $`=`$ $`|\stackrel{~}{}_\theta |^2+|\stackrel{~}{}_\varphi |^2,`$ (8)
$`\stackrel{~}{Q}`$ $`=`$ $`|\stackrel{~}{}_\theta |^2|\stackrel{~}{}_\varphi |^2,`$ (9)
$`\stackrel{~}{U}`$ $`=`$ $`2\mathrm{}\stackrel{~}{}_\theta \stackrel{~}{}_\varphi ^{},`$ (10)
$`\stackrel{~}{V}`$ $`=`$ $`2\mathrm{}i\stackrel{~}{}_\theta \stackrel{~}{}_\varphi ^{}.`$ (11)
Note that the effective Stokes parameters for the beam are defined on the same basis as the incoming radiation, which is responsible for the minus sign in front of the $`V\stackrel{~}{V}`$ term in Eq. 3.

The intensity $`I`$ and circular polarization $`V`$ are invariant under rotations of the polarization basis vectors, whilst $`Q`$ and $`U`$ transform like the components of a second-rank tensor. Introducing the linear polarization tensor for the incident radiation

$`𝒫^{ab}(𝒆)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[Q(𝝈_\theta 𝝈_\theta 𝝈_\varphi 𝝈_\varphi )`$ (13)
$`U(𝝈_\theta 𝝈_\varphi +𝝈_\varphi 𝝈_\theta )],`$
we can write the total power received in the basis-independent form

$$W_{\text{tot}}\frac{1}{2}(I\stackrel{~}{I}V\stackrel{~}{V}+2𝒫_{ab}^{ab})\text{d}\mathrm{\Omega }.$$
(14)
Here, $`_{ab}`$ is the (linear) beam response tensor:

$`^{ab}(𝒆)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[\stackrel{~}{Q}(𝝈_\theta 𝝈_\theta 𝝈_\varphi 𝝈_\varphi )`$ (16)
$`\stackrel{~}{U}(𝝈_\theta 𝝈_\varphi +𝝈_\varphi 𝝈_\theta )]`$
$`=`$ $`\mathrm{}[\stackrel{~}{𝓔}\stackrel{~}{𝓔}^{}]^{\text{TT}},`$ (17)
where $`^{\text{TT}}`$ denotes the transverse, trace-free part. Our main task now is to derive the dependence of the total power received on the pointing direction and orientation of the detector.

### A Co- and cross-polarized basis

The polar basis $`\{𝝈_\theta ,𝝈_\varphi \}`$ is fixed relative to the sky and is singular at the north and south poles. For describing the beam, it is standard practice to use an alternative basis which is fixed relative to the detector, and has its only singularity in the opposite direction to the main beam . We define a set of Cartesian basis vectors $`\{𝝈_x^{},𝝈_y^{},𝝈_z^{}\}`$ which are fixed relative to detector. It is convenient to take $`𝝈_z^{}`$ to be along the (nominal) main beam, and $`𝝈_y^{}`$ along the polarization direction on axis. Using this Cartesian frame we derive a set of polar basis vectors $`\{𝝈_\theta ^{},𝝈_\varphi ^{}\}`$ on the sphere in the standard manner. The co- and cross-polar basis vectors are then derived by parallel-transporting $`𝝈_y^{}`$ and $`𝝈_x^{}`$ respectively from the north pole along great circles through the poles:

$`𝝈_{\text{co}}`$ $`=`$ $`\mathrm{sin}\varphi ^{}𝝈_\theta ^{}+\mathrm{cos}\varphi ^{}𝝈_\varphi ^{}`$ (18)
$`𝝈_{\text{cross}}`$ $`=`$ $`\mathrm{cos}\varphi ^{}𝝈_\theta ^{}\mathrm{sin}\varphi ^{}𝝈_\varphi ^{},`$ (19)
where $`\theta ^{}`$ and $`\varphi ^{}`$ are spherical polar coordinates. A well-defined linearly polarized receiver has $`|\stackrel{~}{}_{\text{cross}}||\stackrel{~}{}_{\text{co}}|`$ along the main beam, where e.g. $`\stackrel{~}{}_{\text{co}}`$ is the component of $`\stackrel{~}{𝓔}`$ along $`𝝈_{\text{co}}`$. Cross-polar contamination arises from a non-zero $`|\stackrel{~}{}_{\text{cross}}|`$. To rotate from the co- and cross-polarized basis to the spherical polar basis we have to rotate through $`\pi /2\varphi ^{}`$ in a right-handed sense about the inward normal to the sphere. Transforming the Stokes parameters for the beam from the co- and cross-polar basis to the spherical polar basis, we have:

$`\stackrel{~}{I}`$ $`=`$ $`|\stackrel{~}{}_{\text{co}}|^2+|\stackrel{~}{}_{\text{cross}}|^2,`$ (20)
$`\stackrel{~}{Q}`$ $`=`$ $`|\stackrel{~}{}_{\text{co}}|^2|\stackrel{~}{}_{\text{cross}}|^2\mathrm{cos}2\varphi ^{}`$ (22)
$`+2\mathrm{}\stackrel{~}{}_{\text{co}}\stackrel{~}{}_{\text{cross}}^{}\mathrm{sin}2\varphi ^{},`$
$`\stackrel{~}{U}`$ $`=`$ $`|\stackrel{~}{}_{\text{co}}|^2|\stackrel{~}{}_{\text{cross}}|^2\mathrm{sin}2\varphi ^{}`$ (24)
$`2\mathrm{}\stackrel{~}{}_{\text{co}}\stackrel{~}{}_{\text{cross}}^{}\mathrm{cos}2\varphi ^{},`$
$`\stackrel{~}{V}`$ $`=`$ $`2\mathrm{}i\stackrel{~}{}_{\text{co}}\stackrel{~}{}_{\text{cross}}^{},`$ (25)
which reflect the spin-2 nature of the linear polarization. For simulation purposes, $`\stackrel{~}{𝓔}`$ is usually determined with physical optics codes, the results being reported on the co- and cross-polar basis. For real experiments the Stokes parameters for the beam must be calibrated using sources with known surface brightness and polarization. Note that we do not assume that the beam response is fully polarized, so the formalism developed here can also be applied to total power experiments (only $`\stackrel{~}{I}`$ non-zero). In practice, the optics and feeds will introduce some beam polarization in the side-lobes even for a nominal total power experiment. Although small, the role of such effects in total power experiments could be quantified with our formalism.

## III Multipole expansions

The dependence of the total power received on the direction and orientation of the telescope is most easily formulated in multipole space. We say that the detector is in its reference orientation when it is oriented so that the basis $`\{𝝈_x^{},𝝈_y^{},𝝈_z^{}\}`$ coincides with the $`\{𝝈_x,𝝈_y,𝝈_z\}`$ basis, which is fixed relative to the sky. We describe the beam via a set of constant multipole coefficients which are extracted on the sky basis when the detector is in its reference orientation. To describe an arbitrary orientation of the detector at some time along the scan, we can rotate the beam (which is most easily performed in multipole space) to obtain the rotated beam response which is convolved with the sky, as in Eq. (14).

In the reference orientation, the total intensity and circular polarization parts of the beam response can be expanded in scalar spherical harmonics, e.g.

$$\stackrel{~}{I}(𝒆)=\underset{lm}{}b_{(lm)}^IY_{(lm)}(𝒆),$$
(26)
where the sum is over $`l0`$ and $`|m|l`$. The multipoles $`b_{(lm)}^V`$ for the circular polarization are defined analogously. For the beam response tensor $`_{ab}`$ we must expand in the transverse, trace-free tensor harmonics. Here we follow the coordinate-dependent approach of Ref.  (although for some applications the coordinate-free approach of Ref.  is more convenient):

$$_{ab}(𝒆)=\underset{\text{P}lm}{}b_{(lm)}^\mathrm{P}Y_{(lm)ab}^\mathrm{P},$$
(27)
where the sum is over $`l2`$, $`|m|l`$, and the two types of transverse trace-free harmonics $`\text{P}=\text{G}`$ (for Gradient, often called electric) or C (Curl, often called magnetic). All multipoles satisfy $`b_{(lm)}^{}=(1)^mb_{(lm)}`$ since the fields are real and we have adopted the Condon-Shortley phase for the spherical harmonics. The Stokes parameters $`I`$ and $`V`$ for the sky can be similarly expanded in multipoles $`a_{(lm)}^I`$ and $`a_{(lm)}^V`$, and the linear polarization in multipoles $`a_{(lm)}^\mathrm{G}`$ and $`a_{(lm)}^\mathrm{C}`$.

The tensor harmonics are derived from the scalar harmonics by covariant differentiation over the sphere  (see also the Appendix). Performing the differentiation gives

$`Y_{(lm)}^{\mathrm{G}ab}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left({}_{2}{}^{}Y_{(lm)}^{}𝒎𝒎+{}_{2}{}^{}Y_{(lm)}^{}𝒎^{}𝒎^{}\right),`$ (28)
$`Y_{(lm)}^{\mathrm{C}ab}`$ $`=`$ $`{\displaystyle \frac{1}{i\sqrt{2}}}\left({}_{2}{}^{}Y_{(lm)}^{}𝒎𝒎{}_{2}{}^{}Y_{(lm)}^{}𝒎^{}𝒎^{}\right),`$ (29)
where $`𝒎(𝝈_\theta +i𝝈_\varphi )/\sqrt{2}`$. The $`{}_{\pm 2}{}^{}Y_{(lm)}^{}`$ are the spin-weight 2 harmonics:

$${}_{\pm 2}{}^{}Y_{(lm)}^{}=\frac{N_l}{\sqrt{2}}(W_{(lm)}\pm iX_{(lm)}),$$
(30)
where $`N_l[2(l2)!/(l+2)!]^{1/2}`$, and

$`W_{(lm)}(\theta ,\varphi )`$ $`=`$ $`{\displaystyle \frac{^2}{\theta ^2}}Y_{(lm)}\mathrm{cot}\theta {\displaystyle \frac{}{\theta }}Y_{(lm)}`$ (32)
$`\mathrm{csc}^2\theta {\displaystyle \frac{^2}{\varphi ^2}}Y_{(lm)},`$
$`X_{(lm)}(\theta ,\varphi )`$ $`=`$ $`2\mathrm{csc}\theta \left({\displaystyle \frac{^2}{\theta \varphi }}Y_{(lm)}\mathrm{cot}\theta {\displaystyle \frac{}{\varphi }}Y_{(lm)}\right).`$ (33)
Explicit expressions for $`W_{(lm)}`$ and $`X_{(lm)}`$ with the derivatives eliminated are given in Ref. . Note that our convention for the spin-weight functions follows Ref. , which differs from Goldberg et al.  by the inclusion of the factor $`(1)^m`$ in the definition of the spherical harmonics.

### A Extracting the beam multipoles

With the detector in the reference orientation, the spherical polar bases fixed relative to the sky and detector coincide. We can extract the beam multipoles from the effective Stokes parameters on this polar basis using the orthonormality of the scalar and tensor harmonics (e.g. Ref. ). For the linear polarization, we have

$$b_{(lm)}^\mathrm{G}\pm ib_{(lm)}^\mathrm{C}=\frac{1}{\sqrt{2}}(\stackrel{~}{Q}i\stackrel{~}{U}){}_{\pm 2}{}^{}Y_{(lm)}^{}\text{d}\mathrm{\Omega },$$
(34)
which is the inverse of the expansion of $`\stackrel{~}{Q}`$ and $`\stackrel{~}{U}`$ in spin-weight 2 harmonics:

$$\frac{1}{\sqrt{2}}(\stackrel{~}{Q}\pm i\stackrel{~}{U})=\underset{lm}{}(b_{(lm)}^\mathrm{G}ib_{(lm)}^\mathrm{C}){}_{2}{}^{}Y_{(lm)}^{}.$$
(35)
(The sum is over $`l2`$ and $`|m|l`$.) Note that with our conventions for the polarization basis vectors, $`\stackrel{~}{Q}\pm i\stackrel{~}{U}`$ is a spin-weight $`2`$ quantity, which differs from some authors (notably Ref. ).

Some care is needed in extracting the beam multipoles at the north and south poles since the Stokes parameters on the polar basis are ill-defined there. However, using Eqs. (22) and (24), it is simple to show that any well-defined polarization field must have $`\varphi `$ dependence going like $`Q\pm iU\mathrm{exp}(\pm 2i\varphi )`$ on the polar basis at the north pole ($`\theta =0`$). This is consistent with Eq. (35) since the spin-weight $`\pm 2`$ harmonics satisfy

$${}_{\pm 2}{}^{}Y_{(lm)}^{}=\delta _{m2}\sqrt{\frac{(2l+1)}{4\pi }}\text{e}^{2i\varphi }$$
(36)
at $`\theta =0`$. With finitely sampled, simulated data for $`\stackrel{~}{}_{\text{co}}`$ and $`\stackrel{~}{}_{\text{cross}}`$, the contribution to $`b_{(lm)}^\mathrm{P}`$ from samples on (or very near to) the north pole can be treated by approximating $`{}_{\pm 2}{}^{}Y_{(lm)}^{}`$ with Eq. (36), and absorbing the $`\mathrm{exp}(\pm 2i\varphi )`$ into $`(\stackrel{~}{Q}i\stackrel{~}{U})\mathrm{exp}(\pm 2i\varphi )`$ which is then well-defined by the data. Similar problems occur at the south pole, but since the beam has virtually no power there this is not problematic.

## IV Beam rotation and convolution

The kinematics of the experiment can be specified by a scan strategy which describes the rotation necessary to take each detector from the reference orientation to its orientation at the specified time in the scan. For simplicity we consider only a single detector, but our approach could easily be generalised to experiments with multiple detectors. The rotation is specified by its Euler angles $`\{\varphi ,\theta ,\psi \}`$, such that first we rotate in a right-handed sense by $`\psi `$ about $`𝝈_z`$, then by $`\theta `$ about $`𝝈_y`$, and finally by $`\varphi `$ about $`𝝈_z`$ again. We denote the rotation by $`D(\varphi ,\theta ,\psi )`$ so that the image of $`𝝈_i`$ is $`D(\varphi ,\theta ,\psi )𝝈_i=𝝈_i^{}`$ for $`i=x`$, $`y`$, $`z`$. For $`\psi =0`$,

$`𝝈_z^{}`$ $`=`$ $`\mathrm{sin}\theta (\mathrm{cos}\varphi 𝝈_x+\mathrm{sin}\varphi 𝝈_y)+\mathrm{cos}\theta 𝝈_z`$ (37)
$`𝝈_y^{}`$ $`=`$ $`\mathrm{sin}\varphi 𝝈_x+\mathrm{cos}\varphi 𝝈_y.`$ (38)
For non-zero $`\psi `$, $`𝝈_x^{}`$ and $`𝝈_y^{}`$ are additionally rotated by $`\psi `$ about $`𝝈_z^{}`$ (see Fig. 1).

Under the rotation $`D(\varphi ,\theta ,\psi )`$ the beam Stokes parameters $`\stackrel{~}{I}`$ and $`\stackrel{~}{V}`$ rotate as scalar fields so that e.g.

$$\stackrel{~}{I}(𝒆)\stackrel{~}{I}[D(\psi ,\theta ,\varphi )𝒆],$$
(39)
where $`D(\psi ,\theta ,\varphi )`$ is the inverse rotation. The beam response tensor $`_{ab}(𝒆)`$ rotates as a rank-two tensor, so that

$$_{ab}(𝒆)\mathrm{\Lambda }_a{}_{}{}^{c_1}\mathrm{\Lambda }_{b}^{}{}_{}{}^{c_2}_{c_1c_2}^{}[D(\psi ,\theta ,\varphi )𝒆],$$
(40)
where $`\mathrm{\Lambda }_a^c`$ is the $`\text{SO}(3)`$ rotation matrix representing $`D(\varphi ,\theta ,\psi )`$ (see Appendix). Since we are describing the beams in multipole space, we must consider the transformation properties of the scalar and tensor harmonics under the rotations given in Eqs. (39) and (40) respectively. The scalar spherical harmonics transform irreducibly under rotations as

$$Y_{(lm)}(𝒆)\underset{|m^{}|l}{}D_{m^{}m}^l(\varphi ,\theta ,\psi )Y_{(lm^{})}(𝒆),$$
(41)
and, as we show in the Appendix, the same is true of the tensor harmonics:

$$Y_{(lm)ab}(𝒆)\underset{|m^{}|l}{}D_{m^{}m}^l(\varphi ,\theta ,\psi )Y_{(lm^{})ab}(𝒆).$$
(42)
Here, the $`D_{m^{}m}^l(\varphi ,\theta ,\psi )`$ are Wigner’s $`D`$-matrices. With our conventions for the Euler angles, we have

$$D_{m^{}m}^l(\varphi ,\theta ,\psi )=\text{e}^{im^{}\varphi }d_{m^{}m}^l(\theta )\text{e}^{im\psi },$$
(43)
where

$`d_{mn}^l(\theta )`$ $`=`$ $`{\displaystyle \underset{t}{}}(1)^t{\displaystyle \frac{[(l+m)!(lm)!(l+n)!(ln)!]^{1/2}}{(l+mt)!(lnt)!(t+nm)!t!}}`$ (45)
$`\times [\mathrm{cos}(\theta /2)]^{2l+mn2t}[\mathrm{sin}(\theta /2)]^{nm+2t}.`$
The sum is over integers $`t`$ such that the arguments of the factorials are non-negative.

Performing the integral over the sphere in Eq. (14) is now straightforward in multipole space using the orthonormality of the harmonics. The final result for the total power as a function of orientation of the detector is

$`W_{\text{tot}}(\varphi ,\theta ,\psi )`$ $``$ $`{\displaystyle \underset{lmm^{}}{}}[{\displaystyle \frac{1}{2}}(a_{(lm)}^Ib_{(lm^{})}^Ia_{(lm)}^Vb_{(lm^{})}^V)`$ (47)
$`+{\displaystyle \underset{\text{P}}{}}a_{(lm)}^\text{P}b_{(lm^{})}^\mathrm{P}]D^l_{mm^{}}(\varphi ,\theta ,\psi ).`$
The sum is over $`l0`$ with $`|m|`$ and $`|m^{}|l`$, where we defined the linear polarization multipoles to be zero for $`l=0`$ and $`1`$. Our result for the dependence of the total power on orientation is quite general; we have made no assumptions about the beam profile and level of cross-polar contamination. Equation (47) is one of the main results of this paper. Note that the function $`W_{\text{tot}}(\varphi ,\theta ,\psi )`$ is expressed as a linear combination of the $`D`$-matrices, which form a complete set for expanding single-valued (square-integrable) functions on the three-sphere.

### A Fast convolution algorithms

The right-hand side of Eq. (47) can be evaluated rapidly by making only minor modifications to the algorithm developed recently by Wandelt and Górski  for the case of an unpolarized detector. The key to the algorithm is to factor the rotation $`D(\varphi ,\theta ,\psi )`$ as follows:

$$D(\varphi ,\theta ,\psi )=D(\varphi \pi /2,\pi /2,\theta )D(0,\pi /2,\psi +\pi /2),$$
(48)
so we may write

$`D_{mm^{}}^l(\varphi ,\theta ,\psi )`$ $`=`$ $`{\displaystyle \underset{|M|l}{}}[D_{mM}^l(\varphi \pi /2,\pi /2,\theta )`$ (50)
$`\times D_{Mm^{}}^l(0,\pi /2,\psi +\pi /2)].`$
The advantage of factoring the rotation in this way is that now the Euler angles only occur in complex exponentials, and we only need evaluate $`d_{mm^{}}^l(\theta )`$ at $`\theta =\pi /2`$ \[since $`d_{mm^{}}^l(\theta )=d_{m^{}m}^l(\theta )`$\]. The full three-sphere of rotations can now be calculated with a three-dimensional Fast Fourier Transform. The $`d_{mm^{}}^l(\theta )`$ can be computed with the accurate recursive method in Ref.  (which can be further enhanced by making use of the symmetries of the $`d_{mm^{}}^l`$). To perform the convolution to a resolution corresponding to multipoles $`l_{\text{max}}`$, for all possible orientations of the system, requires $`O(l_{\text{max}}^4)`$ operations. This should be compared with the $`O(l_{\text{max}}^5)`$ operations required for a brute force computation in pixel space. For experiments such as the Planck mission, where $`l_{\text{max}}`$ is of the order of a few thousand, the saving is considerable. For many experiments, the approximate azimuthal symmetry of the beam limits the sum over $`m^{}`$ in Eq. (47) to $`|m^{}|l_{\text{max}}`$. Since $`W_{\text{tot}}`$ is then a slowly varying function of $`\psi `$, it is possible to sample the $`\psi `$ variation much more sparsely than for $`\theta `$ and $`\varphi `$, which effectively reduces the operations count to $`O(l_{\text{max}}^3)`$ .

## V Axisymmetric beams

In this section we consider the limiting form of the general result, Eq. (47), when the beam is approximated as being pure co-polarized (i.e. having no cross-polar contamination), and axisymmetric. Typically, these approximations hold well across the main beam for a well-defined linear polarimetry system. However, the approximations are not valid for describing the response of the system in the far side-lobes, where the complete expression, Eq. (47), should be used. To quantify the errors introduced in a given experiment by assuming an axisymmetric beam requires detailed simulation with the apparatus developed in the earlier sections of this paper. Techniques for propagating these errors to the band power estimates of the power spectrum in a total power experiment have been developed recently , but the extension to polarized experiments must await the detailed development of the full analysis pipeline for polarized data.

For a pure co-polar beam, $`\stackrel{~}{}_{\text{cross}}`$ vanishes in Eqs. (20)–(25). If we further assume axisymmetry, then $`\stackrel{~}{}_{\text{co}}`$ is a function of $`\theta ^{}`$ alone. Writing $`|\stackrel{~}{}_{\text{co}}|^2=B(\theta ^{})`$, we have $`\stackrel{~}{I}=B(\theta ^{})`$, $`\stackrel{~}{V}=0`$, and

$$\stackrel{~}{Q}\pm i\stackrel{~}{U}=B(\theta ^{})\text{e}^{\pm 2i\varphi ^{}}.$$
(51)
It follows that $`b_{(lm)}^\mathrm{P}=0`$ unless $`|m|=2`$. Furthermore, using Eq. (34) with $`m=2`$, we see that $`b_{(l2)}^\mathrm{C}=ib_{(l2)}^\mathrm{G}`$ and hence $`b_{(l2)}^\mathrm{C}=ib_{(l2)}^\mathrm{G}`$. The beam response tensor for linear polarization, $`_{ab}`$, is now fully specified by $`b_{(l2)}^\mathrm{G}`$. If we use Eqs. (32) and (33) to express the spin-weight harmonics in terms of Legendre functions, we find the following expression for $`b_{(l2)}^\mathrm{G}`$:

$`b_{(l2)}^\mathrm{G}`$ $`=`$ $`{\displaystyle \frac{\pi N_l^2}{\sqrt{2}}}\sqrt{{\displaystyle \frac{2l+1}{4\pi }}}{\displaystyle _1^1}\text{d}xB(x)\{(l+2)(x2)P_{l1}^{\prime \prime }`$ (53)
$`+[2(l1)x{\displaystyle \frac{1}{2}}l(l1)(1x^2)(l4)]P_l^{\prime \prime }\},`$
where primes denote differentiation with respect to $`x\mathrm{cos}\theta ^{}`$, and $`P_l(x)`$ is the $`l`$th Legendre polynomial. Note that $`b_{(l2)}^\mathrm{G}`$ is real, so $`b_{(l2)}^\mathrm{C}`$ is imaginary. For the intensity multipoles we have $`b_{(lm)}^I=0`$ unless $`m=0`$, with

$$b_{(l0)}^I=2\pi \sqrt{\frac{2l+1}{4\pi }}_1^1P_l(x)B(x)\text{d}x.$$
(55)
Given the restrictions on the beam multipoles for an axisymmetric, co-polar beam, it is possible to simplify Eq. (47) for the total power received. Using the relation

$$D_{ms}^l(\varphi ,\theta ,\psi )=(1)^s\sqrt{\frac{4\pi }{2l+1}}{}_{s}{}^{}Y_{(lm)}^{}(\theta ,\varphi )\text{e}^{is\psi }$$
(56)
between the $`D`$-matrices and the spin-weight $`s`$ harmonics for integer $`l`$, $`m`$, and $`s`$, and the reality of $`b_{(l0)}^I`$ and $`b_{(l2)}^\mathrm{G}`$, Eq. (47) can be written in the form

$$W_{\text{tot}}(\theta ,\varphi ,\psi )\frac{1}{2}[I_{\text{eff}}Q_{\text{eff}}\mathrm{cos}2\psi +U_{\text{eff}}\mathrm{sin}2\psi ].$$
(57)
Here, $`I_{\text{eff}}`$ is the usual beam smoothed intensity:

$$I_{\text{eff}}(\theta ,\varphi )=\underset{lm}{}W_la_{(lm)}^IY_{(lm)}(\theta ,\varphi ),$$
(58)
where the window function

$$W_l=\sqrt{\frac{4\pi }{2l+1}}b_{(l0)}^I.$$
(59)
Similarly, $`Q_{\text{eff}}`$ and $`U_{\text{eff}}`$ are beam smoothed Stokes parameters on the spherical-polar basis, given by

$$\frac{1}{\sqrt{2}}(Q_{\text{eff}}\pm iU_{\text{eff}})=\underset{lm}{}{}_{2}{}^{}W_{l}^{}(a_{(lm)}^\mathrm{G}ia_{(lm)}^\mathrm{C}){}_{2}{}^{}Y_{(lm)}^{},$$
(60)
where, following Ref. , we have introduced the spin-weight 2 window function

$${}_{2}{}^{}W_{l}^{}=2\sqrt{2}\sqrt{\frac{4\pi }{2l+1}}b_{(l2)}^\mathrm{G}.$$
(61)
Note that Eq. (57) shows that the signal obtained by convolving the pure co-polar, axisymmetric beam with the sky is equivalent to the response of an idealised co-polar detector, with vanishing beam width, on the smoothed sky. This result does not depend on any assumptions about the angular size of the beam response; for the polarized contribution it is a consequence of the definition of the co-polar vector field, Eq. (18), which is the obvious generalisation of a constant vector field to the surface of sphere. Equations (57)–(61) provide a complete description of the power received in polarimetry experiments with axisymmetric, co-polar beams.

### A Gaussian beams

It is often the case that the axisymmetric beam profile is approximately Gaussian:

$$B(\theta ^{})=B\mathrm{exp}[(1\mathrm{cos}\theta ^{})/\sigma ^2],$$
(62)
where $`\sigma `$ is a measure of the beam width. Note that we follow Bond and Efstathiou  in taking the beam to be Gaussian in $`2\mathrm{sin}(\theta ^{}/2)`$ rather than $`\theta ^{}`$. The former allows us to derive simple analytic results valid for any $`\sigma `$. However, in most cases a Gaussian profile is only appropriate close to the beam axis, in which case the two definitions are almost indistinguishable. For the beam profile in Eq. (62), the window functions $`W_l`$ and $`{}_{2}{}^{}W_{l}^{}`$ can be evaluated analytically:

$$W_l=4\pi B\text{e}^\alpha i_l(\alpha ),$$
(63)
and

$`{}_{2}{}^{}W_{l}^{}`$ $`=`$ $`2\pi BN_l^2\text{e}^\alpha \{2(1)^l[(l+2)(l1)+6\alpha ]\text{e}^\alpha `$ (66)
$`+[(l^24\alpha )(l1)^2+12\alpha ^2]i_l(\alpha )`$
$`+4\alpha (l^2+l+13\alpha )i_{l1}(\alpha )\},`$
where $`i_l(x)`$ is a modified spherical Bessel function, and $`\alpha 1/\sigma ^2`$. For $`\sigma ^21`$, which is always the case for high resolution experiments, the window functions are very well approximated by their asymptotic expansions, which give

$`W_l`$ $``$ $`2\pi B\sigma ^2\mathrm{exp}[l(l+1)\sigma ^2/2],`$ (67)
$`{}_{2}{}^{}W_{l}^{}`$ $``$ $`2\pi B\sigma ^2\mathrm{exp}\{[l(l+1)4]\sigma ^2/2\},`$ (68)
in full agreement with Ng and Liu . For $`l1`$, the polarized and unpolarized window functions in the small-scale limit are approximately equal, which can easily be verified by making a flat sky approximation.

## VI Conclusion

In this paper we have presented a multipole method for describing the response of an arbitrary detector and feed system. We fully include the effects of non-zero beam size, asymmetric beam patterns, and cross-polar contamination. Such inclusions are essential for the accurate modelling and interpretation of precision polarimetry data, such as that expected from the Planck mission. Working in multipole space, we derived a simple expression for the response of the system, when convolved with the sky, as a function of the three Euler angles needed to describe a general orientation of the system. Given the form of this expression, it is straightforward to modify the fast algorithm of Wandelt and Górski  to compute the system response for the entire three-sphere of orientations in $`O(l_{\text{max}}^4)`$ operations. Finally, we showed how, for the case of a pure co-polar, axisymmetric beam, the response can be described by spin-weighted window functions. This extended the results of Ref.  to arbitrary size beams, and we gave the exact form of the window functions for a beam with a Gaussian profile. Although our discussion has been in the context of CMB polarimetry experiments, the formalism introduced here should be useful in any applications involving anisotropic filtering of tensor fields on the sphere.

The techniques described in this paper have now been implemented in the simulation pipeline for polarized channels of the Planck mission.

## Acknowledgments

AC acknowledges a Research Fellowship from Queens’ College, Cambridge. PF is supported by a Research Fellowship from ESA. DM and MA are supported by PPARC. We thank the members of the Cambridge Planck Analysis Centre (CPAC) for many useful discussions.

## A Rotating the tensor harmonics

In this appendix we establish the transformation properties of the transverse, trace-free tensor harmonics under active rotations. The action of the rotation $`D(\varphi ,\theta ,\psi )`$ can be represented by a rotation matrix $`\mathrm{\Lambda }^a_b`$ such that an arbitrary vector $`v^a`$ rotates to $`v^a=\mathrm{\Lambda }^a{}_{b}{}^{}v_{}^{b}`$. Orthogonality of the rotation implies $`\mathrm{\Lambda }^a{}_{b}{}^{}\mathrm{\Lambda }_{a}^{}{}_{}{}^{c}=\delta _b^c`$, so the inverse rotation is $`v^a=v^b\mathrm{\Lambda }_b^a`$. Under an active rotation of the beam, the beam response tensor rotates to $`_{ab}^{}(e^c)`$ which is obtained by forward rotating the original tensor evaluated at the back rotated position:

$$_{ab}^{}(e^c)=\mathrm{\Lambda }_a{}_{}{}^{c_1}\mathrm{\Lambda }_{b}^{}{}_{}{}^{c_2}_{c_1c_2}^{}(e^d\mathrm{\Lambda }_d{}_{}{}^{c}).$$
(A1)
Under this transformation, the transverse nature of the tensor field $`_{ab}(e^c)`$ is preserved. We shall demand the same transformation properties for the tensor harmonics $`Y_{(lm)ab}^\mathrm{P}`$.

It is convenient to write the tensor harmonics in terms of covariant derivatives on the sphere of the scalar harmonics (e.g. Ref. ):

$`Y_{(lm)ab}^\mathrm{G}`$ $`=`$ $`N_l\left(\stackrel{~}{}_a\stackrel{~}{}_bY_{(lm)}{\displaystyle \frac{1}{2}}\stackrel{~}{g}_{ab}\stackrel{~}{}^2Y_{(lm)}\right)`$ (A2)
$`Y_{(lm)ab}^\mathrm{C}`$ $`=`$ $`{\displaystyle \frac{N_l}{2}}(\stackrel{~}{ϵ}^c{}_{b}{}^{}\stackrel{~}{}_{a}^{}\stackrel{~}{}_cY_{(lm)}+\stackrel{~}{ϵ}^c{}_{a}{}^{}\stackrel{~}{}_{b}^{}\stackrel{~}{}_cY_{(lm)}),`$ (A3)
where $`\stackrel{~}{}_a`$ is the covariant derivative on the unit sphere, $`\stackrel{~}{g}_{ab}=g_{ab}e_ae_b`$ is the (induced) metric (with $`g_{ab}`$ the metric of Euclidean 3-space) and $`\stackrel{~}{ϵ}_{ab}=ϵ_{abc}e^c`$ is the projected alternating tensor. The covariant derivative $`\stackrel{~}{}_a`$ is obtained from the 3-dimensional (flat) covariant derivative $`_a`$ by total projection:

$$\stackrel{~}{}_aT_{b\mathrm{}c}=\stackrel{~}{g}_a^d\stackrel{~}{g}_b^e\mathrm{}\stackrel{~}{g}_c^f_dT_{e\mathrm{}f},$$
(A4)
for an arbitrary tensor $`T_{b\mathrm{}c}`$. Making use of the results

$`(\stackrel{~}{}_aY_{(lm)}^{})|_{e^c}`$ $`=`$ $`\mathrm{\Lambda }_a{}_{}{}^{b}(\stackrel{~}{}_bY_{(lm)})|_{e^d\mathrm{\Lambda }_d^c}`$ (A5)
$`(\stackrel{~}{}_a\stackrel{~}{}_bY_{(lm)}^{})|_{e^c}`$ $`=`$ $`\mathrm{\Lambda }_a{}_{}{}^{d_1}\mathrm{\Lambda }_{b}^{}{}_{}{}^{d_2}(\stackrel{~}{}_{d_1}\stackrel{~}{}_{d_2}Y_{(lm)})|_{e^{d_3}\mathrm{\Lambda }_{d_3}^c},`$ (A6)
it is straightforward to prove that the rotated tensor harmonics are obtained by replacing the scalar harmonics by their rotated counterparts $`Y_{(lm)}^{}(e^c)=Y_{(lm)}(e^d\mathrm{\Lambda }_d{}_{}{}^{c})`$ in Eqs. (A2) and (A3). Since the $`l`$th scalar harmonics transform irreducibly under rotations as (e.g. Ref. )

$$Y_{(lm)}^{}(𝒆)=\underset{|m^{}|l}{}D_{m^{}m}^l(\varphi ,\theta ,\psi )Y_{(lm^{})}(𝒆),$$
(A7)
the tensor harmonics inherit the same transformation law:

$$Y_{(lm)}^\text{P}^{}(𝒆)=\underset{|m^{}|l}{}D_{m^{}m}^l(\varphi ,\theta ,\psi )Y_{(lm^{})}^\mathrm{P}(𝒆).$$
(A8)


text:
# Explaining away ambiguity: Learning verb selectional preference with Bayesian networksWe would like to thank the Brown Laboratory for Linguistic Information Processing; Thomas Hofmann; Elie Bienenstock; Philip Resnik, who provided us with training and test data; and Daniel Garcia for his help with the SMILE library of classes for Bayesian networks that we used for our experiments. This research was supported by NSF awards 9720368, 9870676 and 9812169.

## 1 Selectional preference and sense ambiguity

Regularities of a verb with respect to the semantic class of its arguments (subject, object and indirect object) are called selectional preferences (SP) \[Katz and Fodor, 1964, Chomsky, 1965, Johnson-Laird, 1983\]. The verb pilot carries the information that its object will likely be some kind of vehicle; subjects of the verb think tend to be human; and subjects of the verb bark tend to be dogs. For the sake of simplicity we will focus on the verb-object relation although the techniques we will describe can be applied to other verb-argument pairs.

Models of the acquisition of SP are important in their own right and have applications in Natural Language Processing (NLP). The selectional preferences of a verb can be used to infer the possible meanings of an unknown argument of a known verb; e.g., it might be possible to infer that xxxx is a kind of dog from the following sentence: “The *xxxx* barked all night”. In parsing a sentence selectional preferences can be used to rank competing parses, providing a partial measure of semantic well-formedness. Investigating SP might help us to understand the structure of the mental lexicon.

Systems for unsupervised learning of SP usually combine statistical and knowledge-based approaches. The knowledge-base component is typically a database that groups words into classes. In the models we will see, the knowledge base is Wordnet \[Miller, 1990\]. Wordnet groups nouns into classes of synonyms representing concepts, called synsets, e.g., $`\{car,auto,automobile,\mathrm{}\}`$. A noun that belongs to several synsets is ambiguous. A transitive and asymmetrical relation, hyponymy, is defined between synsets. A synset is a hyponym of another synset if the former has the latter as a broader concept; for example, BEVERAGE is a hyponym of LIQUID. Figure 1 depicts a portion of the hierarchy.

The statistical component consists of predicate-argument pairs extracted from a corpus in which the semantic class of the words is not indicated. A trivial algorithm might get a list of words that occurred as objects of the verb and output the semantic classes the words belong to according to Wordnet. For example, if the verb drink occurred with water and $`waterLIQUID`$, the model would learn that drink selects for LIQUID. As ?) and ?) have found, the main problem these systems face is the presence of ambiguous words in the training data. If the word java also occurred as an object of drink, since $`javaBEVERAGE`$ and $`javaISLAND`$, this model would learn that drink selects for both $`BEVERAGE`$ and $`ISLAND`$.

More complex models have been proposed. These models, though, deal with word sense ambiguity by applying an unselective strategy similar to the one above; i.e., they assume that ambiguous words provide equal evidence for all their senses. These models choose as the concepts the verb selects for those that are in common among several words (e.g., BEVERAGE above). This strategy works to the extent that these overlapping senses are also the concepts the verb selects for.

## 2 Previous approaches to learning selectional preference

### 2.1 Resnik’s model

Ours system is closely related to those proposed in \[Resnik, 1997\] and \[Abney and Light, 1999\]. The fact that a predicate p selects for a class c, given a syntactic relation r, can be represented as a relation, $`selects(p,r,c)`$; e.g., that eat selects for FOOD in object position can be represented as $`selects(eat,object,FOOD)`$. In \[Resnik, 1997\] selectional preference is quantified by comparing the prior distribution of a given class c appearing as an argument, $`P(c)`$, and the conditional probability of the same class given a predicate and a syntactic relation $`P(c|p,r)`$, e.g., $`P(FOOD)`$ and $`P(FOOD|eat,object)`$. The relative entropy between $`P(c)`$ and $`P(c|p,r)`$ measures how much the predicate constrains its arguments:

$$S(p,r)=D(P(c|p,r)||P(c))$$
(1)
Resnik defines the selectional association of a predicate for a particular class c to be the portion of the selectional preference strength due to that class:

$$A(p,r,c)=\frac{1}{S(p,r)}P(c|p,r)\mathrm{log}\frac{P(c|p,r)}{P(c)}$$
(2)
Here the main problem is the estimation of $`P(c|p,r)`$. Resnik suggests as a plausible estimator $`\widehat{P}(c|p,r)\stackrel{\mathrm{def}}{=}freq(p,r,c)/freq(p,r).`$ But since the model is trained on data that are not sense-tagged, there is no obvious way to estimate $`freq(p,r,c)`$. Resnik suggests considering each observation of a word as evidence for each of the classes the word belongs to,

$$freq(p,r,c)\underset{wc}{}\frac{count(p,r,w)}{classes(w)}$$
(3)
where $`count(p,r,w)`$ is the number of times the word w occurred as an argument of p in relation r, and $`classes(w)`$ is the number of classes w belongs to.

For example, suppose the system is trained on (eat,object) pairs and the verb occurred once each with meat, apple, bagel, and cheese, and Wordnet is simplified as in Figure 2. An ambiguous word like meat provides evidence also for classes that appear unrelated to those selected by the verb. Resnik’s assumption is that only the classes selected by the verb will be associated with each of the observed words, and hence will receive the highest values for $`P(c|p,r)`$. Using (3) we find that the highest frequency is in fact associated with FOOD: $`freq(eat,object,food)\frac{1}{4}+\frac{1}{2}+\frac{1}{2}+\frac{1}{2}=\frac{7}{4}`$ and $`P(FOOD|eat)=0.44`$. However, some evidence is found also for COGNITION: $`freq(eat,object,cognition)\frac{1}{4}`$ and $`P(COGNITION|eat)=0.06`$.

### 2.2 Abney and Light’s approach

?) pointed out that the distribution of senses of an ambiguous word is not uniform. They noticed also that it is not clear how the probability $`P(c|p,r)`$ is to be interpreted since there is no explicit stochastic generation model involved.

They proposed a system that associates a Hidden Markov Model (HMM) with each predicate-relation pair (p,r). Transitions between synset states represent the hyponymy relation, and $`\epsilon `$, the empty word, is emitted with probability 1; transitions to a final state emit a word w with probability $`0P(w)1`$. Transition and emission probabilities are estimated using the EM algorithm on training data that consist of the nouns that occurred with the verb. Abney and Light’s model estimates $`P(c|p,r)`$ from the model trained for $`(p,r)`$; the distribution $`P(c)`$ can be calculated from a model trained for all nouns in the corpus.

This model did not perform as well as expected. An ambiguous word in the model can be generated by more than one state sequence. Abney and Light discovered that the EM algorithm finds parameter values that associate some probability mass with all the transitions in the multiple paths that lead to an ambiguous word. In other words, when there are several state sequences for the same word, EM does not select one of them over the others.111As a matter of fact, for this HMM there are (infinitely) many parameter values that maximize the likelihood of the training data; i.e., the parameters are not identifiable. The intuitively correct solution is one of them, but so are infinitely many other, intuitively incorrect ones. Thus it is no surprise that the EM algorithm cannot find the intuitively correct solution. Figure 3 shows the parameters estimated by EM for the same example as above. The transition to the COGNITION state has been assigned a probability of $`1/8`$ because it is part of a possible path to meat. The HMM model does not solve the problem of the unselective distribution of the frequency of occurrence of an ambiguous word to all its senses. Abney and Light claimed that this is a serious problem, particularly when the ambiguous word is a frequent one, and caused the model to learn the wrong selectional preferences. To correct this undesirable outcome they introduced some smoothing and balancing techniques. However, even with these modifications their system’s performance was below that achieved by Resnik’s.

## 3 Bayesian networks

A Bayesian network \[Pearl, 1988\], or Bayesian belief network (BBN), consists of a set of variables and a set of directed edges connecting the variables. The variables and the edges define a directed acyclic graph (DAG) where each variable is represented by a node. Each variable is associated with a finite number of (mutually exclusive) states. To each variable $`A`$ with parents $`B_1,\mathrm{},B_n`$ is attached a conditional probability table (CPT) $`P(A|B_1,\mathrm{},B_n)`$. Given a BBN, Bayesian inference can be used to estimate marginal and posterior probabilities given the evidence at hand and the information stored in the CPTs, the prior probabilities, by means of Bayes’ rule, $`P(H|E)=\frac{P(H)P(E|H)}{P(E)}`$, where H stands for hypothesis and E for evidence.

Bayesian networks display an extremely interesting property called explaining away. Word sense ambiguity in the process of learning SP defines a problem that might be solved by a model that implements an explaining away strategy. Suppose we are learning the selectional preference of drink, and the network in Figure 4 is the knowledge base. The verb occurred with java and water. This situation can be represented as a Bayesian network. The variables *ISLAND* and *BEVERAGE* represent concepts in a semantic hierarchy. The variables *java* and *water* stand for possible instantiations of the concepts. All the variables are Boolean; i.e., they are associated with two states, true or false. Suppose the following CPTs define the priors associated with each node.222I, B, j and w abbreviate ISLAND, BEVERAGE, java and water, respectively. The unconditional probabilities are $`P(I=true)=P(B=true)=0.01`$ and $`P(I=false)=P(B=false)=0.99`$, and the CPTs for the child nodes are

|  | $`P(X=x|Y_1=y_1,Y_2=y_2)`$ | | | |
| --- | --- | --- | --- | --- |
|  | $`I,B`$ | $`I,\neg B`$ | $`\neg I,B`$ | $`\neg I,\neg B`$ |
| j = true | 0.99 | 0.99 | 0.99 | 0.01 |
| j = false | 0.01 | 0.01 | 0.01 | 0.99 |
| w = true | 0.99 | 0.99 | 0.01 | 0.01 |
| w = false | 0.01 | 0.01 | 0.99 | 0.99 |

These values mean that the occurrence of either concept is *a priori* unlikely. If either concept is true the word *java* is likely to occur. Similarly, if *BEVERAGE* occurs it is likely to observe also the word *water*. As the posterior probabilities show, if *java* occurs, the beliefs in both concepts increase: $`P(I|j)=P(B|j)=0.3355`$. However, *water* provides evidence for *BEVERAGE* only. Overall there is more evidence for the hypothesis that the concept being expressed is *BEVERAGE* and not *ISLAND*. Bayesian networks implement this inference scheme; if we compute the conditional probabilities given that both words occurred, we obtain $`P(B|j,w)=0.98`$ and $`P(I|j,w)=0.02`$. The new evidence caused the “island” hypothesis to be explained away!

### 3.1 The relevance of priors

Explaining away seems to depend on the specification of the prior probabilities. The priors define the background knowledge available to the model relative to the conditional probabilities of the events represented by the variables, but also about the joint distributions of several events. In the simple network above, we defined the probability that either concept is selected (i.e., that the corresponding variable is true) to be extremely small. Intuitively, there are many concepts and the probability of observing any particular one is small. This means that the joint probability of the two events is much higher in the case in which only one of them is true (0.0099) than in the case in which they are both true (0.0001). Therefore, via the priors, we introduced a bias according to which the hypothesis that one concept is selected will be favored over two co-occurring ones. This is a general pattern of Bayesian networks; the prior causes simpler explanations to be preferred over more complex ones, and thereby the explaining away effect.

## 4 A Bayesian network approach to learning selectional preference

### 4.1 Structure and parameters of the model

The hierarchy of nouns in Wordnet defines a DAG. Its mapping into a BBN is straightforward. Each word or synset in Wordnet is a node in the network. If A is a hyponym of B there is an arc in the network from B to A. All the variables are Boolean. A synset node is true if the verb selects for that class. A word node is true if the word can appear as an argument of the verb. The priors are defined following two intuitive principles. First, it is unlikely that a verb *a priori* selects for any particular synset. Second, if a verb does select for a synset, say FOOD, then it is likely that it also selects for its hyponyms, say FRUIT. The same principles apply to words: it is likely that a word appears as an argument of the verb if the verb selects for any of its possible senses. On the other hand, if the verb does not select for a synset, it is unlikely that the words instantiating the synset occur as its arguments. “Likely” and “unlikely” are given numerical values that sum up to 1. The following table defines the scheme for the CPTs associated with each node in the network; $`p_i(X)`$ denotes the $`ith`$ parent of the node $`X`$.

|  | $`P(X=x|p_1(X),\mathrm{},p_n(X)=true)`$ |
| --- | --- |
| x = true | likely |
| x = false | unlikely |
|  | $`P(X=x|p_1(X),\mathrm{},p_n(X)=false)`$ |
| x = true | unlikely |
| x = false | likely |

For the root nodes, the table reduces to the unconditional probability of the node. Now we can test the model on the simple example seen earlier. $`W^+`$ is the set of words that occurred with the verb. The nodes corresponding to the words in $`W^+`$ are set to true and the others left unset. For the previous example $`W^+=\{meat,apple,bagel,cheese\}`$, and the corresponding nodes are set to true, as depicted in Figure 5. With likely and unlikely respectively equal to 0.99 and 0.01, the posterior probabilities are333F, C, m, a, b and c respectively stand for FOOD, COGNITION, meat, apple, bagel and cheese $`P(F|m,a,b,c)=0.9899`$ and $`P(C|m,a,b,c)=0.0101`$. Explaining away works. The posterior probability of COGNITION gets as low as its prior, whereas the probability of FOOD goes up to almost 1. A Bayesian network approach seems to actually implement the conservative strategy we thought to be the correct one for unsupervised learning of selectional restrictions.

### 4.2 Computational issues in building BBNs based on Wordnet

The implementation of a BBN for the whole of Wordnet faces computational complexity problems typical of graphical models. A densely connected BBN presents two kinds of problems. The first is the storage of the CPTs. The size of a CPT grows exponentially with the number of parents of the node.444Some words in Wordnet have more than 20 senses. For example, line in Wordnet is associated with 25 senses. The size of its CPT is therefore $`2^{26}`$. Storing a table of float numbers for this node alone requires around $`(2^{26})8=537`$ MBytes of memory. This problem can be solved by optimizing the representation of these tables. In our case most of the entries have the same values, and a compact representation for them can be found (much like the one used in the noisy-OR model \[Pearl, 1988\]).

A harder problem is performing inference. The graphical structure of a BBN represents the dependency relations among the random variables of the network. The algorithms used with BBNs usually perform inference by dynamic programming on the triangulated moral graph. A lower bound on the number of computations that are necessary to model the joint distribution over the variables using such algorithms is $`2^{|n|+1}`$, where $`n`$ is the size of the maximal boundary set according to the visitation schedule.

### 4.3 Subnetworks and balancing

Because of these problems we could not build a single BBN for Wordnet. Instead we simplified the structure of the model by building a smaller subnetwork for each predicate-argument pair. A subnetwork consists of the union of the sets of ancestors of the words in $`W^+`$. Figure 6 provides an example of the union of these “ancestral subgraphs” of Wordnet for the words java and drink (compare it with Figure 1).

This simplification does not affect the computation of the distributions we are interested in; that is, the marginals of the synset nodes. A BBN provides a compact representation for the joint distribution over the set of variables in the network. If $`N=X_1,\mathrm{},X_n`$ is a Bayesian network with variables $`X_1,\mathrm{},X_n`$, its joint distribution $`P(N)`$ is the product of all the conditional probabilities specified in the network,

$$P(N)=\underset{j}{}P(X_j|pa(X_j))$$
(4)
where $`pa(X)`$ is the set of parents of $`X`$. A BBN generates a factorization of the joint distribution over its variables. Consider a network of three nodes $`A,B,C`$ with arcs from $`A`$ to $`B`$ and $`C`$. Its joint distribution can be characterized as $`P(A,B,C)=P(A)P(B|A)P(C|A)`$. If there is no evidence for $`C`$ the joint distribution is

$`P(A,B,C)`$ $`=`$ $`P(A)P(B|A){\displaystyle \underset{C}{}}P(C|A)`$ (5)
$`=`$ $`P(A)P(B|A)`$
$`=`$ $`P(A,B)`$
The node $`C`$ gets marginalized out. Marginalizing over a childless node is equivalent to removing it with its connections from the network. Therefore the subnetworks are equivalent to the whole network; i.e., they have the same joint distribution.

Our model computes the value of $`P(c|p,r)`$, but we did not compute the prior $`P(c)`$ for all nouns in the corpus. We assumed this to be a constant, equal to the unlikely value, for all classes. In a BBN the values of the marginals increase with their distance from the root nodes. To avoid undesired bias (see table of results) we defined a balancing formula that adjusted the conditional probabilities of the CPTs in such a way that we got all the marginals to have approximately the same value.555More details can be found in an extended version of the paper: www.cog.brown.edu/$``$massi/.

## 5 Experiments and results666For these experiments we used values for the likely and unlikely parameters of 0.9 and 0.1, respectively.

### 5.1 Learning of selectional preferences

When trained on predicate-argument pairs extracted from a large corpus, the San Jose Mercury Corpus, the model gave very good results. The corpus contains about 1.3 million verb-object tokens. The obtained rankings of classes according to their posterior marginal probabilities were good. Table 1 shows the top and the bottom of the list of synsets for the verb maneuver. The model learned that maneuver “selects” for members of the class VEHICLE and of other plausible classes, hyponyms of VEHICLE. It also learned that the verb does not select for direct objects that are members of classes, like CONCEPT or PHILOSOPHY.

### 5.2 Word sense disambiguation test

A direct evaluation measure for unsupervised learning of SP models does not exist. These models are instead evaluated on a word-sense disambiguation test (WSD). The idea is that systems that learn SP produce word sense disambiguation as a side-effect. Java might be interpreted as the island or the beverage, but in a context like “the tourists flew to Java” the former seems more correct, because fly could select for geographic locations but not for beverages. A system trained on a predicate p should be able to disambiguate arguments of *p* if it has learned its selectional restrictions.

We tested our model using the test and training data developed by Resnik (see Resnik, 1997). The same test was used in \[Abney and Light, 1999\]. The training data consists of predicate-object counts extracted from 4/5 of the Brown corpus (about 1M words). The test set consists of predicate-object pairs from the remaining 1/5 of the corpus, which has been manually sense-annotated by Wordnet researchers. The results are shown in Table 2. The baseline algorithm chooses at random one of the multiple senses of an ambiguous word. The “first sense” method always chooses the most frequent sense (such a system should be trained on sense-tagged data).

Our model performed better than the state of the art models for unsupervised learning of SP. It seems to define a better estimator for $`P(c|p,r)`$.

It is remarkable that the model achieved this result making only a limited use of distributional information. A noun is in $`W^+`$ if it occurred at least once in the training set, but the system does not know if it occurred once or several times; either it occurred or it didn’t. The model did not suffer too much from this limitation during this task. This is probably due to the sparseness of the training data for the test. For each verb the average number of object types is $`3.3`$, for each of them the average number of tokens is $`1.3`$; i.e., most of the words in the training data only occurred once. For this training set we also tested a version of the model that built a word node for each observed object token and therefore integrated the distributional information. On the WSD test it performed exactly the same as the simpler version. When trained on the San Jose Mercury Corpus the model performed worse on the WSD test ($`35.8\%`$). This is not too surprising considering the differences between the SJM and the Brown corpora: the former, a recent newswire corpus; the latter, an older, balanced corpus. Another important factor is the different relevance of distributional information. The training data from the SJM Corpus are much richer and noisier than the Brown data. Here the frequency information is probably crucial; however, in this case we could not implement the simple scheme above.

### 5.3 Conclusion

Explaining away implements a cognitively attractive and successful strategy. A straightforward improvement would be for the model to make full use of the distributional information present in the training data; we only partially achieved this. Bayesian networks are usually confronted with a single presentation of evidence. Their extension to multiple evidence is not trivial. We believe the model can be extended in this direction. Possibly there are several ways to do so (multinomial sampling, dedicated implementations, etc.). However, we believe that the most relevant finding of this research might be that “explaining away” is not only a property of Bayesian networks but of Bayesian inference in general and that it might be implementable in other kinds of graphical models. We observed that the property seems to depend on the specification of the prior probabilities. We found that the HMM model of \[Abney and Light, 1999\] was unidentifiable; that is, there are several solutions for the parameters of the model, including the desired one. Our intuition is that it should be possible to implement “explaining away” in a HMM with priors, so that it would prefer only one or a few solutions over many. This model would have also the advantage of being computationally simpler.


text:
# Semiempirical Shell Model Masses with Magic Number 𝐙=𝟏𝟐𝟔 for Translead Elements with 𝐍≤𝟏𝟐𝟔

## Abstract

A semiempirical shell model mass equation based on magic number $`Z=126`$ and applicable to translead elements with $`N126`$ is presented. For $`\alpha `$-decay energies the equation is shown to have a high predictive power and an rms deviation from the data of about 100 keV. The rms deviations for masses and other mass differences are between about 200 and 300 keV.

PACS numbers: 21.10.Dr, 21.60.Cs, 27.80.+w

Recent progress in superheavy elements (SHE) research reaching to $`{}_{}{}^{293}118`$ and its $`\alpha `$-decay products  makes it necessary to find an appropriate substitute for the semiempirical shell-model mass equation (SSME)  (see also ref. ) for nuclei in the neighbourhood of $`Z=114`$ and beyond . The $`\alpha `$ energies of the decaying chain vary smoothly from $`{}_{}{}^{293}118`$ to 269Sg $`(Z=106)`$, with no indication of magicity at $`Z=114`$ in these nuclei , whereas the SSME assumes that $`Z=114`$ is the next spherical proton magic number after lead and it stops there. Furthermore, the SSME becomes unsuitable for extrapolation already earlier, beyond Hs $`(Z=108)`$, as shown by its increasing deviations from the data when $`Z`$ increases. (Like in fig. 4 of ref. .)

Recent phenomenological studies of BE(2) systematics  and of the persistence of the Wigner term in masses of heavy nuclei  indicate $`Z=126`$ as the next spherical proton magic number after lead, and this is consistent with considerations based on nuclear diffuseness . Recent self-consistent and relativistic mean-field calculations  variously predict proton magicities for $`Z=114,120,124`$ and 126.

During the early stages of the SSME , when it was adjusted separately in individual shell regions in the $`NZ`$ plane, both $`Z=114`$ and $`Z=126`$, which were at the time considered possible candidates for the postlead proton magic number (see, e.g., ref. ), were tried as a shell region boundary in each of the two heaviest regions with $`Z82`$ and respective $`N`$-boundaries $`82N126`$ (called here region A) and $`126N184`$ (called region B). The agreement with the data was about the same for both choices, and the prevailing view in the mid nineteen-seventies led to the choice of $`Z=114`$ for the SSME mass table . In a recent communication  we showed that the early $`Z=126`$ results have a high predictive power in the interior of region B and proposed their use there as a predictive tool in SHE research. In the present note we study the predictive power or extrapolatability of the early $`Z=126`$ results in region A and propose a mass equation based on them as a substitute for the SSME  in the interior of the region. The study is based on the newer data that became available after the adjustments were made, as was done in refs. .

In the SSME the total nuclear energy $`E`$ is a sum of pairing, deformation and Coulomb energies:

$$E(N,Z)=E_{pair}(N,Z)+E_{def}(N,Z)+E_{Coul}(N,Z).$$
(1)
The form of $`E_{Coul}`$ is the same in all shell regions:

$$E_{Coul}(N,Z)=\left(\frac{2Z_0}{A}\right)^{1\mathrm{/}3}[\alpha ^c+\beta ^c\left(ZZ_0\right)+\gamma ^c\left(ZZ_0\right)^2],$$
(2)
and that of $`E_{pair}`$ is the same separately in all diagonal shell regions, where the major valence shells are the same for neutrons and for protons, and in all non-diagonal regions, where the neutron and proton valence shells are different. Unlike in , with $`Z=126`$ rather than 114 as an upper proton boundary, region A becomes a diagonal region with

$`E_{pair}(N,Z)=`$
$`\left({\displaystyle \frac{A_0}{A}}\right)\left[\alpha +\beta \left(AA_0\right)+\gamma \left(AA_0\right)^2+\epsilon T\left(T+1\right)+{\displaystyle \frac{1\left(1\right)^A}{2}}\mathrm{\Theta }+{\displaystyle \frac{1\left(1\right)^{NZ}}{2}}\kappa \right].`$ (3)
The part $`E_{def}`$ for region A with $`Z=126`$ as upper proton boundary is

$$E_{def}(N,Z)=\left(\frac{A_0}{A}\right)\left[\phi _{11}\mathrm{\Phi }_{11}(N,Z)+\psi _{20}\left[\mathrm{\Psi }_{20}(N,Z)+\mathrm{\Psi }_{20}(Z,N)\right]\right]$$
(4)
with

$$\mathrm{\Phi }_{11}(N,Z)=\left(N82\right)\left(126N\right)\left(Z82\right)\left(126Z\right)$$
(5)
$$\mathrm{\Psi }_{20}(N,Z)=\left(N82\right)^2\left(126N\right)^2\left(N104\right).$$
(6)
In eqs. (2)$``$(4) $`A=N+Z`$ and $`T=|T_z|=\frac{1}{2}|NZ|`$222In the as yet unknown odd-odd $`N=Z`$ translead nuclei the g.s. is expected to have $`T=|T_z|+1`$ and seniority zero, whereas eq. (1) with $`T=|T_z|`$ gives the energy of a low excited seniority two state . (See also ref. .). The respective values of $`Z_0`$ and $`A_0`$ are 82 and 164. The coefficients multiplying the functions of $`N`$ and $`Z`$ are adjustable parameters which were determined by a least-squares adjustment to the data . Their numerical values resulting from that adjustment are given in the second column of table I. The mass excesses $`\mathrm{\Delta }M(N,Z)`$ are obtained by adding to eq. (1) the sum of nucleon mass excesses $`N\mathrm{\Delta }M_n+Z\mathrm{\Delta }M_H`$.

The experimental data used in the adjustment included 29 masses and 62 $`Q_\alpha `$ values connecting unknown masses (Ref.  augmented by data from the literature up to Spring 1973). Presently there are 150 known experimental masses in region A and 3 $`Q_\alpha `$ values connecting unknown masses (Refs.  (excluding values denoted “systematics” (#)) and , augmented by data from the recent literature). There are 121 new masses that were not used in the adjustments.

Fig. 1 shows the deviations from the data of the predictions of eq. (1) (with the definitions (2)$``$(4) and the coefficients from the second column of table I) for all the 150 known masses. The deviations are plotted as function of the distance from the line of $`\beta `$-stability, denoted “neutrons from stability” ($`NFS`$) and defined by $`NFS=NZ\frac{0.4A^2}{A+200}`$ . Empty circles denote the deviations of the 29 originally adjusted masses and full circles mark the deviations of the 121 new masses. The original deviations are relatively small. The new deviations are mostly considerably larger and are almost all negative. There is a large scatter of the points, superposed on an overall decreasing trend when NFS increases in the negative direction. For closer scrutiny we show in fig. 2 lines connecting deviations of Pb, Bi, Po and At isotopes as function of $`T`$ rather than $`NFS`$. Isotopic lines of other elements follow closely. There is an odd-even staggering of the points, where for even-even and odd-odd nuclei the deviations are respectively higher and lower than for odd-mass ones, and there is an overall increasing overbinding when $`T`$ decreases away from stability. The increasing staggering when $`T`$ decreases as compared to the old data indicates increasing pairing parameters $`\mathrm{\Theta }`$ (see also ref. ) and (even more so) $`\kappa `$ towards the proton drip line. Similarly, the increasing binding when $`T`$ decreases indicates that the symmetry-energy parameter $`\epsilon `$ decreases towards the drip line.

Table II, patterned after similar more elaborate ones  (see also ref. ), shows the values of $`\delta _{av}`$ and $`\delta _{rms}`$, the respective average and rms deviations of eq. (1) from the data, for $`\mathrm{\Delta }M,S_n,S_p,Q_\beta ^{}`$ and $`Q_\alpha `$. The deviations are shown separately for the older data that were used in the adjustment and for the newer data. The last column shows the error ratios $`\delta _{rms}^{new}:\delta _{rms}^{old}`$.

For the older data the $`\delta _{av}`$ are few tens keV at most, and the $`\delta _{rms}`$ are in the range 100-250 keV. For the newer nuclei the agreement with the data of predicted $`\mathrm{\Delta }M`$ (fig. 1), $`S_n,S_p`$ and $`Q_\beta ^{}`$ values has much deteriorated. On the other hand, the $`\delta _{rms}`$ of the predicted $`Q_\alpha `$ values has even very slightly improved. The $`Q_\alpha `$ deviations, both old and new, are remarkably small.

The high degree of extrapolatability of the $`Q_\alpha `$ values as compared to the poorer predictions of $`\mathrm{\Delta }M,S_n,S_p`$ and $`Q_\beta ^{}`$ is presumably due to the composition of the old data set used in the adjustments. The coefficients $`\beta ,\gamma ,\beta ^c,\gamma ^c,\phi _{11}`$ and $`\psi _{20}`$, which contribute to both mass and $`Q_\alpha `$ values, and also $`\alpha ^c`$, which was determined using $`\beta ^c`$ and $`\gamma ^c`$  (see also ref. ), were determined by all the 91 available mass and $`Q_\alpha `$ data. On the other hand, neglecting their $`A`$-dependence the coefficients $`\alpha ,\epsilon ,\mathrm{\Theta }`$ and $`\kappa `$ cancel in $`Q_\alpha `$ and they were determined essentially by the smaller group of 29 masses found in the nearest-to-stability corner of the shell region, with too-small values of $`\mathrm{\Theta }`$ and $`\kappa `$ and too large value of $`\epsilon `$ as compared to more proton-rich nuclei. These values, which to a large extent cancel in $`Q_\alpha `$, are responsible for the large deviations of the predicted new masses.

One would like to restore to the new $`\mathrm{\Delta }M,S_n,S_p`$ and $`Q_\beta ^{}`$ predictions the same degree of agreement with the data as the old predictions had, while at the same time retaining the high quality of $`Q_\alpha `$ predictions. For eq. (1) this goal might most simply be approached by making a least-squares adjustment of eq. (1) to all the 150 known masses, with only four adjustable parameters $`\alpha ,\epsilon ,\mathrm{\Theta }`$ and $`\kappa `$, while the other seven coefficients are held fixed on their old values from the second column of table I. The re-adjusted parameters $`\alpha `$ and $`\epsilon `$ would correct the systematic overbinding when $`T`$ decreases, and the re-adjusted $`\mathrm{\Theta }`$ and $`\kappa `$ would decrease the odd-even staggering of the deviations.

The re-adjusted values of the coefficients $`\alpha ,\epsilon ,\mathrm{\Theta }`$ and $`\kappa `$ are given in the third column of table I. Their changes compared to the second column are in the anticipated directions. The respective average and rms mass deviations obtained in the adjustment for all the 150 mass data are 2 and 246 keV.

Fig. 3 shows the deviations from the data of the predicted mass values resulting from the new adjustment, similarly to fig. 1. For ease of comparison the same nuclei are denoted by the same kind of circle (empty or full) in the two figures. The deviations for the post 1973 measured nuclei (denoted by full circles) and their odd-even staggering are considerably smaller than in fig. 1, and the overall decreasing trend when $`T`$ decreases has largely disappeared, as was expected. The deviations for the older data have worsened, though, and there is a new overall oscillatory trend not observed in fig. 1, indicating that enlarging eq. (1) by additional particle-hole symmetric $`E_{def}`$ terms (eq. (4)) might reduce the deviations further.

Table III shows the resulting $`\delta _{av}`$ and $`\delta _{rms}`$ values of the predicted deviations. Like in table II they are shown separately for the old 1973 data nearer to stability and for the new data extending into the interior of region A towards the proton drip line. For $`S_n,S_p`$ and $`Q_\alpha `$ the new deviations are similar to those of the original nuclei in table II. For $`\mathrm{\Delta }M`$ and $`Q_\beta ^{}`$ the $`\delta _{rms}`$ are respectively 1.5 and 1.1 times larger.

The new mass deviations shown in table III are about one half to two thirds of the corresponding deviations of several recent mass models , and for $`Q_\alpha `$ they are smaller. The smaller values of the deviations in table III are presumably due mainly to the inclusion of the particle-hole symmetric configuration-interaction terms $`E_{def}`$ (eq. (4)) in eq. (1) . Until a new adjustment of the $`Z=126`$ SSME to the data in both regions A and B is undertaken, which would further reduce the deviations in both regions, we propose the use of eq. (1), with the new values of $`\alpha ,\epsilon ,\mathrm{\Theta }`$ and $`\kappa `$ and the old values of the other coefficients, given in table I, as a predictive tool for masses and their differences , and particularly for the extrapolatable-proven $`Q_\alpha `$ values, in the interior of region A.

It should be emphasized, though, that the above results are not a proof of superior magicity of $`Z=126`$ in region A as compared to other recently proposed prediction , because no comparative studies of this kind were made. (See also ref. .)

We thank the Atomic Mass Data Center in Orsay for ref.  and Stelian Gelberg and Dietmar Kolb for help with the calculations.


text:
# Simulations on Infinite Size Lattices

\[

## Abstract

We introduce a Monte Carlo method, as a modification of existing cluster algorithms, which allows simulations directly on systems of infinite size, and for quantum models also at $`\beta =\mathrm{}`$. All two-point functions can be obtained, including dynamical information. When the number of iterations is increased, correlation functions at larger distances become available. Limits $`q0`$ and $`\omega 0`$ can be approached directly. As examples we calculate spectra for the d=2 Ising model and for Heisenberg quantum spin ladders with 2 and 4 legs.

\]

Standard Monte Carlo simulations are limited to systems of finite size. Physical results for infinite systems have to be obtained by finite size scaling, assuming that one knows the correct scaling laws, and assuming that the data are already in a suitable asymptotic regime. It is therefore very desirable to obtain results also directly at infinite system size. We will show how to do so, with only a small modification of existing cluster algorithms, by using the cluster representation of the models to calculate two-point functions. In the quantum case we can then also simulate directly at $`\beta =\mathrm{}`$ and calculate correlation functions and dynamical greens functions. As examples we will show calculations for the classical Ising model and for quantum Heisenberg ladder systems.

The Swendsen Wang cluster method  for the classical Ising model is based on the Edwards-Sokal-Fortuin-Kasteleyn  representation

$`e^{\beta J(s_is_j\mathrm{\hspace{0.33em}1})}`$ $`=`$ $`{\displaystyle \underset{b_{ij}=0,1}{}}p\delta _{s_is_j}\delta _{b_{ij},1}+(1p)\delta _{b_{ij},0}`$ (1)
$`=:`$ $`{\displaystyle \underset{b_{ij}=0,1}{}}W_{ij}(s_i,s_j,b_{ij})`$ (2)
for the Boltzmann weight of a spin-pair, with $`p=1e^{2\beta J}`$, which enlargens the phase space of spin variables $`s_i`$ by additional bond variables $`b_{ij}`$. The partition function $`Z=_{\{s_i\}}_{\{b_{ij}\}}W(s,b)`$ with total weight $`W(s,b)=_{ij}W_{ij}(s_i,s_j,b_{ij})`$ is then simulated efficiently by switching between representations: given a spin-configuration $`s:=\{s_i\}`$, one generates a bond configuration $`b:=\{b_{ij}\}`$ with the conditional probability $`p(b|s)W(s,b)`$, thus creating a configuration of clusters of sites connected by bonds $`b_{ij}=1`$. Given a bond configuration, a new spin configuration is generated with probability $`p(s|b)W(s,b)`$. Because of the factor $`\delta _{s_is_j}\delta _{b_{ij},1}`$ in $`W_{ij}`$, this amounts to setting all spins of each cluster randomly to a common new value, independent of other clusters.

Observables $`\widehat{O}`$ can be computed either in spin-representation as $`O(s)`$ or in bond representation as $`O(b)=\left(_sO(s)W(s,b)\right)/\left(_sW(s,b)\right)`$ (so called “improved estimators”) . For two-point functions $`O(s)=s_is_j`$ the bond representation is particularly simple:

$$O(b)=\delta (\text{sites }i\text{ and }j\text{ are in the same cluster}).$$
(3)
Thus two-point functions and derived quantities (including susceptibility, energy, specific heat) can be computed from the properties of individual clusters.

A variant of the Swendsen-Wang method is Wolff’s single cluster method . Given a spin configuration, only one of the bond-clusters is constructed, namely the cluster which contains a randomly chosen intial starting site $`x_0`$. All spins in this cluster are then flipped. One advantage is that the single cluster will on average be larger than the average Swendsen-Wang cluster, so that its flip results in a bigger move in phase space. Using eq. (3), the correlation function can then be measured as

$$C(r)=s_is_{i+r}=\frac{1}{V_{cl}}\underset{i\text{ in cluster}}{}\delta (i+r\text{ in cluster}),$$
(4)
where $`V_{cl}`$ is the number of sites in the single cluster, and the brackets $`\mathrm{}`$ on the right refer to the Monte Carlo average. The susceptibility is then $`\chi =\beta V_{cl}`$, the average energy can be measured either as $`E=J\mathrm{\hspace{0.17em}2}dC(1)`$ or from the average number of bonds in the cluster , and the specific heat can be obtained either as a numerical derivative or directly from the bond representation .

Our new infinite system method now employs the single cluster method, except that it starts each cluster not from a randomly chose site, but always from the same site $`x_0`$ (e.g. the origin of the coordinate system). Our method satisfies detailed balance and ergodicity (for any finite region including the origin) similarly to Wolff’s method. For the ferromagnetic Ising model, we begin with an initial staggered spin configuration of unlimited size. (Only a finite part if this configuration will have to be stored). We iterate the following two steps: (1) For the current spin configuration, a cluster containing site $`x_0`$ is constructed with Swendsen-Wang bond-probabilities, and (2) all spins in this cluster are flipped, resulting in a new spin configuration. The current cluster is discarded after each iteration.

After a sufficient total number $`N(r)`$ of iterations, this process will “equilibrate” all spins within a radius $`r`$ around $`x_0`$, and two-point functions can then be measured within this area. With increasing number of iterations, the cluster will occasionally reach larger and larger distances. Because of eq. (4), the probability to do so is (roughly) proportional to the correlation function $`C(r)`$. A region of the lattice at distance $`r`$ from $`x_0`$ will become “equilibrated” after the cluster has reached that region a sufficient number $`n_{eq}`$ of times, thus $`N(r)\frac{n_{eq}}{C(r)}`$. Since the cluster can grow without bounds, we are calculating correlation functions for the infinite lattice, while we only need to store the spin configuration within the area actually reached by a cluster during a finite run.

Fig. 1 shows as an example the correlation function for the two dimensional Ising model at $`\beta =0.42<\beta _c=0.44068..`$, calculated according to eq. (4).

For comparison we show results for finite size lattices, with periodic boundary conditions. Usually, the infinite lattice result can only be obtained as the extrapolation of the finite lattice ones. In contrast, with our new method we can obtain the infinite latttice result in a single simulation, producing the lower straight line in fig. 1. We checked that the method works just as well in three dimensions.

The calculation of error bars for $`C(r)`$ needs special care. For each distance $`|r|`$, we demand that the cluster contains pairs of sites at that distance $`|r|^{d1}\times n_{eq}`$ times, i.e. it reaches each site at that distance $`O(n_{eq})`$ ($`=O(10)`$) times, before we consider data at that distance thermalized. Only then do we start to accumulate measurements for $`C(r)`$ according to eq. (4) for all sites $`i,j`$ in the cluster. Alternative similar procedures are possible.

Note that the overall computational time is by far dominated by measurements, both in usual cluster algorithms and in our method: In order to obtain a relative error $`ϵ`$ for $`C(r)`$, a large number of about $`1/ϵ^2`$ clusters with sites at distance $`r`$ need to occur, each of which is generated with a probability of about $`1/C(r)`$, in both cases. Thus our method is at least as efficient as usual cluster algorithms, while providing results for the infinite lattice directly.

The average computational time per cluster is proportional to the average cluster size, which is proportional to the susceptibility (see below eq. (4)). The effective simulated system size grows with the number of iterations.

When will the new method work ? The computational effort and thus the susceptibility needs to be finite. In addition, there must not be a finite probability for a percolating cluster, which implies $`\beta <\beta _c`$. The present implementation of the new method will thus work anywhere in the unbroken phase, right up to the critical point. Systems with long range order may become accessible with a suitable modification of the method. We note that indeed points away from the transition are often the region of interest, when comparing results of simulations to experimental measurements.

For nonrelativistic quantum systems, the loop algorithm  provides a cluster method. It is based on an enlargened representation in terms of the original spin operators and additional loop operators, similar in spirit to the Edwards-Sokal-Fortuin-Kasteleyn representation, eq. (1). We will in the following specialize to the spin $`\frac{1}{2}`$ quantum Heisenberg model

$$=J\underset{ij}{}\frac{1}{2}\left(S_i^+S_j^{}+S_i^{}S_j^+\right)+\lambda S_i^zS_j^z.$$
(5)
The two-point function for the single cluster version of the loop algorithm then reads

$$\begin{array}{cccc}\frac{1}{2}\left(S_i^+S_j^{}+S_i^{}S_j^+\right)=\delta (\text{i and j on the loop})\hfill & & & \\ S_i^zS_j^z=S_i^zS_j^z\delta (\text{i and j on the loop})\hfill & & & \end{array}.$$
(6)
Therefore our approach to simulate infinite size systems can be used in the same way as for the Ising model. “Infinite size” here can be applied to the spatial directions, as well as, independently, to the direction of imaginary time, yielding simulations directly at $`\beta =\mathrm{}`$, while retaining all dynamical information.

Note that our method is not a projector method. We can simulate directly at $`L=\mathrm{}`$ and $`\beta =\mathrm{}`$. This also enables us to obtain the limits $`q0`$ and $`\omega 0`$ directly from the simulations. Contributions from $`q0`$ or $`\omega 0`$, which are in general different from those limits, and are present in normal simulations with periodic boundary conditions, are completely avoided in our approach, as well as the double limit $`L\mathrm{}`$ and $`q0`$ that is usually required. When the number of iterations increases, clusters will reach larger distances in spatial and/or imaginary time direction, thereby improving the resolution in $`q`$ and $`\omega `$.

Again, our present method will be applicable provided that there is no percolation and that the unsubtracted correlation function drops off sufficiently quickly (i.e. that the corresponding susceptibility is finite).

As an example we studied spin ladder systems  with $`N=2`$ and $`4`$ legs for the isotropic antiferromagnet ($`\lambda =1`$). We used the discrete time version of the loop algorithm ($`\mathrm{\Delta }\tau =1/16`$). The method can be applied just as well in continuous time . Fig. 2 shows results for the equal time staggered spatial correlation functions along the chains, which behave similarly to the Ising case. A fit to the infinite lattice result gives $`\xi =2.93(2)`$ for $`N=2`$ and $`\xi =8.2(1)`$ for $`N=4`$. Fig. 3 shows greens functions for $`L=\mathrm{}`$, the infinite size system. Whereas finite temperature calculations give results periodic in imaginary time, which have to be extrapolated, the new approach provides the $`\beta =\mathrm{}`$ ($`T=0`$) result directly. With increased number of iterations, the infinite system results become available at larger distances both in spatial and in imaginary time direction. A fit to the exponential decay $`G(\tau )e^{\tau \mathrm{\Delta }}`$ of our data directly provides estimates for the gaps $`\mathrm{\Delta }=0.5059(4)`$ at $`N=2`$ and $`\mathrm{\Delta }=0.19(1)`$ at $`N=4`$, consistent with previous results. We also show results for $`L=20`$ and $`\beta =\mathrm{}`$ to exemplify the effect of finite size systems.

Continuing the imaginary time greens function to real frequencies by the maximum entropy method provides the spectra of fig. 4, in which the gaps, the single magnon peaks, and higher excitations for $`N=4`$ are clearly visible. This appears to be the first time that the spectrum for $`N=4`$ has been calculated.

Let us compare our method to other approaches. Exact diagonalization provides ground state results but is limited to small systems. Projector methods are also limited in system size and rely critically on proper convergence to the ground state. DMRG can achieve $`\beta =\mathrm{}`$ (mostly without dynamical information) for fairly small systems, or infinite size for large temperatures and finite $`\mathrm{\Delta }\tau `$. The most powerful method to extrapolate to infinite system size is the Finite Size Scaling method of Kim  and Caracciolo et al. , which allows extrapolation at correlation lengths far larger than the system size, but requires knowledge about scaling and corrections to scaling of the model.

In summary, we have introduced a new method, as a small modification of existing cluster methods, to simulate both classical and quantum systems at infinite system size and/or at zero temperature, while obtaining all two-point functions and derived quantities. Larger distances (in space and/or imaginary time) are accessed for longer simulations. The method in its present form is applicable to all cases where existing single-cluster algorithms can be used, as long as the unsubtracted two-point correlation function decays sufficiently quickly, including any $`\beta <\beta _c`$ in the unbroken phase. Such points away from a transition are often the region of interest when comparing simulation results to experimental measurements.