text:
# 1 Gauge Quantum Numbers and Yukawa Couplings of Scalar Diquarks (Q_{em}=T₃+{Y/2}).

MRI-P-000801

hep-ph/0008165

$`bs\gamma `$ confronts $`B`$-violating scalar couplings: $`R`$-parity violating supersymmetry or diquarks

Debrupa Chakraverty1 and Debajyoti Choudhury2

Mehta Research Institute, Chhatnag Road, Jhusi, Allahabad 211 019, India.
E-mail: 1rupa@mri.ernet.in, 2debchou@mri.ernet.in

Abstract

We investigate the possible role that baryon number violating Yukawa interactions may take in the inclusive decay $`BX_s\gamma `$. The constraints, derived using the experimental results of the CLEO collaboration, turn out, in many cases, to be more stringent than the existing bounds.

The Standard Model (SM) of the strong, weak and electromagnetic interactions is in very good agreement with almost all present experimental data, even though a few important predictions have not yet been tested. Still, most physicists would readily admit that the SM cannot be the final theory, both on aesthetic grounds as well as on account of certain well-founded technical objections. As a result, numerous attempts have been and are being made in the quest of a more fundamental theory. Experimentally, there have been two main strategies to probe new physics. On the one hand, we attempt to directly produce, and observe, new particles at high energy colliders. On the other, we look for virtual effects of such particles and/or interactions in various low and intermediate energy processes. The decay $`bs\gamma `$ is an excellent candidate for the latter option . Experimentally, the branching ratio for the inclusive decay $`BX_s\gamma `$ have been measured by CLEO  and ALEPH  to be

$$BR(BX_s\gamma )=\{\begin{array}{ccc}\hfill (3.15\pm 0.93)\times 10^4& & (\mathrm{CLEO})\hfill \\ \hfill (3.11\pm 1.52)\times 10^4& & (\mathrm{ALEPH}).\hfill \end{array}$$
(1)
The above are in good agreement with each other and with the SM prediction  of $`BR(BX_s\gamma )=(3.29\pm 0.33)\times 10^4`$. While a small window for the contribution of new physics does remain, this agreement can obviously be used to constrain deviations from the SM.

In this paper, we investigate the influence that a scalar diquark may have on the above decay111A brief discussion on the sensitivity of the branching ratio $`BX_s\gamma `$ to scalar diquark-top contribution has been presented in .. Diquarks abound in many grand unified theories (with or without supersymmetry) and even in composite models . While vector diquarks are constrained to be superheavy222We do not consider the case of non-gauged vector diquarks as such theories are nonrenormalizable. (with masses of the scale of breaking of the additional gauge symmetry), no such restrictions apply to the masses of scalar diquarks. Consequently, such particles can be as light as the electroweak scale. For example, a diquark like behaviour can be found even in a low energy theory like the Minimal Supersymmetric Standard Model (MSSM), albeit in the version with broken $`R`$-parity.

A generic diquark is a scalar or vector particle that couples to a quark current with a net baryon number $`B=\pm 2/3`$. Clearly, this may transform as either a $`SU(3)_c`$ triplet or sextet. Concentrating on the scalars (for reasons mentioned above), the generic Yukawa term in the Lagrangian can be expressed as

$$_Y^{(A)}=h_{ij}^{(A)}\overline{q}_i^cP_{L,R}q_j\mathrm{\Phi }_A+h.c.,$$
(2)
where $`i,j`$ denote quark flavours, $`A`$ denotes the diquark type and $`P_{L,R}`$ reflect the quark chirality. Standard Model gauge invariance demands that a scalar diquark transform either as a triplet or as a singlet under $`SU(2)_L`$ and that it have a $`U(1)`$ hypercharge $`|Y|=\frac{2}{3},\frac{4}{3},\frac{8}{3}`$. The full list of quantum numbers is presented in Table 1. It is clear that the couplings $`h_{ij}^{(1)}`$, $`h_{ij}^{(4)}`$, $`h_{ij}^{(5)}`$ and $`h_{ij}^{(7)}`$ must be symmetric under the exchange of $`i`$ and $`j`$ while $`h_{ij}^{(2)}`$, $`h_{ij}^{(3)}`$, $`h_{ij}^{(6)}`$ and $`h_{ij}^{(8)}`$ must be antisymmetric. For the other two, viz. $`\stackrel{~}{h}_{ij}^{(3)}`$ and $`\stackrel{~}{h}_{ij}^{(4)}`$, there is no particular symmetry property. Note that the quantum numbers of $`\mathrm{\Phi }_{2,4,6}`$ allow them to couple to a leptoquark (i.e. a quark-lepton) current as well. This implies that these particular diquarks could also mediate lepton-number ($`L`$) violating processes. Clearly, such leptoquark couplings need to be suppressed severely so as to prevent rapid proton decay.

We make a brief interlude here to discuss the MSSM . Whereas $`B`$ and $`L`$ are (accidentally) preserved in the SM (at least in the perturbative context), it is not so within the MSSM. Supersymmetry and gauge invariance, together with the field content, allow terms in the superpotential that violate either $`B`$ or $`L`$ . Catastrophic rates for proton decay can be avoided though by imposing a global $`Z_2`$ symmetry  under which the quark and lepton superfields change by a sign, while the Higgs superfields remain invariant. Representible as $`R(1)^{3BL+2S}`$, where $`S`$ is the spin of a field, this “$`R`$-parity” is positive for the SM fields and negative for all the supersymmetric partners. However, while this symmetry is useful in preventing phenomenologically unacceptable terms, it has no theoretical foundation and is entirely ad hoc in nature. Hence, it is of interest to examine the consequences of violating this symmetry, not in the least because it plays a crucial role in the search for supersymmetry. In our study, we shall restrict ourselves to the case where only the $`B`$-violating terms are non-zero. Such scenarios can be motivated from a class of supersymmetric GUTs as well . The corresponding terms in the superpotential can be parametrized as

$$W_{R/}=\lambda _{ijk}^{\prime \prime }\overline{U}_R^i\overline{D}_R^j\overline{D}_R^k$$
(3)
where $`U_R^i`$ and $`D_R^i`$ denote the right-handed up-quark and down-quark superfields respectively. The couplings $`\lambda _{ijk}^{\prime \prime }`$ are antisymmetric under the exchange of the last two indices. The corresponding Lagrangian can then be written in terms of the component fields as:

$$_{R/}=\lambda _{ijk}^{\prime \prime }\left(u_i^cd_j^c\stackrel{~}{d}_k^{}+u_i^c\stackrel{~}{d}_j^{}d_k^c+\stackrel{~}{u}_i^{}d_j^cd_k^c\right)+\mathrm{h}.\mathrm{c}.$$
(4)
Thus, a single term in the superpotential corresponds to three different diquark interactions, namely two of type $`\stackrel{~}{h}_{ij}^{(4)}`$ and one of type $`h_{ij}^{(8)}`$.

The best direct bound on diquark type couplings comes from the analysis of dijet events by the CDF collaboration . Considering the process $`q_iq_j\mathrm{\Phi }_Aq_iq_j`$, an exclusion curve in the ($`m_{\mathrm{\Phi }_A},h_{ij}^{(A)}`$) plane can be obtained from this data. Two points need to be noted though. At a $`p\overline{p}`$ collider like the Tevatron, the $`uu`$ and $`dd`$ fluxes tend to be small and hence the bounds are relatively weak. This is even more true for quarks of the second or third generation (which are relevant for the couplings that we are interested in). Secondly, such an analysis needs to make assumptions regarding the branching fraction of $`\mathrm{\Phi }_A`$ into quark pairs, a point that is of particular importance in the context of $`R`$-parity violating supersymmetric models.

There also exist some constraints derived from low energy processes. Third generation couplings, for example, can be constrained from the precision electroweak data at LEP  or, to an extent, by demanding perturbative unitarity to a high scale . Couplings involving the first two generations, on the other hand, are constrained333Although many of these analyses have been done for the case of $`R`$-parity violating models, clearly similar bounds would also apply to nonsupersymmetric diquark couplings as well. by the non-observance of neutron-antineutron oscillations or from an analysis of rare nucleon and meson decays . While many of these individual bounds are weak, certain of their products are much more severely constrained by the data on neutral meson mixing and $`CP`$–violation in the $`K`$–sector . It is our aim, in this article, to derive analogous but stronger bounds.

Within the SM, the quark level transition $`bs\gamma `$ is mediated, at the lowest order, by electromagnetic penguin diagrams shown in Fig. 1($`a`$$`d`$). While only the top-quark diagrams have been shown, for consistency’s sake, other charge $`2/3`$ quarks should also be included. However, these contributions are negligible on two counts: ($`i`$) the small mixing angles and ($`ii`$) the corresponding loop integrals being suppressed to a great extent due to the smallness of the light quark masses. The matrix element for this process is then governed by the dipole operator:

$$\frac{4G_F}{\sqrt{2}}V_{ts}^{}V_{tb}(\frac{e}{32\pi ^2})C_7^{SM}(m_W)\overline{s}\sigma _{\mu \nu }F^{\mu \nu }[m_b(1+\gamma _5)]b.$$
(5)
The QCD corrections to this process are calculated via an operator product expansion based on the effective Hamiltonian

$$H_{\mathrm{eff}}=\frac{4G_F}{\sqrt{2}}V_{ts}^{}V_{tb}\underset{i=1}{\overset{8}{}}C_i(\mu )O_i(\mu )$$
(6)
which is then evolved from the electroweak scale down to $`\mu =\mu _b`$ through renormalisation group (RG) equations. A large correction owes itself to the chromomagnetic operator $`bsG`$ ($`G`$ being a gluon)

$$\frac{4G_F}{\sqrt{2}}V_{ts}^{}V_{tb}(\frac{g_s}{32\pi ^2})C_8^{SM}(m_W)\overline{s}_\alpha \sigma _{\mu \nu }G_{\alpha \beta }^{\mu \nu }[m_b(1+\gamma _5)]b_\beta ,$$
(7)
which arises from the diagrams of Fig. 2(a-c)), The Wilson coefficients $`C_7^{SM}(m_W)`$ and $`C_8^{SM}(m_W)`$ can be evaluated perturbatively  at the $`W`$ scale where the matching conditions are imposed. The explicit expressions are

$$\begin{array}{ccc}\hfill C_7^{SM}(m_W)& =& x\left[\frac{75x8x^2}{24(x1)^3}+\frac{x(3x2)}{4(x1)^4}\mathrm{ln}x\right],\hfill \\ & & \\ \hfill C_8^{SM}(m_W)& =& x\left[\frac{2+5xx^2}{8(x1)^3}\frac{3x}{4(x1)^4}\mathrm{ln}x\right],\hfill \end{array}$$
(8)
where $`x=m_t^2/m_W^2`$.

The leading order results for the Wilson coefficients at $`\mu _b`$, the $`B`$-meson scale, is given by

$$C_7(\mu _b)=\eta ^{16/23}\left[C_7^{SM}(m_W)\frac{8}{3}C_8^{SM}(m_W)[1\eta ^{2/23}]+\frac{232}{513}[1\eta ^{19/23}]\right]$$
(9)
with $`\eta \alpha _s(m_W)/\alpha _s(\mu _b)`$, calculated using the leading $`\mu `$ dependence of $`\alpha _s`$, and the present world average value of the strong coupling constant viz. $`\alpha _s(m_Z)=0.118\pm 0.005`$. To this order, then,

$$\mathrm{\Gamma }(bs\gamma )=\frac{\alpha G_F^2m_b^5}{32\pi ^4}\left|V_{tb}V_{ts}^{}C_7^{SM}(\mu _b)\right|^2,$$
(10)
where $`\alpha `$ is the fine structure constant. As the above decay rate suffers from large uncertainties due to $`m_b`$ and the CKM matrix elements, it is prudent to normalise it against the measured semileptonic decay rate of the $`b`$ quark

$$\mathrm{\Gamma }(bce\overline{\nu }_e)=\frac{G_F^2m_b^5}{192\pi ^3}\kappa (z)g(z)|V_{cb}|^2,$$
(11)
where $`z=m_c^2/m_b^2`$ and

$$g(z)18z+8z^2z^412z^2\mathrm{ln}z$$

is the phase space factor. The analytic expression for $`\kappa (z)`$, the one loop QCD correction to the semileptonic decay, can be found in Ref.. The explicit dependence on $`m_b^5`$ is thus removed, while the ratio of the CKM elements in the scaled decay rate viz.,

$$\left|\frac{V_{ts}^{}V_{tb}}{V_{cb}}\right|=0.976\pm 0.010.$$
(12)
is much better known than the individual elements.

An updated next to leading order (NLO) analysis  of the $`BX_s\gamma `$ branching ratio with QED corrections has been presented in Ref.. Incorporating both the NLO QCD and the resummed QED corrections, the Wilson coefficient $`C_7^{\mathrm{eff}}(\mu _b)`$ in SM can be expanded as

$$C_7^{\mathrm{eff}}(\mu _b)=C_7(\mu _b)+\frac{\alpha _s(\mu _b)}{4\pi }C_7^{(1)}(\mu _b)+\frac{\alpha }{\alpha _s(\mu _b)}C_7^{(em)}(\mu _b).$$
(13)
For brevity’s sake, we do not give here the expressions for $`C_7^{(1)}(\mu _b)`$ and $`C_7^{(em)}(\mu _b)`$ as these can be found in Ref. and Ref. respectively. The inclusion of the NLO and QED corrections in the $`bs\gamma `$ decay rate has significantly reduced the large uncertainty present in the previous LO calculation. From the quark level $`bs\gamma `$ decay rate, it is possible to infer the $`B`$ meson inclusive branching ratio $`BR(BX_s\gamma )`$ by including the nonperturbative $`1/m_b`$ and $`1/m_c`$ corrections. These bound state corrections also have been taken into account in Ref..

Having delineated the formalism, it now remains to calculate the additional contributions due to the possible presence of non-zero diquark couplings. At the one-loop level, the only new contributions to $`bs\gamma `$ and $`bsG`$ arise from the diagrams of Fig. 1($`e`$$`h`$) and Fig. 2($`d`$$`g`$) respectively444Clearly, to this order, none of $`\mathrm{\Phi }_{5,6}`$ can mediate either of these processes and hence we shall not consider such fields any further.. In the generic case, apart from modifications in the coefficients of $`O_{7,8}`$, two additional operators arise. Denoted as

$$\begin{array}{ccc}\hfill \stackrel{~}{O}_7& =& \frac{e}{32\pi ^2}\overline{s}\sigma _{\mu \nu }F^{\mu \nu }[m_b(1\gamma _5)]b\hfill \\ \hfill \stackrel{~}{O}_8& =& \frac{g_s}{32\pi ^2}\overline{s}_\alpha \sigma _{\mu \nu }G_{\alpha \beta }^{\mu \nu }[m_b(1\gamma _5)]b_\beta ,\hfill \end{array}$$
(14)
they differ from their SM counterparts in their chirality structure555In all of these four operators, contributions proportional to the strange quark mass have been neglected..

To keep the analysis simple, we shall assume that only one diquark multiplet is light and that all the fields within a multiplet are degenerate666Large splittings within a multiplet is disfavoured by the electroweak precision data.. With this simplifying assumption, the new contributions, at the electroweak scale, are given by

$$\begin{array}{ccc}\hfill C_7^D(m_W)& =& \frac{N_c}{𝒜}\left[\left\{Q_\mathrm{\Phi }F_1(y)+Q_tF_3(y)\right\}l_bl_s^{}+\frac{m_t}{m_b}\left\{Q_tF_4(y)Q_\mathrm{\Phi }F_2(y)\right\}r_bl_s^{}\right],\hfill \\ \hfill \stackrel{~}{C}_7^D(m_W)& =& \frac{N_c}{𝒜}\left[\left\{Q_\mathrm{\Phi }F_1(y)+Q_tF_3(y)\right\}r_br_s^{}+\frac{m_t}{m_b}\left\{Q_tF_4(y)Q_\mathrm{\Phi }F_2(y)\right\}l_br_s^{}\right],\hfill \\ \hfill C_8^D(m_W)& =& 𝒜^1\left[\left\{C_\mathrm{\Phi }F_1(y)+C_tF_3(y)\right\}l_bl_s^{}+\frac{m_t}{m_b}\left\{C_tF_4(y)C_\mathrm{\Phi }F_2(y)\right\}r_bl_s^{}\right],\hfill \\ \hfill \stackrel{~}{C}_8^D(m_W)& =& 𝒜^1\left[\left\{C_\mathrm{\Phi }F_1(y)+C_tF_3(y)\right\}r_br_s^{}+\frac{m_t}{m_b}\left\{C_tF_4(y)C_\mathrm{\Phi }F_2(y)\right\}l_br_s^{}\right],\hfill \\ \hfill 𝒜& & 4\sqrt{2}G_FV_{ts}^{}V_{tb}m_\mathrm{\Phi }^2,\hfill \end{array}$$
(15)
where

$$\begin{array}{ccc}\hfill F_1(y)& =& \frac{1}{12(y1)^4}\left[6y^2\mathrm{ln}y2y^33y^2+6y1\right],\hfill \\ \hfill F_2(y)& =& \frac{1}{2(y1)^3}\left[1y^2+2y\mathrm{ln}y\right],\hfill \\ \hfill F_3(y)& =& \frac{1}{12(y1)^4}\left[2+3y6y^2+y^3+6y\mathrm{ln}y\right],\hfill \\ \hfill F_4(y)& =& \frac{1}{2(1y)^3}\left[34y+y^2+2\mathrm{ln}y\right],\hfill \end{array}$$
(16)
with $`y=m_t^2/m_\mathrm{\Phi }^2`$. The color factor $`N_c`$ is $`1`$ and $`2`$ for triplet and $`\overline{6}`$ scalar respectively. $`Q_t`$ and $`Q_\mathrm{\Phi }`$ are the charges of top quark and the diquark respectively. The color factors $`C_t`$ and $`C_\mathrm{\Phi }`$— for the diagrams of Fig.2$`f`$ and Fig.2$`g`$—are given in Table 2.

Note that, once again, we consider only such contributions, as involve the top quark. As is easy to ascertain from eq.(15), for other quarks in the loop, the corresponding integrals are too small to be of any consequence. Thus, any coupling to the diquarks $`\mathrm{\Phi }_{7,8}`$, for example, would not be constrained to an appreciable degree by radiative $`b`$ decays. In Table 2, we also display the relevant chiral Yukawa couplings (to the $`b`$\- and $`s`$-quarks) for different choices of the diquark.

In estimating the effects of scalar diqaurk couplings, it is useful to consider the ratios  $`\xi _{7,8}`$ with

$$\xi _71+\frac{C_7^D(m_W)}{C_7^{SM}(m_W)},$$
(17)
and similarly for $`\xi _8`$. For the new operators $`\stackrel{~}{O}_7`$ and $`\stackrel{~}{O}_8`$, we define,

$$\stackrel{~}{\xi }_7=\frac{\stackrel{~}{C}_7^D(m_W)}{C_7^{SM}(m_W)}$$
(18)
and $`\stackrel{~}{\xi }_8`$ in an analogous fashion. With these definitions, the $`BX_s\gamma `$ branching ratio can be written as,

$`BR(BX_s\gamma )`$ $`=`$ $`B_{22}(\delta )+(\xi _7^2+\stackrel{~}{\xi }_{7}^{}{}_{}{}^{2})B_{77}(\delta )+(\xi _8^2+\stackrel{~}{\xi }_{8}^{}{}_{}{}^{2})B_{88}(\delta )+\xi _7B_{27}(\delta )`$ (19)
$`+`$ $`\xi _8B_{28}(\delta )+(\xi _7\xi _8+\stackrel{~}{\xi }_7\stackrel{~}{\xi }_8)B_{78}(\delta ).`$
In a parton level analysis, the photon would be monochromatic, with $`E_\gamma =E_\gamma ^{max}=m_b/2`$. However, once the gluon Bremsstrahlung contribution is included, the photon spectrum becomes nontrivial and, for experimental purposes, one needs to make an explicit demand on the photon energy, namely

$$E_\gamma >(1\delta )E_\gamma ^{max},$$
(20)
where $`\delta `$ is the fraction of the spectrum above the cut. The values of $`B_{ij}(\delta )`$ are listed in ref. for different choices of the renormalisation scale $`\mu _b`$ and the cut off parameter on the photon energy $`\delta `$. As is well known, some ambiguities exist in the choice of $`\mu _b`$ which should, typically, lie in the region $`m_b/2`$ to $`2m_b`$. For our analysis, we used $`\mu _b=m_b`$ and $`\delta =0.9`$. We have checked that other values of $`\delta `$ do not change the bound significantly.

Since the new physics becomes operative only above the electroweak scale, the additional contributions to the operators $`O_7`$ and $`O_8`$ will only serve to change the Wilson coefficients at $`m_W`$. Of course, the additional operators $`\stackrel{~}{O}_7`$ and $`\stackrel{~}{O}_8`$ would influence the RG equations for $`C_7`$ and $`C_8`$ as well. However, since we are primarily interested in small $`C_{7,8}^D(m_W)`$, it is safe to neglect any term in the RG equations involving these coefficients.

In the absence of an $`L`$-violating coupling, these diquarks clearly do not influence the semileptonic decay modes of the $`B`$-meson. Thus, we may continue to normalize the radiative $`b`$-decay against $`bce\overline{\nu }_e`$ in order to avoid the severe dependence on $`m_b`$. In Fig. 3, we plot the branching ratio $`BR(BX_s\gamma )`$ in presence of a diquark (multiplet) of mass $`100\mathrm{GeV}`$. We continue to work under the assumption that only one pair of couplings are non-zero. Furthermore, we assume the said couplings to be real777The extension to complex couplings is straightforward. The imaginary parts, however, can be better constrained from an analysis of the $`CP`$ violating decay modes.. That the curves should be parabolic in the product of the two couplings in question is obvious. To appreciate the fact that many of these curves have their minima lying on the SM value, one needs to consider the chirality structure of the corresponding diquark couplings (see Table 2). For example, combinations involving either of $`\stackrel{~}{h}_{32}^{(3)}`$ or $`\stackrel{~}{h}_{32}^{(4)}`$, imply that there is no left-handed coupling to the strange quark. From eqs.(15), it is then easy to see that the SM Wilson coefficients $`C_{7,8}`$ remain unaffected. Consequently, the new contribution adds incoherently with the SM amplitude. For the rest of the combinations, though, the interference term is non-negligible leading to a shift in the minimum. Hence, unlike those for the first set of combinations (those involving $`\stackrel{~}{h}_{32}^{(3,4)}`$), the branching ratios corresponding to these sets are in agreement with the experimental numbers for two non-contiguous ranges of the product.

In each of eqs.(15), the second term, whenever allowed, is clearly the dominant piece. This enhancement by the factor of $`m_t/m_b`$ comes into play only when the diquark couplings to the bottom and the strange quarks have pieces with opposite chirality. In other words, for a diquark of the type $`\mathrm{\Phi }_3`$ (or $`\mathrm{\Phi }_4`$), the simultaneous presence of both the allowed types of couplings is severely constrained by the data on $`BX_s\gamma `$.

In Table 3, we capture the essence of Fig. 3 in the form of actual limits that can be set on such products of couplings, for a diquark mass of $`100\mathrm{GeV}`$. Understandably, the $`2\sigma `$ bounds are weaker than the $`1\sigma `$ ones. Similarly, the color-sextet couplings are more severely constrained than the color-triplet ones. It should be noted that the structure of the interaction terms $`h_{ij}^{(1)}`$ are the same as those for $`h_{ij}^{(3)}`$. Consequently, the bounds are exactly the same. A similar story obtains for $`h_{ij}^{(2)}`$ and $`h_{ij}^{(4)}`$. As discussed above, for a few of the products there are two non-contiguous bands allowed. For the combinations $`h_{33}^{(2)}h_{32}^{(2)}`$ and $`h_{33}^{(4)}h_{32}^{(4)}`$, though, the second window (both at the $`1\sigma `$ and $`2\sigma `$ levels) lies beyond the perturbative limit and, hence, are phenomenologically uninteresting.

As discussed earlier, $`\stackrel{~}{h}_{ij}^{(4)}`$ is analogous to the trilinear $`R`$-parity violating coupling $`\lambda _{ijk}^{\prime \prime }`$. Thus the constraints on $`\stackrel{~}{h}_{33}^{(4)}\stackrel{~}{h}_{32}^{(4)}`$ are equivalent to those on the product $`\lambda _{3j2}^{\prime \prime }\lambda _{3j3}^{\prime \prime }`$. For each of these couplings, the best individual bound comes from the precision measurements at the $`Z`$ pole , and amounts to $`\lambda _{3j2}^{\prime \prime },\lambda _{3j3}^{\prime \prime }<0.50`$ at the $`1\sigma `$ level. We thus do not do very well as far as this particular combination is concerned. This can be attributed to both the chirality and the color structure of the operator, each of which is “unfavourable” as far the $`bs\gamma `$ decay is concerned. For most of the other combinations though, we do significantly better than the product of individual bounds .

In our effort to compare with the results available in the literature, we have, until now, held the diquark mass to be $`100\mathrm{GeV}`$ and varied the strength of its coupling. In reality, though a diquark is more likely to be somewhat heavier. For the sake of completeness, we next investigate the dependence on the diquark mass (Figs. 4), while holding the product fixed. As is expected, the extra contribution falls off with $`m_\mathrm{\Phi }`$. The fall-off is somewhat slower than $`m_\mathrm{\Phi }^2`$ (see the expressions for $`F_i(y)`$ in eq.(16)) and the effects persist till $`3TeV`$. The different rates of fall-off are governed by the dominant $`F_i(y)`$ in each case. The case for the combination $`\stackrel{~}{h}_{33}^{(3)}h_{32}^{(3)}`$ looks somewhat nontrivial. However, the shape is just a consequence of accidental cancellations between various terms of eq.(19). As the exact nature of these cancellations depend crucially on the value of the diquark couplings, not much should be read into the shape in general or the minimum in particular.

The two dependences ($`m_\mathrm{\Phi }`$ and coupling strength) that we have studied can be combined to rule out parts of the phase space. In Figs.5, we exhibit this for two particular combinations. In each, the shaded regions of the parameter space are in agreement with the experimental results at the designated level. For $`h_{33}^{(1)}h_{32}^{(1)}`$, the second allowed region is beyond the perturbative limit.

In summary, we have studied the effects of the scalar diquark and/or R prity violating coupling to the branching ratio $`BX_s\gamma `$. Among the possible new contributions, the scalar diquark mediated diagram yield promising effects. The precise measurement of this branching ratio at the upcoming B factories in near future and the reduction of theoretical uncertainty will improve the limits on the product combination of different scalar diquark and/or R parity violating couplings, we obtained.

Acknowledgement

D. Choudhury acknowledges the Department of Science and Technology, India for the Swarnajayanti Fellowship grant.


format: text

text:
# HST Spectroscopy of Spot 1 on the Circumstellar Ring of SN 1987A

## 1 INTRODUCTION

After the discovery of the circumstellar ring around Supernova (SN) 1987A, several authors pointed out that the supernova blast wave should strike the ring some 10 – 20 years after the supernova explosion (Luo & McCray, 1991; Luo, McCray, & Slavin, 1994; Chevalier & Dwarkadas, 1995; Borkowski, Blondin, & McCray, 1997b). The estimated time of first contact was uncertain, mainly because it depended on the unknown density of circumstellar gas between the supernova and the ring.

This event is now underway. In 1995 an unresolved spot at position angle $`29^{}`$ (Spot 1) on the near (North) side of the ring began to brighten (Garnavich, Kirshner, & Challis, 1997). This spot continues to brighten, and within the last year several new spots have appeared (Lawrence & Crotts, 2000; Garnavich, Kirshner, & Challis, 2000; Lawrence et al., 2000). A shock interpretation for these spots was suggested by Space Telescope Imaging Spectrograph (STIS) observations of the ring taken in 1997 (Sonneborn et al., 1998). The spectral image of the ring in H$`\alpha `$ showed a streak at the position of Spot 1 indicating blueshifted emission with velocities up to $`250\text{ km s}^1`$. Evidently, at the positions of the spots, the blast wave is striking inward protrusions of the dense equatorial ring.

Here (§2) we report on early optical and ultraviolet STIS spectra of SNR 1987A which show several broad lines associated with Spot 1. In §3 we present the radiative shock model for the emission from the spots, and in §4 we discuss what we can infer from the spectrum and its evolution.

## 2 OBSERVATIONS

The Supernova INtensive Study (SINS) collaboration obtained STIS spectra of Spot 1 at ultraviolet and optical wavelengths in 1997 September and 1998 March, respectively. The ultraviolet spectrum (not shown) was taken on September 27, 1997 (3869.3 days since explosion) with the G140L low dispersion grating (1150 – 1736 Å, $`\mathrm{\Delta }v300\text{ km s}^1`$) and a $`0\stackrel{}{\mathrm{.}}5`$ slit. Five spectral images totaling $`11,200\text{ s}`$ were combined. Emission lines from nearly stationary gas in the ring are apparent at N v $`\lambda \lambda 1239,1243`$, O iv\] $`\lambda 1400`$, N iv\] $`\lambda 1487`$, and He ii $`\lambda 1640`$. Emission from Spot 1 is also visible in these lines, as well as at Si iv $`\lambda \lambda 1394,1403`$, \[N iv\] $`\lambda 1483`$, and C iv $`\lambda \lambda 1548,1551`$. The low-velocity Si iv and C iv emission from the ring is blocked by interstellar absorption (Fransson et al., 1989). Due to geocoronal Ly$`\alpha `$ emission no information about Ly$`\alpha `$ from the spot was obtained.

The optical spectrum was taken on March 7, 1998 (4030.0 days since explosion), with a $`0\stackrel{}{\mathrm{.}}2`$ slit and the G750M grating (6295 – 6867 Å, $`\mathrm{\Delta }v50\text{ km s}^1`$). Three spectral images totaling $`8,056\text{ s}`$ were combined. A section of the optical spectrum is shown in Fig. 2. The central horizontal streak (which extends to $`\pm 3,000\text{ km s}^1`$) is H$`\alpha `$ emission from the radioactively heated supernova debris. Pairs of bright spots where the slit intersects the stationary inner ring are evident in \[N ii\] $`\lambda \lambda 6548,6584`$ and H$`\alpha `$ (not displayed but also present are \[O i\] $`\lambda \lambda 6300,6364`$, He i $`\lambda 6678`$, and \[S ii\] $`\lambda \lambda 6717,6731`$). Three fainter spots due to the outer rings are also visible in H$`\alpha `$ and \[N ii\]. Broad streaks of emission from Spot 1 appear on the upper segments of the ring lines. These streaks are inside the ring emission indicating that Spot 1 lies at the inner edge of the ring.

The spectra of Spot 1 are extracted from the STIS spectral images. For the optical emission lines, the major source of contamination is the emission from the nearly stationary ($`v=10.5\pm 0.3\text{ km s}^1`$, Crotts & Heathcote 1999) circumstellar ring filling the $`0\stackrel{}{\mathrm{.}}2(100\text{ km s}^1)`$ slit. The H$`\alpha `$ line profile is shown in Figure 2. The parts of the line profiles that lie outside of the ring emission are fit well by Gaussians. The velocity shifts of the Gaussian fits are poorly constrained since the spot’s spatial position in the slit is uncertain. In the ultraviolet the ring emission is much weaker and therefore its contribution to the extracted Spot 1 spectrum is neglected. The UV lines are poorly resolved at the resolution of the G140L grating. It is possible to account for the detector line spread function in order to obtain intrinsic widths for the UV lines (Pun et al., 2000). While this method introduces errors in the determined width, we find that all the UV lines have intrinsic widths $`<400\text{ km s}^1`$.

All the lines from Spot 1 were fit with Gaussians in order to extract their fluxes and full widths at half maximum (FWHM) (see Table 1). The errors stated in Table 1 are $`1\sigma `$ statistical errors to the fits. Close multiplets (N v $`\lambda \lambda 1239,1243`$ and C iv $`\lambda \lambda 1548,1551`$) were fit simultaneously by requiring a set wavelength separation, identical widths, and line ratios dictated by their oscillator strengths. For the N v doublet, the emission from Spot 1 overlapped with very broad Ly$`\alpha `$ emission ($`\pm 15,000\text{ km s}^1`$) from the reverse shock (Sonneborn et al., 1998; Michael et al., 1998). This Ly$`\alpha `$ background was removed using a quadratic fit and is the source of the large uncertainty in the N v $`\lambda \lambda 1239,1243`$ flux. At the G140L resolution, the Si iv $`\lambda 1403`$ emission of Spot 1 is mixed with that from the O iv\] $`\lambda 1400`$ multiplet. Using the measured Si iv $`\lambda `$1394 flux and assuming $`I(1394)/I(1403)=2`$, the Si iv $`\lambda 1403`$ flux was determined and removed from the O iv\] $`\lambda 1400`$ multiplet.

Emission Lines from Spot 1

| Emission Line | FWHMa | Observed Fluxb |
| --- | --- | --- |
| N v $`\lambda \lambda 1239,1243`$ | $`\mathrm{}`$ | $`2.5\pm 0.9`$ |
| Si iv $`\lambda \lambda 1394,1403`$ | $`\mathrm{}`$ | $`0.6\pm 0.2`$ |
| O iv\] $`\lambda 1400`$ | $`\mathrm{}`$ | $`1.3\pm 0.4`$ |
| $`[`$N iv\] $`\lambda 1483`$ | $`\mathrm{}`$ | $`1.4\pm 0.4`$ |
| N iv\] $`\lambda 1487`$ | $`\mathrm{}`$ | $`1.5\pm 0.4`$ |
| C iv $`\lambda \lambda 1548,1551`$ | $`\mathrm{}`$ | $`2.9\pm 0.8`$ |
| He ii $`\lambda 1640`$ | $`\mathrm{}`$ | $`3.5\pm 0.9`$ |
| $`[`$O i\] $`\lambda 6300`$ | $`154\pm 7`$ | $`8.9\pm 0.7`$ |
| $`[`$O i\] $`\lambda 6364`$ | $`179\pm 14`$ | $`3.3\pm 0.4`$ |
| $`[`$N ii\] $`\lambda 6548`$ | $`165\pm 11`$ | $`5.6\pm 0.8`$ |
| H$`\alpha `$ $`\lambda 6563`$ | $`253\pm 4`$ | $`46.7\pm 1.2`$ |
| $`[`$N ii\] $`\lambda 6584`$ | $`204\pm 7`$ | $`13.1\pm 0.9`$ |
| He i $`\lambda 6678`$ | $`177\pm 38`$ | $`0.8\pm 0.2`$ |
| $`[`$S ii\] $`\lambda 6717`$ | $`113\pm 40`$ | $`0.9\pm 0.9`$ |
| $`[`$S ii\] $`\lambda 6731`$ | $`149\pm 48`$ | $`0.9\pm 0.4`$ |
| a in km s-1 | | |
| --- | --- | --- |
| b in units of $`10^{16}\text{ ergs cm}^2\text{ s}^1`$ | | |

In addition to interstellar extinction, the observed C iv and Si iv fluxes are further reduced by line absorption along the line of sight to the supernova. The absorption profiles for these lines are known from IUE high-dispersion spectra of SN 1987A taken at the time of outburst (Welty et al., 1999). Since the intrinsic profiles of the lines are poorly known, correcting for this absorption is a somewhat uncertain endeavor. For now we just present the observed fluxes and note that under reasonable assumptions of the intrinsic profiles we find correction factors of $`2`$ for these lines (Pun et al., 2000).

## 3 RADIATIVE SHOCK MODEL

Figure 3 illustrates the scenario that we believe accounts for the emission from Spot 1. The freely expanding supernova ejecta are slowed from $`15,000\text{ km s}^1`$ to $`4,000\text{ km s}^1`$ by a reverse shock. The shocked gas drives a blast wave (forward shock), with velocity $`v_\mathrm{b}4,000\text{ km s}^1`$, into an H ii region having density $`n_{\mathrm{II}}100\text{ cm}^3`$ (Chevalier & Dwarkadas, 1995; Borkowski, Blondin, & McCray, 1997a; Lundqvist, 1999). The bright spots occur where the blast wave encounters protrusions of the dense ($`n10^4\text{ cm}^3`$, Lundqvist & Fransson 1996) equatorial ring. As it enters a protrusion, the blast wave is slowed to a value $`v_\mathrm{s}v_\mathrm{b}(n_{\mathrm{II}}/n)^{1/2}f(\theta )`$, where $`\theta `$ is the angle between the surface of the protrusion and the propagation direction of the blast wave. The function $`f(\theta )`$, which accounts for shock obliquity and the pressure increase due to the reflected shock, varies from $`f(0)2`$ at the head of the protrusion to $`f(\pi /2)0.7`$ along its sides (Borkowski, Blondin, & McCray, 1997b). We therefore expect a range of shock velocities, $`280n_4^{1/2}\text{ km s}^1<v_\mathrm{s}<800n_4^{1/2}\text{ km s}^1`$, to be present (where $`n_4`$ is the pre-shock density in units of $`10^4\text{ cm}^3`$).

An important parameter for the shocks considered here is the radiative cooling time $`t_\mathrm{c}`$, defined as the time for the shocked gas to cool from the post-shock temperature, $`T_\mathrm{s}=3\mu m_\mathrm{p}v_\mathrm{s}^2/16k1.2\times 10^6(v_\mathrm{s}/300\text{ km s}^1)^2\text{ K}`$, to $`T_\mathrm{c}=10^4\text{ K}`$. A fit to the results from a plane-parallel steady-state (1D-SS) shock code (e.g., Cox & Raymond 1985) with abundances appropriate to SN 1987A (Lundqvist & Fransson, 1996; Russell & Dopita, 1992) shows that, for $`100\text{ km s}^1<v_\mathrm{s}<600\text{ km s}^1`$, $`t_\mathrm{c}4.6n_4^1(v_\mathrm{s}/300\text{ km s}^1)^{3.7}\text{ yr}`$.

When these spectra were taken, the age of the spot, $`t_{\mathrm{spot}}3\text{ yr}`$, was comparable to the cooling time for shocks with the range of densities and velocities expected in the protrusion. Shocks with velocities $`>250n_4^{0.27}\text{ km s}^1`$ have $`t_\mathrm{c}>t_{\mathrm{spot}}`$ and would not yet have radiated away their energy, while slower (radiative) shocks will have had time to cool and will convert all their thermal energy into radiation. Non-radiative shocks produce far less ultraviolet and optical emission than radiative shocks. Therefore, although faster shocks might be present, we expect the ultraviolet and optical emission from Spot 1 to be dominated by shocks with speeds $`<250n_4^{0.27}\text{ km s}^1`$.

The radiative shock model does predict the presence of both strong UV and optical lines. The ultraviolet lines produced in a radiative shock come mainly from the cooling region, where the shocked gas’s thermal energy is radiated away. The optical lines, however, are formed in the cool layer of dense ($`nT_\mathrm{s}/T_\mathrm{c}`$, where $`T_c10^4\text{ K}`$) gas which forms behind the cooling region. This photoabsorption layer quickly becomes thick enough to absorb the downstream ionizing radiation produced in the cooling region and converts the flux into optical emission lines.

The optical lines are formed with thermal widths characteristic of $`10^4\text{ K}`$ gas. The actual profiles are much broader and must be due to the macroscopic motion of fluid parcels in the photoabsorption layer. In a steady shock, the gas in the photoionization zone will have the same velocity as the shock front, and therefore the optical lines would be expected to have velocity profiles representative of the hydrodynamics of the shock (specifically, those parts of the shock that have developed a radiative layer). The presence of both redshifted and blueshifted emission can be explained by shocks traveling both into and out of our line of sight (see Figure 3).

## 4 DISCUSSION

It will not be an easy task to develop a quantitative model for the spectrum and line profiles of the hot spot. As Figure 3 illustrates, these result from a superposition of shocks having a range of velocities and aspect angles that depend on the unknown geometry and density structure of the protrusion. Our task is further complicated by the fact that radiative shocks are subject to violent thermal and dynamical instabilities (e.g., Innes, Giddings, & Falle 1987; Klein, McKee, & Colella 1994; Walder & Folini 1996, 1998), so steady shock models may be inadequate to interpret the spectrum.

Despite the complexity of the actual situation, we can gain some insight into the physical conditions in the spot by comparing the actual spectra with the results from the 1D-SS shock code. For example, the code shows that the ratios of the flux of N v $`\lambda \lambda 1239,1243`$ to the flux of lines from lower ionization states (O iv\] $`\lambda 1400`$, N iv\] $`\lambda 1487`$, and C iv $`\lambda \lambda 1548,1551`$) are sensitive functions of shock velocity. N v is not produced at temperatures $`T_\mathrm{s}<2\times 10^5\text{ K}`$ ($`v_\mathrm{s}<120\text{ km s}^1`$), above which its emissivity rises rapidly with shock velocity. On the other hand, O iv, N iv, and C iv all form at lower temperatures and their emissivities are less sensitive to velocity. To reproduce the observed ratios with the shock code we must constrain the shock velocity to the range $`120\text{ km s}^1<v_\mathrm{s}<150\text{ km s}^1`$. But such a velocity range conflicts with the observed optical line widths, which require radiative shocks with velocities up to at least $`250\text{ km s}^1`$.

However, only $`20\%`$ of the emission in the broad H$`\alpha `$ line profile comes at high velocity ($`|v|>135\text{ km s}^1`$, measured from line center). We may then argue that only $`20\%`$ of the H$`\alpha `$ emission must come from shocks faster than $`135\text{ km s}^1`$. If this is true, the line ratios will be dominated by emission from the slower shocks. We find that if the slower shocks ($`v_s<135\text{ km s}^1`$) cover about four times the surface area than the faster shocks ($`v_s>135\text{ km s}^1`$), then the net emission from these shocks produced ratios near those observed.

The presence of both fast and slow radiative shocks indicates that the protrusion must have a slightly higher density than the average ring density. Using hydrodynamical arguments presented in §3 we find that preshock densities $`n>4\times 10^4\text{ cm}^3`$ are required to produce shocks as slow as $`135\text{ km s}^1`$. A less stringent lower limit is provided by an independent argument based on the radiative cooling time. Preshock densities of $`n>10^4\text{ cm}^3`$ are required for shocks as fast as $`250\text{ km s}^1`$ (as indicated by the H$`\alpha `$ line width) to cool within $`t_s3`$ years.

We can also use the spectroscopic data to estimate the surface area of the shock interaction. Most of the emission comes from shocks with velocities $`135\text{ km s}^1`$, for which the shock code gives an H$`\alpha `$ surface emissivity of $`0.4n_4\text{ ergs cm}^2\text{ s}^1`$. Correcting the observed H$`\alpha `$ flux for reddening (Scuderi et al., 1996) and assuming a distance of $`50\text{ kpc}`$, we find that the shock must have an H$`\alpha `$ emitting surface area of $`A_s10^{33}\text{ cm}^2`$ (assuming $`n=4\times 10^4\text{ cm}^3`$). Then, assuming a hemispherical shape for the emitting surface, we estimate that Spot 1 should have a characteristic dimension $`D_{\mathrm{spot}}(2A_s/\pi )^{1/2}3\times 10^{16}\text{ cm}`$ (or $`0\stackrel{}{\mathrm{.}}04`$). This is smaller than HST’s resolution so we are not surprised that the spot is spatially unresolved.

As Table 1 indicates, different emission lines have different widths. This is expected, since the emissivities of the lines are sensitive functions of shock velocities, and a range of shock velocities is clearly in evidence. Moreover, the differences among the optical line widths may be explained by partial cooling in fast shocks. Emission from ions like O i and S ii comes from regions downstream of the gas that produces H$`\alpha `$ emission. Faster shocks take longer to develop these downstream regions. The fastest shocks that produce H$`\alpha `$ will not yet have developed their \[O i\] and \[S ii\] emitting zones. Therefore, the \[O i\] and \[S ii\] lines are narrower than the H$`\alpha `$ line.

Now that several new bright spots have appeared on the circumstellar ring of SN 1987A (Lawrence & Crotts, 2000; Garnavich, Kirshner, & Challis, 2000; Lawrence et al., 2000), we have new opportunities to elucidate the development of radiative shocks. By observing the spectra of each of the spots, we can follow the development of shocks having different ages and velocities. The spots will continue to brighten and their spectra will change. The pressure that drives the shocks will increase, the shocks will engulf more of the protrusions, and faster shocks will develop radiative layers. The most powerful constraint on the hydrodynamics will be to resolve the line ratios as a function of Doppler velocity. For example, the wings of the N v $`\lambda \lambda 1239,1243`$ line profile should be more prominent than those of C iv $`\lambda \lambda 1548,1551`$. As the optical lines brighten, their wings should broaden.

For these reasons, it is imperative to observe the development of the spectra of these spots on a regular basis. These observations will need to be performed by STIS in order to obtain UV line profiles as well as spatially distinguish the spectra of the different spots.

E.M. thanks S. Zhekov for fruitful discussions. This research was supported by NASA through grants NAG 5-3313 and NTG 5-80 to the University of Colorado and grants GO-2563 and GO-7434 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. C.S.J.P. acknowledges support by the STIS IDT through NOAO by NASA.


text:
# Topology and metastability in the lattice Skyrme model

## I Introduction

The Skyrme model  is a popular model of the dynamics of pions and nucleons, incorporating the former as its fundamental, pseudo-Goldstone field and the latter as topological solitons. The continuum theory has been widely studied via semiclassical techniques, giving a satisfactory phenomenology of low-momentum and low-temperature physics. In this paper we develop a lattice formulation of the model.***For previous work see .

Our motivation in employing a lattice cutoff is to overcome the limitations of the continuum theory. In point of fact, the term “continuum theory” is misleading. The Skyrme Lagrangian is a non-linear sigma model with a four-derivative term, which makes it non-renormalizable in perturbation theory. This means that calculations of quantum effects must involve a short-distance cutoff. In working with chiral Lagrangians coupled to a local nucleon field, this is not a serious problem . One absorbs divergences into an ever-lengthening list of counterterms; as long as external momenta are kept small, dimensional analysis limits the contributions of these counterterms to the results. When the nucleon is a soliton rather than a fundamental field, however, things are more difficult. The soliton’s collective degrees of freedom are quantized separately, and the systematic development of higher-order quantum effects involves disentangling the pion field from these collective quanta . In any case, there is no way to limit the addition of higher-dimension terms to the Lagrangian as long as one must consider energy scales approaching the nucleon mass.

We propose to turn these points to advantage, by considering the Skyrme model as an effective field theory. The cutoff, and the cutoff scheme (the lattice), are part and parcel of the specification of the theory. The free choice of terms in the Lagrangian is now a virtue, constrained only by phenomenological necessity. The procedure will be to pick a Lagrangian; to fix the lattice spacing; and finally to calculate any physical quantity of interest, going beyond perturbation theory and beyond semiclassical methods, since the lattice offers many more direct techniques. From a lattice point of view, the fixed cutoff is a virtue as well, since no continuum limit is necessary (nor is it possible).

The most interesting feature of the Skyrme model is the stability conferred on the soliton by the topological conservation law and the four-derivative coupling. We wish in fact to isolate a single skyrmion in order to see how its properties change with temperature, density, etc. A recurring problem in lattice studies of solitons is the tendency of lattice dislocations to destroy topological stability or, in other words, to allow topology to slip through the lattice. We avoid this problem by choosing a lattice action that enforces continuity at short distances. In order to keep our numerical evolution within a single topological sector, we insist on a local updating scheme that should not nucleate smooth skyrmions. We find that these two ingredients suffice to stabilize the lattice skyrmion.

The main technical development in this paper is the construction of an exact topological density that can be measured on a lattice field configuration. By “exact” we mean that the winding number, the sum of the density over the lattice, is always an integer. We cut the lattice into fundamental tetrahedra and map each tetrahedron into a curved tetrahedron in the $`S^3`$ target space. The winding number is then the sum of the signed volumes of the tetrahedra in $`S^3`$. We calculate the volume of a spherical tetrahedron via a formula due to Coxeter , derived as a solution of differential equations first written down by Schläfli .

We summarize the continuum Skyrme model in Sec. II in order to establish notation. General considerations regarding continuity and topology on the lattice are presented in Sec. III, which concludes with formulae for the lattice action we employ. In Sec. IV we present our definition of the local topological density via Coxeter’s formula for the volume of a quadrirectangular tetrahedron in $`S^3`$. As a first application, we present in Sec. V the results of Monte Carlo simulation for the classical Skyrme model in three dimensions. Easier to simulate than the full four-dimensional path integral, this model may be regarded as a Ginzburg-Landau theory for the Skyrme model at finite temperature—a sort of effective theory for the effective theory. We use the topological density to calculate the mean-square radius of a skyrmion as a function of the couplings $`\beta _1`$ and $`\beta _2`$ of the lattice action. If we fix $`\beta _2/\beta _1`$, so that the form of the action is fixed, we find a multitude of metastable configurations of the skyrmion that are smaller than the ground state solution. As the “temperature” is raised by decreasing $`\beta _1`$ and $`\beta _2`$, these metastable states are made accessible to fluctuations, with the result that the skyrmion shrinks as it is “heated.” Presumably the “temperature” of the 3d theory is an increasing function of the real temperature in the 4d theory that it approximates. Thus we reach the result that a quantized skyrmion shrinks as it is heated.

## II Continuum Skyrme model

The Euclidean action of the Skyrme model is

$$S=d^4x\left[\frac{f_\pi ^2}{16}\text{Tr}|_\mu U|^2+\frac{1}{32e^2}\text{Tr}([L_\mu ,L_\nu ])^2\right],$$
(1)
where the non-linear chiral field $`U`$ is an $`SU(2)`$ matrix, and we have defined

$$L_\mu =iU^{}_\mu U.$$
(2)
We have omitted a mass term of the form $`m_\pi ^2\text{Tr}U`$ and thus $`S`$ is invariant under the $`SU(2)\times SU(2)`$ group of chiral rotations,

$$UAUB^{},\text{ with }A,BSU(2).$$
(3)
Classically, the symmetry is spontaneously broken. Finite-energy field configurations must tend to a constant at infinity, and one can use a symmetry rotation to make this constant the unit matrix, i.e.,

$$U(𝐱)𝟙\text{ as }|𝐱|\mathrm{}.$$
(4)
Using the Pauli matrices $`\tau _i`$, $`i=1,2,3`$, we can write $`U`$ in terms of new fields $`𝝈=(\sigma _0,\sigma _i)`$ via

$$U=\sigma _0+i\sigma _i\tau _i,$$
(5)
where $`𝝈𝝈=1`$. By considering small fluctuations about $`U=𝟙`$, we can identify $`\sigma _i`$, $`i=1,2,3`$ with the Goldstone pion field.

Equation (4) means that 3-dimensional space is compactified to the 3-sphere $`S^3`$. Since $`USU(2)`$ also takes values in $`S^3`$ \[see Eq. (5)\], configurations $`U(𝐱)`$ can be classified according to the homotopy group $`\pi _3(S^3)=Z`$. Thus there is an integer winding number $`n`$ that denotes how many times $`U(𝐱)`$ covers the 3-sphere in field space as x is varied over its 3-sphere. This winding number is topologically conserved, meaning that it cannot change under continuous deformation of the field $`U(𝐱)`$.

The $`n=0`$ sector includes the vacuum $`U(𝐱)=𝟙`$ and perturbations around it. An example of an $`n=1`$ configuration is the spherically symmetric skyrmion,

$$U(𝐱)=\mathrm{exp}\left[if(r)\frac{x_i\tau _i}{r}\right],$$
(6)
with

$`f(0)`$ $`=`$ $`\pi ,`$ (7)
$`f(\mathrm{})`$ $`=`$ $`0.`$ (8)
$`f(r)`$ should be determined so as to minimize the static energy

$$E=𝑑𝐱\left[\frac{f_\pi ^2}{16}\text{Tr}|_iU|^2+\frac{1}{32e^2}\text{Tr}([L_i,L_j])^2\right],$$
(9)
but a simple choice with the right topology is

$$f(r)=\pi \left(1\mathrm{tanh}\frac{r}{r_0}\right).$$
(10)
Skyrme identified $`n`$ with the baryon number of a field configuration, and the lowest soliton configuration with the nucleon.

Given an arbitrary field configuration $`U(𝐱)`$ satisfying Eq. (4), its winding number may be calculated with the formula

$$n=\frac{i}{24\pi ^2}𝑑𝐱ϵ_{ijk}\text{Tr}L_iL_jL_k.$$
(11)
The geometric meaning of Eq. (11) will be apparent in its lattice counterpart below.

## III Lattice topology and continuity

In this section we define a topological density for a lattice field configuration that is unambiguous and that sums exactly to an integer. The further demand of conservation of the winding number will lead us to choosing a lattice action that constrains discontinuities in the field.

A lattice configuration is specified by the field $`U_𝐧SU(2)`$ or, equivalently, by the 4-vector $`𝝈_𝐧S^3`$. In order to define the winding number, we begin  by cutting the cubic lattice $`L^3`$ into tetrahedra, five tetrahedra per cubic cell (see Fig. 1).

The four vertices $`𝐧^{(i)}`$, $`i=1`$–4 of each tetrahedron map to four unit 4-vectors $`𝝈^{(i)}`$, which are vertices of a spherical tetrahedron in $`S^3`$. This tetrahedron is defined via its vertices; its edges are arcs of great circles and its faces are spherical triangles drawn on great spheres.

Since three (non-collinear) points in $`S^3`$ determine a great sphere, it is clear that two adjacent tetrahedra in $`L^3`$ map into adjacent tetrahedra in $`S^3`$, with the common face in $`L^3`$ mapping into a common face in $`S^3`$. Thus the field configuration $`𝝈_𝐧`$ gives a triangulation of some volume in $`S^3`$. We impose periodic boundary conditions on the lattice, and thus no tetrahedron possesses a face that is not shared by another tetrahedron (a free face). This implies that the complex of tetrahedra in $`S^3`$ possesses no free faces either, and thus the volume covered by the complex must be an integer multiple of the volume The three-dimensional volume of the unit 3-sphere is $`2\pi ^2`$. We will redefine this to be unity, and thus measure volumes in units of $`2\pi ^2`$. of $`S^3`$. If the field configuration is smooth, the complex in $`S^3`$ will be composed of tetrahedra that are small compared to $`S^3`$, and it will wrap around the sphere much as a smooth mapping of the 3-torus must. In any case, the volume in $`S^3`$ of the five tetrahedra corresponding to a cube in $`L^3`$ gives a definition of the topological density $`\rho _𝐧`$ contained in the cube.

The topological density thus defined is not unambiguous, however. If a tetrahedron in $`S^3`$ is specified by its faces, then there are two volumes in $`S^3`$ that are bounded by these faces. One of the volumes includes the north pole (for example), and the other does not. If one of the volumes is measured to be a positive $`V`$ (with $`V<1`$), then the other volume will be $`V1`$. (Careful attention to the orientation will make the latter negative.) In order to assign a unique topological density to a lattice field configuration, we define a tetrahedron’s volume to satisfy $`|V|<\frac{1}{2}`$.

We have solved the problem of uniqueness, but not that of conservation of the winding number. Consider a field configuration wherein one tetrahedron in $`S^3`$ has volume $`\frac{1}{2}ϵ`$, where $`ϵ`$ is small. Under a fluctuation of one of the vertices of the tetrahedron, its volume may shift to $`\frac{1}{2}+\delta `$, while the neighboring tetrahedra change their volumes by $`(ϵ+\delta )`$ so that the winding number is unchanged. Unfortunately, our algorithm will now redefine the volume of the first tetrahedron to $`\frac{1}{2}+\delta `$, resulting in a loss of 1 in the winding number. This is often called “topology dropping through the lattice.”

The procedure to be followed at this point depends on the physics to be investigated. One might want to study, for example, the thermodynamics of the Skyrme model by fixing a chemical potential $`\mu `$ coupled to the winding number. The grand partition function would then be

$$Z(\mu )=\underset{n}{}e^{\mu n}Z_n,$$
(12)
where $`Z_n`$ is a sum over all field configurations with winding number $`n`$. A Monte Carlo simulation should then be allowed to wander freely among the sectors of different $`n`$, subject to an acceptance/rejection test that enforces the relative probabilities $`e^{\mu n}`$. A local updating algorithm can nucleate “point skyrmions” as in the preceding paragraph, which change $`n`$ by a unit and then spread out into smoother skyrmions; in addition, a non-local update could be permitted that creates smooth skyrmions directly.

Our interest, however, is in the properties of a single skyrmion, which are accessible through a canonical ensemble at fixed $`n=1`$. The simulation must be constrained so as not to change $`n`$. This requires that we prevent the nucleation of both point skyrmions and smooth skyrmions. The latter can be prevented by choosing a local updating algorithm. For the former, we choose a lattice action that excludes the possibility of tetrahedra in $`S^3`$ with $`|V|\frac{1}{2}`$. We adopt the kinetic term

$$S_1=(\alpha 1)\underset{𝐧\mu }{}\mathrm{log}(𝝈_𝐧𝝈_{𝐧+\widehat{\mu }}\alpha )$$
(13)
in order to constrain $`𝝈_𝐧𝝈_{𝐧+\widehat{\mu }}>\alpha `$ and thus to reject fluctuations that put large angles between neighboring field variables. With some Monte Carlo exploration, we can find a value of $`\alpha `$ that will keep the tetrahedral volumes far enough from $`\frac{1}{2}`$.

For smooth configurations, we expand $`𝝈_𝐧𝝈_{𝐧+\widehat{\mu }}`$ around 1 and Eq. (13) becomes

$`S_1`$ $``$ $`{\displaystyle \underset{𝐧\mu }{}}(1𝝈_𝐧𝝈_{𝐧+\widehat{\mu }})`$ (14)
$`=`$ $`{\displaystyle \underset{𝐧\mu }{}}\left[1{\displaystyle \frac{1}{2}}\text{Tr}\left(U_𝐧U_{𝐧+\widehat{\mu }}^{}\right)\right]`$ (15)
$`=`$ $`{\displaystyle \underset{𝐧\mu }{}}{\displaystyle \frac{1}{4}}\text{Tr}\left(U_{𝐧+\widehat{\mu }}U_𝐧\right)\left(U_{𝐧+\widehat{\mu }}^{}U_𝐧^{}\right),`$ (16)
which approaches the kinetic term in Eq. (1) in the continuum limit.

A technical point remains. We have discussed the ambiguity in fixing the volume of a tetrahedron if its faces are given, which led us to require $`|V|<\frac{1}{2}`$. A lattice field configuration $`𝝈_𝐧`$ gives us only the vertices of each tetrahedron, however, not its faces. Given three vertices that determine a great sphere, the face that connects them can be chosen to be either of two triangles that together make up the sphere. A little thought shows that the difference between the volumes enclosed is $`\frac{1}{2}`$. This gives an ambiguity between a value $`V>0`$ and $`V\frac{1}{2}`$ for the volume of the tetrahedron.

This last ambiguity involves choosing which triangle constitutes the face of the tetrahedron. But this face is shared between two adjacent tetrahedra, and the ambiguity will be immaterial if the same choice is made for both. This is hard to program, however, if each tetrahedron is to be handled on its own. We prefer to resolve the ambiguity for each tetrahedron separately, by requiring that $`|V|<\frac{1}{4}`$. We impose this via a more restrictive choice of $`\alpha `$ in Eq. (13). We find that setting $`\alpha =0.1`$ is adequate for the purpose.

A local updating algorithm based on the action (13) will conserve winding number. Even so, a skyrmion will tend to shrink down to a point in accordance with the scaling argument that gives Derrick’s Theorem. This is because the stabilization is being done at the scale of a single lattice spacing; only here will the scaling argument fail and the skyrmion run into a repulsive potential. In order to have a stable skyrmion of appreciable size, we cannot avoid adding a Skyrme term as in Eq. (1). The most straightforward lattice transcription uses symmetric derivatives ,

$$S_2^{\text{SYM}}=\underset{𝐧}{}\underset{\mu >\nu }{}\left\{\left(𝝈_{𝐧+\widehat{\mu }}𝝈_{𝐧\widehat{\mu }}\right)^2\left(𝝈_{𝐧+\widehat{\nu }}𝝈_{𝐧\widehat{\nu }}\right)^2\left[\left(𝝈_{𝐧+\widehat{\mu }}𝝈_{𝐧\widehat{\mu }}\right)\left(𝝈_{𝐧+\widehat{\nu }}𝝈_{𝐧\widehat{\nu }}\right)\right]^2\right\}.$$
(17)
This couples each site to sites two links away, effectively via plaquette couplings across plaquettes of side $`\sqrt{2}`$, and fails to couple the odd and even sublattices. We prefer to use a discretization  that couples only across single plaquettes, of side 1:

$$S_2=4\underset{𝐧}{}\underset{\mu >\nu }{}\left\{\left(𝝈_{𝐧+\widehat{\mu }}𝝈_{𝐧+\widehat{\nu }}\right)^2\left(𝝈_{𝐧+\widehat{\mu }+\widehat{\nu }}𝝈_𝐧\right)^2\left[\left(𝝈_{𝐧+\widehat{\mu }}𝝈_{𝐧+\widehat{\nu }}\right)\left(𝝈_{𝐧+\widehat{\mu }+\widehat{\nu }}𝝈_𝐧\right)\right]^2\right\}.$$
(18)
Our lattice action is a combination of the kinetic and Skyrme terms,

$$S=\beta _1S_1+\beta _2S_2.$$
(19)
In the naive continuum limit, we can compare with Eq. (1) to deduce that

$$\beta _1=\frac{f_\pi ^2}{4},\beta _2=\frac{1}{32e^2}.$$
(20)
## IV Topological density on the lattice

As shown in Fig. 1, we cut each cube on the lattice into five tetrahedra. Each tetrahedron maps into a curved tetrahedron in $`S^3`$, the volume of which gives the local topological density. We give here a practical, exact formula for this volume.

In the 2-sphere $`S^2`$, the area of a triangle is given simply by Girard’s theorem, $`A=\alpha +\beta +\gamma \pi `$, where the right hand side is the sum of the triangle’s angles minus $`\pi `$, known as the angular excess. The easiest proof of the theorem is an argument based on overlapping lunes. This argument may be generalized  to a simplex in $`S^n`$, but it only gives a formula for the volume when $`n`$ is even; for odd $`n`$, one obtains instead a constraint on the angles of the simplex.

L. Schläfli  attacked the general problem of the volume of convex polytopes in $`S^n`$. He was led to consider (for $`n=3`$) the special case of the quadrirectangular tetrahedron, which we shall abbreviate as qrt. A qrt (see Fig. 2) is constructed as follows .

Choose four points $`𝐏_0`$, $`𝐏_1`$, $`𝐏_2`$, $`𝐏_3`$ such that the line segments $`𝐏_0𝐏_1`$, $`𝐏_1𝐏_2`$, $`𝐏_2𝐏_3`$ are all mutually perpendicular. Then connect all four points together to form a tetrahedron. Every face of this tetrahedron is a right triangle. Three of the dihedral angles—the ones at $`𝐏_1𝐏_2`$, $`𝐏_1𝐏_3`$, and $`𝐏_0𝐏_2`$—are right angles. The dihedral angles at $`𝐏_2𝐏_3`$, $`𝐏_0𝐏_3`$, and $`𝐏_0𝐏_1`$ are not right angles, and they are denoted $`\alpha `$, $`\beta `$, $`\gamma `$.

This construction works equally well in Euclidean space and in $`S^3`$. A simple example of a spherical qrt is constructed by taking $`𝐏_0`$, $`𝐏_1`$, $`𝐏_2`$, $`𝐏_3`$ to be mutually perpendicular unit 4-vectors. This is a qrt that covers 1/16 of $`S^3`$. Its dihedral angles are all $`\pi /2`$.

In Euclidean space, the dihedral angles of a qrt satisfy

$$\mathrm{sin}\alpha \mathrm{sin}\gamma =\mathrm{cos}\beta .$$
(21)
In $`S^3`$, they satisfy instead the inequality

$$\mathrm{sin}^2\alpha \mathrm{sin}^2\gamma >\mathrm{cos}^2\beta .$$
(22)
The quantity

$$D=\sqrt{\mathrm{sin}^2\alpha \mathrm{sin}^2\gamma \mathrm{cos}^2\beta }$$
(23)
may thus be taken to be the generalization of the angular excess to this case. It vanishes as the qrt becomes small, which is the Euclidean limit.

Let us denote the volume of a spherical qrt as $`V(\alpha ,\beta ,\gamma )`$, normalized to $`2\pi ^2`$ for the entire 3-sphere. Schläfli derived formulas for the derivatives $`V/\alpha `$, $`V/\beta `$, $`V/\gamma `$ (see also ). The integral of these formulas was found by Coxeter  in the form of a Fourier series valid in the restricted domain

$`0`$ $`\alpha `$ $`{\displaystyle \frac{1}{2}}\pi `$ (24)
$`0`$ $`\beta `$ $`\pi `$ (25)
$`0`$ $`\gamma `$ $`{\displaystyle \frac{1}{2}}\pi .`$ (26)
Expressing the volume in terms of the complements of $`\alpha `$ and $`\gamma `$,

$$V(\alpha ,\beta ,\gamma )=\frac{1}{4}S(\frac{\pi }{2}\alpha ,\beta ,\frac{\pi }{2}\gamma ),$$
(27)
his solution for the Schläfli function $`S`$ is

$$S(x,y,z)=\underset{m=1}{\overset{\mathrm{}}{}}\left(\frac{D\mathrm{sin}x\mathrm{sin}z}{D+\mathrm{sin}x\mathrm{sin}z}\right)^m\frac{\mathrm{cos}2mx\mathrm{cos}2my+\mathrm{cos}2mz1}{m^2}x^2+y^2z^2.$$
(28)
Armed with Eq. (28), we can calculate the volume of any tetrahedron in $`S^3`$ by cutting it up into six qrt’s (see Fig. 3).

This is done by picking a vertex A of the tetrahedron and dropping a perpendicular to the opposite face (the base) at E, then dropping perpendiculars from E to the edges of the base. Simple trigonometry gives the vertices of the qrt’s in terms of those of the original tetrahedron; from the vertices of each qrt we calculate its dihedral angles $`\alpha ,\beta ,\gamma `$. Some technical problems, as well as further discussion of the Schläfli function, are relegated to the Appendix.

The topological density $`\rho _𝐧`$ in the lattice cube at n is the sum of the volumes of the tetrahedra in $`S^3`$ that correspond to the Euclidean tetrahedra making up the cube. It goes without saying that summing $`\rho _𝐧`$ over the lattice must give an integer (in units of $`2\pi ^2`$) for any configuration, to high precision.

## V Classical Skyrme model in 3 dimensions

The minima of the action (19) in 3 dimensions are the static, classical solutions of the 4-dimensional Skyrme model. When we include fluctuations about these minima, we can think of the 3d action that governs them as of an approximation to an effective action derived via dimensional reduction. Thus the parameters $`\beta _1,\beta _2`$ are (unknown) functions of the parameters of the full 4d theory, namely, $`f_\pi ,e`$, and the temperature $`T`$. Naturally, any reduction scheme that is precisely defined will give an effective action that is far more complex than (19), and further arguments will be necessary to justify its simplification.

As an application of our algorithm for the topological density we calculate the mean-square radius $`R^2`$ of a single skyrmion in the three-dimensional theory, that is, of an equilibrium distribution of configurations with $`n=1`$ at fixed $`\beta _1,\beta _2`$. We fix initial conditions of the form (6) and (10), suitably discretized. Monte Carlo updates are done with the usual local Metropolis algorithm, which preserves winding number as discussed above. The volume of the lattice is $`16^3`$.

The observable $`R^2`$ should be defined carefully. For each configuration used in the average one might calculate the barycenter of the topological density according to

$$𝐑_c=\underset{𝐧}{}𝐧\rho _𝐧,$$
(29)
and then the second moment via

$$R^2=\underset{𝐧}{}\text{min}(𝐧𝐑_c)^2\rho _𝐧.$$
(30)
(The notation “min” means that one should take note of the periodic boundary conditions and always calculate the shortest distance between n and $`𝐑_c`$.) Unfortunately, the boundary conditions make $`𝐑_c`$ calculated via Eq. (29) ill-defined. This may be illustrated by considering a compact skyrmion with barycenter located on the face of the lattice at $`z=0`$. Equation (29) fixes $`𝐑_c`$ in this case to be in the middle of the lattice, between the two half-skyrmions on opposite faces, and $`R^2`$ will turn out to be on the order of the lattice size. The solution lies in regarding the lattice as a 3-torus on which the choice of origin (and of “faces”) is entirely arbitrary. We evaluate Eqs. (29) and (30) for all choices of origin and take the minimal value of $`R^2`$. The $`𝐑_c`$ that corresponds to this minimum is thus defined to be the location of the barycenter; this is as good a definition as any, considering the periodic boundary conditions.

We show in Fig. 4 the calculated skyrmion radius for $`\beta _1=3`$ and various values of $`\beta _2`$.

The multiple values of $`R^2`$ at a given value of $`\beta _2`$ reflect metastability. These metastable skyrmion configurations are accessible from different initial states, specified by different values of $`r_0`$ in Eq. (10). We have observed tunneling from metastable states to lower-energy states of larger radius that are apparently the global minima of the action. In all cases the lowest-energy states appear to be those of largest radius; these are shown connected by the upper curve in Fig. 4. The lower curve in the figure connects metastable states that were found to be mutually accessible by annealing. We have not explored the full set of metastable states for any given coupling; except for the annealing curve, the metastable states shown in Fig. 4 are effectively chosen at random by the initial conditions of the simulations.

In Fig. 5 we present a more extensive set of results, this time showing variation with $`\beta _1`$ as well as with $`\beta _2`$.

The curves connect data points for what we believe to be the equilibrium skyrmion configuration at each coupling. Points not lying on curves are metastable configurations of higher energy; again, we make no effort to show all such configurations at any coupling. We have not run across metastable configurations for $`\beta _12`$.

The abscissa in Fig. 5 is the ratio $`\beta _2/\beta _1`$ that fixes the form of the action. The remaining overall coefficient $`\beta _1`$ then acts as an inverse temperature that governs the fluctuations. Regarding the 3d classical model as a Ginzburg-Landau theory for the fully quantized 4d theory, the lowering of $`\beta _1`$ can be interpreted as enhancing the thermal and quantum fluctuations in the latter. It is remarkable that heating the skyrmion by lowering $`\beta _1`$ causes it to shrink. We ascribe this surprising property to the many metastable, excited states that are smaller than the ground-state configuration and that become accessible as the temperature is raised. The absence of metastability at the smallest values of $`\beta _1`$ merely reflects the fact that the equilibrium skyrmion already averages over the various local minima.

## VI Discussion

The qualitative conclusion of the preceding section must be taken with a large grain of salt. As mentioned above, we have no idea what the connection is between the couplings $`\beta _1`$ and $`\beta _2`$ of the 3d model on the one hand and the couplings and temperature of the 4d theory on the other. It is possible that raising the physical temperature of the latter will lead to simultaneous changes in $`\beta _1`$ and $`\beta _2`$ that cause the skyrmion to expand after all. In any case, the result is not quantitative and cannot be made quantitative until the dimensional reduction to the 3d theory is explored in detail. The metastable solutions, for that matter, might be artifacts of the specific lattice action we study that are absent in the true effective 3d theory.

Our study of the 3d model was motivated by its simplicity rather than by any fundamental obstacle to studying the full 4d theory at finite temperature. Numerical work on the 4d theory will merely require greater computer resources. In order to extract physics from lattice calculations, we will have to renormalize the theory by fixing the bare couplings $`\beta _1`$ and $`\beta _2`$, and the lattice spacing $`a`$, in terms of physical quantities such as the $`\pi \pi `$ scattering length and the skyrmion mass and radius. The temperature can be varied by changing the time extent $`L_t`$ of the lattice in the time direction, or by varying the lattice spacing by changing $`\beta _1`$ and $`\beta _2`$ along lines of constant physics. The 4d quantum theory will offer as well the possibility of projecting the skyrmion to definite spin and isospin. Finally, it will be imperative to study the sensitivity of any physical quantities to the choice of lattice action.

## Acknowledgments

We thank M. Kugler and U.-J. Wiese for their help. Conversations with the latter took place at the Aspen Center for Physics. We also thank the Weizmann Institute of Science for its hospitality. A. J. S. thanks the Theoretical Elementary Particle Physics group at UCLA for its hospitality during his sabbatical leave. The work of B. S. was supported in part by the Israel Science Foundation under Grant No. 255/96-1. The work of A.J.S. was supported by the Research Corporation.

##

We collect in this appendix some practical details concerning the use of the Schläfli function (28) in calculating the topological density in lattice field configurations.

### 1 Small qrt’s

Defining

$$X\frac{\mathrm{sin}x\mathrm{sin}zD}{\mathrm{sin}x\mathrm{sin}z+D},$$
(31)
the formula for the Schläfli function is

$$S(x,y,z)=\underset{m=1}{\overset{\mathrm{}}{}}(X)^m\frac{\mathrm{cos}2mx\mathrm{cos}2my+\mathrm{cos}2mz1}{m^2}x^2+y^2z^2.$$
(32)
The sum converges absolutely when $`|X|1`$; since $`D>0`$, this means $`\mathrm{sin}x\mathrm{sin}z>0`$, which includes the entire domain (26) considered by Coxeter. In the Euclidean limit, where the qrt is small, we have $`D0`$ and hence $`X1`$. Then the sum in Eq. (32) becomes the Fourier series for $`x^2y^2+z^2`$, and thus $`S=0`$ in this limit as expected.

The region of small $`D`$ is where the series converges most slowly. We can avoid evaluating Eq. (32) in this region as follows . The vertices of a small tetrahedron in the 3-sphere are marked by 4-vectors $`𝝈^{(i)}`$ that are nearly equal. These 4-vectors are edges of a spindly 4-pyramid with its apex at the center of the 4-ball bounded by the 3-sphere, the volume of which is

$$V_4=\frac{1}{24}ϵ_{\mu \nu \rho \sigma }\sigma _\mu ^{(1)}\sigma _\nu ^{(2)}\sigma _\rho ^{(3)}\sigma _\sigma ^{(4)}.$$
(33)
The ratio of $`V_4`$ to the total volume $`\pi ^2/2`$ of the 4-ball is approximately the ratio of the 3-volume $`V`$ of the tetrahedron to the total volume $`2\pi ^2`$ of the 3-sphere, so

$$V4V_4.$$
(34)
We calculate $`V_4`$ for every tetrahedron obtained for each lattice cube (see Fig. 1). If it is small, we use Eq. (34) in lieu of cutting it into qrt’s. We also calculate $`V_4`$ for each qrt in turn, and if it small we use Eq. (34) in lieu of evaluating the Schläfli function. Evaluation of $`V_4`$ also gives an easy way to keep track of the sign of a qrt’s volume, which is needed since $`S(x,y,z)`$ is defined always to be positive.

### 2 Large tetrahedra

Another problem region for the Schläfli function arises when $`\alpha \frac{\pi }{2}`$ or $`\gamma \frac{\pi }{2}`$, giving $`X1`$ and slow convergence. In the limit $`X=1`$, we have

$$S(x,y,z)=\pi (x+yz)=\pi \left(\alpha +\beta +\gamma \pi \right),$$
(35)
a useful approximate formula in this regime.

As we have noted, Coxeter’s formula can only be used in the domain (26). One can easily devise larger qrt’s, however, and large qrt’s arise regularly in calculating the topological density. We attack this problem by cutting such large qrt’s into smaller tetrahedra. If either $`\alpha `$ or $`\gamma `$ is larger than $`\frac{\pi }{2}`$, but not both, we determine the largest of $`\alpha ,\beta ,\gamma `$ and draw a plane through its edge and the bisection point of the opposite edge (see Fig. 6).

This creates two new, smaller tetrahedra. If both $`\alpha `$ and $`\gamma `$ are larger than $`\frac{\pi }{2}`$, we cut the tetrahedron by connecting all its edge centers; this gives four tetrahedra at the original vertices, plus a central octahedron. The octahedron can then be cut into four tetrahedra around any line connecting opposite vertices, of which we choose the shortest.

Having cut the original qrt into two or eight pieces, we apply the algorithm (cut into qrt’s and apply the Schläfli function) to each piece. Of course, any qrt that results might still be too large for Coxeter’s formula. We then apply the cutting algorithm recursively. We have found that ten levels of recursion might be necessary for particularly rough lattice configurations.

The amount of recursion can be cut down dramatically by optimizing the operation of cutting the original tetrahedron into qrt’s. There are four ways to choose the vertex A in Fig. 3. For each choice of A, we determine the longest edge of the generated qrt’s, measured by the largest pairwise angle among the 4-vectors $`𝝈^{(i)}`$ of each qrt. This maximal edge is a measure of the size of the qrt’s generated by the cutting process. We choose A in order to minimize the maximal edge.