Semiclassical Study of Baryon and Lepton Number Violation in High-Energy Electroweak Collisions

TL;DR

This paper uses a semiclassical approach to calculate the suppression of baryon and lepton number violation in high-energy electroweak collisions, finding it remains exponentially suppressed up to very high energies and that low-energy approximations are valid only up to the sphaleron energy.

Contribution

It introduces a computational semiclassical method to evaluate suppression exponents for topology-changing processes in electroweak collisions, extending understanding to higher energies.

Findings

Baryon and lepton number violation remains exponentially suppressed up to at least 250 TeV.

Low energy expansion approximations are valid only up to 8 TeV.

The method provides results for s-wave scattering across a wide initial data range.

Abstract

We make use of a semiclassical method for calculating the suppression exponent for topology changing transitions in high-energy electroweak collisions. In the Standard Model these processes are accompanied by violation of baryon and lepton number. By using a suitable computational technique we obtain results for s-wave scattering in a large region of initial data. Our results show that baryon and lepton number violation remains exponentially suppressed up to very high energies of at least 30 sphaleron masses (250 TeV). We also conclude that the known analytic approaches inferred from low energy expansion provide reasonably good approximations up to the sphaleron energy (8 TeV) only.

Equations

$${\sigma }_{tot}(E)\sim \exp \left\{-\frac{16{\pi }^2}{g^2}{F}_{HG}(E/E_{\mathrm{sph}})\right\},$$

$$\sigma (E,N)\sim \exp \left\{-\frac{16{\pi }^2}{g^2}F(E/E_{\mathrm{sph}},\tilde{N})\right\}.$$

$$\underset{\tilde{N}\to 0}{\lim }F(E/E_{\mathrm{sph}},\tilde{N})=F_{HG}(E/E_{\mathrm{sph}}).$$

$$S=\frac{1}{4\pi {\alpha }_W}\int d^{4}x\{$$

$$+(D_{\mu }\Phi {)}^{†}{D}^{\mu }\Phi -\lambda ({\Phi }^{†}\Phi -1{)}^2\},$$

$$F_{\mu \nu }$$

$$D_{\mu }\Phi$$

$$M_W=\frac{1}{\sqrt{2}},$$

$$M_H=\sqrt{8\lambda }{M}_W.$$

$$\sigma (E,N)=\sum _{i,f}|⟨f|\hat{S}{\hat{P}}_E{\hat{P}}_N|i⟩{|}^2,$$

$$\sigma (E,N)\sim \exp \left\{-\frac{4\pi }{{\alpha }_W}F(\tilde{E},\tilde{N})\right\}$$

$$\displaystyle\frac{4\pi}{\alpha_{W}}F(\tilde{E},\tilde{N})=2\mathop{\mathrm{missing}}{Im}\nolimits S_{ABCD}(\varphi)-N\theta-ET-\mathop{\mathrm{missing}}{Re}\nolimits\mathcal{B}_{i}\;.$$

$${\mathcal{B}}_i=\frac{1}{2}\int d\mathbf{k}(f_{\mathbf{k}}^{}{f}_{-\mathbf{k}}^{}{\mathrm{e}}^{-2i{\omega }_{\mathbf{k}}(T_i-iT/2)}-g_{\mathbf{k}}^{\ast }{g}_{-\mathbf{k}}^{\ast }{\mathrm{e}}^{2i{\omega }_{\mathbf{k}}(T_i-iT/2)})$$

$$\displaystyle\varphi({\bm{x}},t)\big{|}_{t\to-\infty+iT/2}=$$

$$\displaystyle\int\frac{d{\bm{k}}}{(2\pi)^{3/2}{\sqrt{2\omega_{{\bm{k}}}}}}\Big{(}$$

$$\displaystyle+g_{{\bm{k}}}^{*}\mathrm{e}^{i\omega_{\bm{k}}(t-iT/2)-i{\bm{kx}}}\Big{)}\;.$$

$$\frac{\delta S}{\delta \phi }=0$$

$$f_{\mathbf{k}}^{}={\mathrm{e}}^{-\theta }{g}_{\mathbf{k}}.$$

$$\mathop{\mathrm{missing}}{Im}\nolimits\dot{\varphi}(\bm{x},T_{f}\to\infty)\to 0\;,\quad\mathop{\mathrm{missing}}{Im}\nolimits\varphi(\bm{x},T_{f}\to\infty)\to 0$$

$$E$$

$$N$$

$$\frac{4\pi }{{\alpha }_W}F(E)=\frac{4\pi }{{\alpha }_W}\left[1-\frac{9}{8}{\left(\frac{E}{E_0}\right)}^{4/3}+\frac{9}{16}{\left(\frac{E}{E_0}\right)}^2\right],$$

$$\sigma(E,N)=\int\!d\vartheta d\mathcal{T}da^{\vphantom{*}}_{\bm{k}}da^{*}_{\bm{k}}db^{\vphantom{*}}_{\bm{k}}db^{*}_{{\bm{k}}}d\varphi(x)d\varphi^{\prime}(x)\,\exp\Bigl{\{}-iN\vartheta-iE\mathcal{T}-\int d\bm{k}a^{\vphantom{*}}_{\bm{k}}a^{*}_{{\bm{k}}}\mathrm{e}^{-i\vartheta-i\omega_{\bm{k}}\mathcal{T}}-\int d\bm{k}b^{\vphantom{*}}_{\bm{k}}b^{*}_{{\bm{k}}}\\ +B_{i}(a_{{\bm{k}}},\varphi_{i})+B_{f}(b^{*}_{{\bm{k}}},\varphi_{f})+B^{*}_{i}(a^{*}_{-{\bm{k}}},\varphi_{i}^{\prime})+B^{*}_{f}(b_{-{\bm{k}}},\varphi_{f}^{\prime})+iS(\varphi)-iS(\varphi^{\prime})\Bigr{\}}\;.$$

$$B_i(a_{\mathbf{k}},{\phi }_i)$$

$$B_f(b_{\mathbf{k}}^{\ast },{\phi }_f)$$

$$\sigma(E,N)=\int d\vartheta d\mathcal{T}d\varphi(x)d\varphi^{\prime}(x)\prod\limits_{{\bm{k}}}\delta(\varphi_{f}(\bm{k})-\varphi^{\prime}_{f}(\bm{k}))\times\exp\biggl{\{}-iN\vartheta-iE\mathcal{T}+iS(\varphi)-iS(\varphi^{\prime})\\ -\frac{1}{2}\int d\bm{k}\frac{\omega_{{\bm{k}}}}{1-\gamma_{{\bm{k}}}^{2}}\Bigl{(}(1+\gamma_{{\bm{k}}}^{2})[\varphi_{i}(\bm{k})\varphi_{i}(-\bm{k})+\varphi^{\prime}_{i}(\bm{k})\varphi^{\prime}_{i}(-\bm{k})]-4\gamma_{{\bm{k}}}\varphi_{i}(\bm{k})\varphi^{\prime}_{i}(-\bm{k})\Bigr{)}\biggl{\}}\;,$$

$${\gamma }_{\mathbf{k}}={\mathrm{e}}^{i\vartheta +i{\omega }_{\mathbf{k}}\mathcal{T}}.$$

$$\frac{\delta S}{\delta \phi }=\frac{\delta S}{\delta {\phi }^{\prime }}=0,$$

$$\dot{\phi }(T_f,\mathbf{x})$$

$$\phi (T_f,\mathbf{x})$$

$$i{\dot{\phi }}_i(\mathbf{k})+{\omega }_{\mathbf{k}}{\phi }_i(\mathbf{k})$$

$$-i{\dot{\phi }}_i(\mathbf{k})+{\omega }_{\mathbf{k}}{\phi }_i(\mathbf{k})$$

$$\phi (x)$$

$${\phi }^{\prime }(x)$$

$$f_{\mathbf{k}}^{}$$

$$g_{\mathbf{k}}^{\ast }$$

Full text

††preprint: INR/TH-2003-5, BUHEP-03-09

I Introduction

Non-perturbative phenomena related to tunneling are often encountered in quantum field theory. Well known examples are false vacuum decay and instanton-like transitions (the latter are accompanied by non-conservation of fermion quantum numbers). When these phenomena are governed by a small coupling constant they can generally be studied by semiclassical methods. This is certainly the case at low energy or in situations which involve large number of quanta in the initial state. At low energy, collision processes can be well described by a semiclassical approximation relying on the existence of classical Euclidean time solutions to the equations of motion interpolating between initial and final states. In the examples mentioned above these are bounce Coleman:1977py and instanton Belavin:1975fg , respectively. The probability of the process is then proportional to the exponent of the Euclidean action of these solutions. As the action is inversely proportional to the (small) coupling constant, the processes are highly suppressed. The effect of low energy excitations in the initial state (colliding particles) gives only a pre-exponential factor and is inessential.

The situation changes at high energy, namely, at energy of the order of the tunneling barrier height which separates initial and final states. In general, there exists a static unstable solution to the equations of motion that lies on top of the potential barrier Klinkhamer:1984di (properly speaking, at a saddle point of the potential). In field theory this solution is often referred to as the “sphaleron”, a name which we will use throughout this paper. The minimum height of the barrier is precisely the sphaleron energy $E_{\mathrm{sph}}$. Naively, from the analogy with quantum mechanics of one degree of freedom, one would expect that at energy higher than the sphaleron energy the exponential suppression disappears. This is indeed what happens at finite temperature Kuzmin:1985mm ; Arnold:1987mh ; Arnold:1988zg ; Bochkarev:1987wg ; Khlebnikov:1988sr ; Grigorev:1989je ; Kuznetsov:1997sf ; Frost:1999eh ; Bonini:2000mb , finite fermion density Rubakov:1985am ; Rubakov:1986nk ; Matveev:1986bw ; Matveev:1987gq ; Diakonov:1992mp , or in the presence of heavy fermions in the initial state Rubakov:1985ix ; Ambjorn:1985bb ; Rubakov:1985it . But in high energy particle collisions this is not necessarily the case, due to the fact that the characteristic size of the sphaleron configuration is much larger than the wavelength of the incoming particles. At the same time the application of a semiclassical technique becomes problematic because the initial state no longer involves a large number of quanta.

As was first noted in Refs. Ringwald:1990ee ; Espinosa:1990qn , at relatively low energy the corrections to the collision-induced tunneling rate can be calculated by perturbative expansion in the background of the instanton. Further studies showed that the actual expansion parameter in most models, including electroweak theory, is $E/E_{\mathrm{sph}}$ McLerran:1990ab ; Khlebnikov:1991ue ; Yaffe:1990iy ; Arnold:1990va and the total cross section of the induced tunneling has an exponential form

[TABLE]

where $g$ is the small coupling constant and the function111The subscript $HG$ here stands for “holy grail” Mattis:1992bj . $F_{HG}(E/E_{\mathrm{sph}})$ is a series in fractional powers of $E/E_{\mathrm{sph}}$ (for a review see Mattis:1992bj ; Tinyakov:1993dr ; Rubakov:1996vz ).

While the perturbation theory in $E/E_{\mathrm{sph}}$ is limited to small $E$, the general exponential form of the total cross section implies that there might exist a semiclassical-type procedure which would allow, at least in principle, to calculate $F_{HG}(E/E_{\mathrm{sph}})$ at $E\gtrsim E_{\mathrm{sph}}$. However, since the initial state of two highly energetic particles is not semiclassical, the standard semiclassical procedure does not apply and a suitable generalization is needed, which was proposed in Refs. Rubakov:1992fb ; Tinyakov:1992fn ; Rubakov:1992ec and further developed in Refs. Kuznetsov:1997az ; Bezrukov:2001dg . The corresponding formalism reduces the calculation of the exponential suppression factor to a certain classical boundary value problem, whose analytical solution is not usually possible.

This approach is based on the conjecture that, with exponential accuracy, the two-particle initial state can be substituted by a multiparticle one provided that the number of particles is not parametrically large (although not proven rigorously, this conjecture was checked in several orders of perturbation theory in $E/E_{\mathrm{sph}}$ in gauge theory Tinyakov:1992fn ; Mueller:1993sc and explicitly in quantum mechanics with two degrees of freedom Bonini:1999cn ; Bonini:1999kj ). The few-particle initial state, in turn, can be considered as a limiting case of a truly multiparticle one with the number of particles $N=\tilde{N}\cdot(4\pi/g^{2})$, when the parameter $\tilde{N}$ is sent to zero. For the multiparticle initial state the transition rate is explicitly semiclassical and has the form

[TABLE]

According to the above conjecture, the function $F_{HG}(E/E_{\mathrm{sph}})$, corresponding to the two-particle incoming state, is reproduced in the limit $\tilde{N}\to 0$,

[TABLE]

Therefore, albeit indirectly, the function $F_{HG}(E/E_{\mathrm{sph}})$ is also calculable semiclassically.

Within the semiclassical framework, the function $F(E/E_{\mathrm{sph}},\tilde{N})$ is determined by the action evaluated at a particular solution to the classical field equations Rubakov:1992ec on a certain contour in complex time. In this formulation, the problem, at least in principle, is amenable to a computational solution. Namely, one has to solve the corresponding classical boundary value problem numerically and calculate the function $F(E/E_{\mathrm{sph}},\tilde{N})$, which then can be used to extract information about $F_{HG}(E/E_{\mathrm{sph}})$.

The implementation of this technique is nevertheless highly non-trivial. The differential equations one encounters are partially of the hyperbolic type (along the Minkowskian parts of the time contour) and partially of elliptic type (along the Euclidean part), which makes their numerical solution particularly challenging. In the electroweak theory, additional difficulties arise from the need to deal with the large number of internal degrees of freedom, unphysical modes due to gauge invariance, and time translational symmetry which cause an unwelcome degeneracy in the numerical procedure used to find the semiclassical solutions.

Moreover, at high energy (roughly, at energy higher than the sphaleron energy) tunneling solutions, interpolating between vicinities of different vacua in finite time, cease to exist. This subtlety turned out to be a general problem in the description of tunneling in systems with many degrees of freedom, and it has to do with the nontrivial way tunneling occurs at high energy—the system prefers to create a state close to the sphaleron, which then decays into the correct vacuum. To find the corresponding suppression exponent numerically one has to use a properly regularized version of the boundary value problem, developed in Ref. Bezrukov:2003yf .

In a long program of investigations we have been able to gradually overcome all of these hurdles. Preliminary results for energies below the sphaleron energy were reported in Ref. Bezrukov:2001dg . Now we are in a position of presenting what we are confident is the full solution to the numerical problem for a wide range of energy and incoming particle number, including energies above the sphaleron. The field configurations we analyse in this paper are restricted to be spherically symmetric in space. Hence our results apply, strictly speaking, to s-wave scattering only.

In this paper we will concentrate on obtaining the suppression exponent for collision-induced tunneling in $SU(2)$ gauge model with the Higgs mechanism, corresponding to the electroweak sector of the Standard Model at zero $\theta_{W}$. The problem is particularly interesting because of the baryon number violation which accompanies such processes 'tHooft:1976fv and the relatively low sphaleron energy $E_{\mathrm{sph}}\simeq 8\mathrm{TeV}$. Though computational limitations do not allow to reach literally zero value of the rescaled number of particles $\tilde{N}$, corresponding to particle collisions, we were able to extrapolate the results to zero $\tilde{N}$ and get a bound on the suppression exponent (strictly speaking, for s-wave scattering) and also provide an estimate for this exponent.

In section II we present the detailed formulation of the problem, outline the method and present the main physical results. In section III we give the derivation of the semiclassical method for the gauge model. The lattice formulation of the equations and subtleties appearing in the discretized version are given in section IV. Application of the regularization method of Ref. Bezrukov:2003yf is described in section V. Detailed numerical results are presented in section VI. Our conclusions are in section VII.

II Formulation of the problem and main results

Non-Abelian gauge models have an infinite number of topologically distinct vacua, labeled by an integer topological number. Processes changing the topological number are accompanied by violation of fermion (baryon and lepton) numbers 'tHooft:1976fv , a phenomenon of great interest for cosmology and particle physics. The topologically distinct vacua are separated by a potential barrier, whose height, in models with the Higgs mechanism, is given by the sphaleron energy. Topology changing transition may occur via tunneling at low energies or, at sufficiently high energy and suitable initial state, via classical evolution over the sphaleron.

In this paper we study a four-dimensional model which captures all the important features of the Standard Model—an $SU(2)$ gauge theory with the Higgs doublet. This model corresponds to the bosonic sector of the Standard Model with $\theta_{W}=0$. To the leading order in the coupling constant, the effect of fermions on the gauge and Higgs fields dynamics can be ignored Espinosa:1992vq . The action of the model is

[TABLE]

where

[TABLE]

with $A_{\mu}=A_{\mu}^{a}\sigma^{a}/2$ and $\alpha_{W}=g^{2}/4\pi$. Here we have eliminated inessential constants by an appropriate choice of units. The dimensional parameters can be restored noting that in the normalization (2), the gauge boson mass is

[TABLE]

and the Higgs boson mass is

[TABLE]

In most of our calculations the Higgs self-coupling $\lambda$ was set equal to $\lambda=0.125$, which corresponds to $M_{H}=M_{W}$. The dependence on the Higgs boson mass is very weak, so this is a reasonable approximation. Also we often omit the omnipresent overall factor $1/\alpha_{W}$.

Any vacuum configuration in this model can be obtained from the trivial vacuum $A_{\mu}=0$, $\Phi=\Phi_{v}=\begin{pmatrix}0\\ 1\end{pmatrix}$ by a certain gauge transformation $U(x)$. We will be using the temporal gauge $A_{0}=0$, where the vacuum configurations are described by time-independent $U(\bm{x})$, corresponding to residual gauge invariance. In this gauge, field values at spatial infinity cannot change during the evolution (otherwise the kinetic term becomes infinite) and thus one considers only those $U(\bm{x})$ which have some fixed asymptotics at spatial infinity. Often the asymptotic $U(\bm{x}\to\infty)\to 1$ is used, so any vacuum configuration corresponds to a mapping from space $R^{3}$ with identified infinity, which is homotopically equivalent to $S^{3}$, to the gauge group $SU(2)\sim S^{3}$. The degree of this mapping is precisely the topological number of the corresponding vacuum. A gauge choice of this form is convenient for analysis of the excitations about the trivial vacuum. For other purposes it may however be useful to choose an alternative behavior of the gauge function at spatial infinity, like $U(\bm{x})\to\exp\{i\boldsymbol{\bm{\sigma x/|x|}}\}$, which maps the $S^{2}$ of spatial infinity to the equatorial $S^{2}$ of the $SU(2)$. The two neighboring vacua then map the space $R^{3}$ either to north or south hemisphere of the $SU(2)$. In this gauge, the sphaleron configuration takes the simplest form, and we will use this gauge everywhere in this paper, except for the analysis of the mode expansion in the initial state.

Numerous perturbative attempts were made to find the probability of the collision-induced topology changing transitions in this model (see Refs. Mattis:1992bj ; Tinyakov:1993dr ; Rubakov:1996vz for reviews), giving reliable results only for relatively low energies. A non-perturbative study of classically allowed over-barrier transitions was presented in Ref. Rebbi:1996zx . All solutions found in Ref. Rebbi:1996zx are configurations with large numbers of particles in the initial state and thus they do not correspond to realistic collisions. Another approach, pursued in this paper, is to use the semiclassical method of Refs. Rubakov:1992fb ; Tinyakov:1992fn ; Rubakov:1992ec ; Kuznetsov:1997az adapted for theories with gauge degrees of freedom. This method was implemented in Ref. Bezrukov:2001dg , where the results were obtained for energies below $E_{\mathrm{sph}}$, what suggested that at the sphaleron energy the suppression is still strong. However, a straightforward application of the technique of Refs. Rubakov:1992fb ; Tinyakov:1992fn ; Rubakov:1992ec ; Kuznetsov:1997az fails for energy above the sphaleron due to the problems one encounters as the energy approaches the height of the barrier in systems with many degrees of freedom. These problems were studied in detail in the context of a quantum mechanical model in Ref. Bezrukov:2003yf , where a regularization technique was suggested to overcome them.

The basic idea in the proposal of Refs. Rubakov:1992fb ; Tinyakov:1992fn ; Rubakov:1992ec ; Kuznetsov:1997az is that, instead of a process with exclusive, two-particle initial state, one considers a topology changing process with inclusive initial state characterized by definite energy $E$ and incoming particle number $N$. The transition probability $\sigma(E,N)$ can then be used to provide a bound on the exclusive two-particle cross-section, while the two-particle transition exponent is obtained in the limit $\alpha_{W}N\to 0$.

The inclusive probability of tunneling from a state with fixed energy and number of particles is

[TABLE]

where $\hat{S}$ is the $S$-matrix, $\hat{P}_{E,N}$ are projectors onto subspaces of fixed energy $E$ and fixed number of particles $N$, and the states $|i\rangle$ and $|f\rangle$ are perturbative excitations about topologically distinct vacua. This matrix element can be written in double path integral representation. For large $N=\tilde{N}/\alpha_{W}$ and $E=\tilde{E}/\alpha_{W}$ the path integral can be calculated in the semiclassical approximation, and this leads to the problem of solving the equations of motion of the system on a special contour in complex time plane, which detours around singularities, as shown in Fig. 1. The presence of branch cut singularities can be inferred from the following argument. One notices that, for energy below the sphaleron energy, if one continues the solution along a line parallel to the real axis, be this via a forward integration of the equation of motion from the AB part of the contour or a backward integration from the CD part of the contour, the field must fall back to the original topological sector. On the other hand, by construction, on the AB and CD parts of the contour the solution must be in different topological sectors. Thus the solution must also be in different topological sectors on the AB part of the contour and on the negative real axis and, likewise, on the positive real axis and the continuation of the AB segment to positive time. This may happen only if two branch cut singularities exist on the two sides of the BC part of the contour (see Fig. 1).

Eventually the semiclassical approximation produces the following result ($\varphi$ here stands for all physical fields in the model),

[TABLE]

Here $S_{ABCD}(\varphi)$ is the action along the time contour, the parameters $T$ and $\theta$ are Legendre conjugate to $E$ and $N$; the parameter $T$ is the same as in Fig. 1; we will have to say more about $\theta$ later on. In what follows we will usually drop the tilde over the rescaled energy and incoming particle number, and the overall $1/\alpha_{W}$ factor, restoring it only in the final results. The boundary term

[TABLE]

is written using frequency components $f^{\vphantom{*}}_{\bm{k}}$ and $g^{\vphantom{*}}_{\bm{k}}$ of the field on the part A of the contour:

[TABLE]

The field $\varphi$ satisfies the field equation

[TABLE]

On the final part of the contour (CD), the field must satisfy the reality condition

[TABLE]

(for complex fields, such as $\Phi$ in (2), this means that both $(\Phi+\Phi^{\dagger})/2$ and $(\Phi-\Phi^{\dagger})/2i$ must be real).

Equations (9a)–(9c) specify the boundary value problem corresponding to the induced topological transition.

The equations obtained by variation over the auxiliary parameters $T$ and $\theta$ are

[TABLE]

These equations indirectly fix values of $T$ and $\theta$ for given energy and number of particles. Alternatively, one can fix $T$ and $\theta$, solve the boundary value problem (9) and obtain the corresponding values of $E$ and $N$ using (10) and (11). This is especially convenient in numerical calculations.

The interpretation of the solutions to the boundary value problem (9) is as follows. On the part CD of the contour, the saddle-point field is real asymptotically; it describes the evolution of the system after tunneling. On the contrary, it follows from boundary conditions (9b) that in the initial asymptotic region the saddle-point field is complex whenever $\theta\neq 0$. Thus, the initial state which maximizes the probability (7) is not described by a real classical field, i.e. this stage of the evolution is essentially quantum even at $N\sim 1/\alpha_{W}$.

There is a subtle point concerning the boundary condition (9c). It can be satisfied in two different ways. Either the solution is exactly real on the whole CD part of the contour and is close to vacuum, or it has an exponentially decaying imaginary part and approaches the sphaleron along the complexified unstable direction. This subtlety is important for the analysis at high energies (section V), $E\gtrsim E_{\mathrm{sph}}$.

The solutions to the boundary value problem can be found numerically for different values of $E$ and $N$. In this paper we study solutions that have spherical symmetry in space. One expects that these are most important for large enough $N$; perturbative calculations about the instanton suggest that spatial spherical symmetry is relevant at relatively low energies and all $N$. We do not have a convincing argument in favor of spherical symmetry for few particle collisions at very high energies; in any case, our results as they stand, are valid for s-wave scattering.

Our numerical analysis shows that the $E-N$ plane is divided into several different regions (see Fig. 2). Values of $E<N\cdot\min(M_{W},M_{H})$ are trivially excluded by kinematics. For relatively low energies (region A) the transitions between the topologically distinct vacua can occur only via tunneling. At the sphaleron energy $E_{\mathrm{sph}}$ the situation changes. A slight excitation of the sphaleron along the unstable direction gives origin to a solution of the classical equations of motion which evolves towards different topological sectors at large negative and positive times. Since the sphaleron has exactly one negative mode, there is only one infinitesimal deformation of this type, and thus the corresponding solution has definite number of particles $N_{\mathrm{sph}}$ in the initial state. At higher energies one may add excitations of the positive modes above the sphaleron to obtain over-barrier solutions with different, and, in particular, smaller initial particle number. These solutions belong to the domain of classically allowed transitions (region B in Fig. 2), where the topology changing processes are unsuppressed. The boundary between region A and region B corresponds to configurations staying for an infinite time close to the sphaleron (since there are no bound states in the sphaleron background Akiba:1989xu , all excitations about the sphaleron fly away at finite time, so the field relaxes to the sphaleron solution).

In the classically forbidden region A there is a special family of solutions, corresponding to $\theta=0$ in the boundary value problem. These are represented by the line $E_{\mathrm{PI}}(N)$ in Fig. 2. In this case, the boundary condition (9b) reduces to the reality condition imposed at \mathop{\mathrm{missing}}{Im}\nolimits t=T/2. The solution to the resulting boundary value problem is the periodic instanton of Ref. Khlebnikov:1991th . The periodic instanton is a real periodic solution to the Euclidean field equations with period $T$ and turning points at $t=0$ and $t=iT/2$ $(\mathrm{mod}T)$. When analytically continued in the Minkowskian direction through the turning points, the periodic instanton is real on the lines \mathop{\mathrm{missing}}{Im}\nolimits t=0 and \mathop{\mathrm{missing}}{Im}\nolimits t=T/2 and therefore satisfies the boundary value problem (9) with $\theta=0$. Like any other solution linearizing at large negative times at part A of the contour of Fig. 1, the periodic instanton has a certain number of incoming particles, Eq. (11). For given energy $E$ below the sphaleron, this number is such that the suppression exponent $F(E,N)$ has a minimum, i.e., the transition occurs at maximum rate.

The classically forbidden region A is further subdivided into two regions. For low energies (region A.I) the system is close to the vacuum on the final part of the evolution, so the boundary condition (9c) leads to the exact reality of the fields on the part CD of the time contour. At energies higher than the sphaleron energy (precisely, on the right of the line $E_{1}(N)$) the system ends up close to the sphaleron (with extra outgoing waves in the sphaleron background). In this case Eq. (9c) is truly asymptotic. So, the system tunnels “on top” of the barrier, creating an unstable sphaleron configuration, which then decays with probability of order 1 to any of the two neighboring vacua. This situation is realized in the region A.II. This new qualitative feature of the tunneling at high energies emerges from the existence of the bifurcation of the solutions and is not seen in any order of perturbative expansion around the instanton. Non-perturbative approaches, however, capture this feature (see also Ref. Voloshin:1994dk for similar results in the context of false vacuum decay).

Our numerical results for the suppression exponent in the whole classically forbidden region are presented in Figs. 3, 4. The almost vertical line in Fig. 3 separates the two regions (denoted by A.I and A.II in our earlier discussion) where the tunneling process assumes characteristically distinct features. It also represents the frontier beyond which numerical calculations based on a straightforward implementation of the method of Refs. Rubakov:1992fb ; Tinyakov:1992fn ; Rubakov:1992ec ; Kuznetsov:1997az appear to fail. It is clear from Fig. 3 that our improved numerical technique can go well beyond that frontier. Reference Bezrukov:2003zn presents a comparison between our results and the analytic predictions for the suppression exponent $F(E,N)$ in the limit of small energy. The two are in remarkable agreement which provides a gratifying check of the numerical calculations.

Another interesting comparison can be made with the results of Ref. Rebbi:1996zx , where the real-time overbarrier solutions close to the boundary of the classically allowed region were searched via Monte-Carlo techniques222In Ref. Rebbi:1996zx the coupling constant $\lambda$ was chosen $0.1$, while we use $\lambda=0.125$. We performed a set of calculations for $\lambda=0.1$. The dependence on $\lambda$ is so weak, that the difference for the results would be invisible in the graph. Much larger discrepancies appear because of the different lattice parameters used in the two calculations. (In Ref. Rebbi:1996zx , having only to solve for the real time evolution of the fields, it was possible to use a larger lattice and a finer lattice spacing than in the present calculations.). In this way, an approximation (which, at the same time, is an upper bound) for the boundary of the classically allowed region was obtained. It is seen that the results of Ref. Rebbi:1996zx are reasonably close to the boundary $E_{0}(N)$, found in our calculations.

Our results by themselves do not reach the physically interesting $N=0$ limit (corresponding to particle collisions). Studying lower $N$ in numerical calculations would need lattices with larger number of lattice points and would require quite substantial amounts of time even on powerful present day supercomputers. Therefore, some extrapolation must still be used to get insight on the suppression factor for actual particle collisions. As we seek such extrapolations, we notice first that it is quite straightforward to obtain a lower bound on the suppression exponent $F$. Insofar as $\theta$ increases as $N\to 0$, and $(4\pi)\partial F/\partial N=-\theta$, by simply continuing $F$ with a linear function of $N$ for each energy one obtains a lower bound on $F$. This bound is shown in Figs. 5, 6, dashed line. It indicates that up to the energy $8M_{W}/\alpha_{W}\simeq 20\mathrm{TeV}$ the suppression is still high: the suppression factor is smaller than $\mathrm{e}^{-60}\sim 10^{-26}$ for $\alpha_{W}\sim 1/30$.

For very high energies a bound may be constructed by exploiting the observation that the lines of constant $F$ in $E-N$ plane have positive curvature (see Fig. 3). So, by extrapolating these lines linearly to $N=0$ one obtains another lower bound on the suppression exponent $F(E,N=0)$. This bound is displayed in Fig. 5, dashed-dotted line. One can see that exponential suppression continues up to an energy of at least $250\mathrm{TeV}$.

One may also attempt to estimate the function $F(E)$ itself. As we discuss in Section VI, a good estimate is obtained by extrapolating, instead of $F(E,N)$, the function $T(N)$ at fixed energy, as $T(N)$ is approximately linear in $N$. Up to the sphaleron energy, the estimate obtained in this way is close to the one loop analytic result Khoze:1991bm ; Arnold:1991cx ; Diakonov:LINPSchool1991 ; Mueller:1991fa , which gives three terms in the low-energy expansion,

[TABLE]

where $E_{0}=\sqrt{6}\pi M_{W}/\alpha_{W}$. Below the sphaleron, our estimate is also consistent with the analytic estimate of Refs. Ringwald:2002sw ; Ringwald:2003px . On the other hand, the behavior of $F_{HG}(E)$ changes dramatically at $E\gtrsim E_{\mathrm{sph}}$. We attribute this to the change in the tunneling behavior—at $E\gtrsim E_{\mathrm{sph}}$ the system tunnels “on top of the barrier”. Our numerical data show that the suppression exponent $F_{HG}(E)$ flattens out, and topology changing processes are in fact much heavier suppressed at $E\gtrsim E_{\mathrm{sph}}$ as compared to the estimate (12) and the estimate of Refs. Ringwald:2002sw ; Ringwald:2003px . We show our estimate, together with analytical estimates and our lower bound, in Fig. 6.

It is worth noting that similar effects of dramatic change of the behavior of the system at high energies were observed in lattice calculations of instanton distribution in QCD in Ringwald:1999ze ; Schrempp:2002kd .

Thus, our numerical results, albeit covering a limited range of energies and initial particle numbers, enable us to obtain both lower bound for and actual estimate of the suppression exponent for the topology changing two-particle cross-section in the electroweak theory well above the sphaleron energy. This cross section remains exponentially suppressed up to very high energies of at least $250\mathrm{TeV}$. In fact, the energy, if any, at which the exponential suppression disappears, is most likely much higher, as suggested by comparison of our lower bound and actual estimate at energies exceeding significantly $E_{\mathrm{sph}}$, see Fig. 6.

III The method of RST

III.1 General formulation

The quantity we wish to calculate is $\sigma(E,N)$, the probability of transition from a state with fixed energy $E$ and number of particles $N$ about one vacuum to any state about another vacuum. The method of semiclassical calculation of this inclusive multiparticle probability was formulated in Refs. Rubakov:1992fb ; Tinyakov:1992fn ; Rubakov:1992ec ; Kuznetsov:1997az . We call it the RST method for brevity, and here we review this prescription in brief.

The inclusive multiparticle probability (6) can be written in functional integral form, where the semiclassical approximation is equivalent to the saddle-point integration. The double path integral representation for $\sigma(E,N)$ reads Rubakov:1992fb

[TABLE]

Here $\varphi$ stands for all physical fields of the theory. The boundary terms $B_{i}$ and $B_{f}$ are

[TABLE]

where $\varphi_{i,f}(\bm{k})$ are the spatial Fourier transforms of the field at initial and final times $T_{i}$ and $T_{f}$, respectively. The limit $T_{i,f}\to\mp\infty$ is assumed at the end of the calculation. The complex integration variables $a_{{\bm{k}}}$ and $b^{*}_{{\bm{k}}}$ come from the coherent state representation of initial and final states; they are the classical counterparts of annihilation and creation operators. The integration over these variables implements the summation over initial and final states in (6). The functional integrals over $\varphi(x)$ and $\varphi^{\prime}(x)$ come from the amplitude and complex conjugate amplitude, respectively. The integrations include the boundary values $\varphi_{i,f}$ and $\varphi^{\prime}_{i,f}$. Integration over $\mathcal{T}$ and $\vartheta$ serve to project onto the subspaces of fixed $E$ and $N$, respectively.

The integral (13) can be evaluated in the saddle point approximation, as long as the exponent is proportional to $1/\alpha_{W}$, implicitly present in the expression, and $N,E\sim 1/\alpha_{W}$.

Let us now discuss the saddle-point equations for the integral (13). We will see that these equations reduce to a certain boundary value problem for the fields $\varphi$ and $\varphi^{\prime}$. The variables $a^{\vphantom{*}}_{\bm{k}}$, $a^{*}_{{\bm{k}}}$, $b^{\vphantom{*}}_{\bm{k}}$ and $b^{*}_{{\bm{k}}}$ enter the exponent quadratically and can be integrated out, yielding

[TABLE]

where

[TABLE]

An important feature of the representation (15) is that the exponent in the r.h.s. contains only the action and the boundary terms. Thus, the discretization of this exponent is relatively straightforward.

Let us turn to the saddle point equations. Varying the exponent with respect to the fields $\varphi(x)$ and $\varphi^{\prime}(x)$ we find

[TABLE]

i.e. the usual field equations. The boundary conditions for these equations come from the variation with respect to the boundary values of the fields. At $t=T_{f}$, because of the $\delta$–function, the variations are subject to the constraint $\delta\varphi_{f}(\bm{x})=\delta\varphi^{\prime}_{f}(\bm{x})$ (at $T_{f}\to\infty$). Since $\delta S/\delta\varphi(T_{f},\bm{x})=\dot{\varphi}(T_{f},\bm{x})$ we obtain

[TABLE]

Thus, in the final asymptotic region the saddle-point fields $\varphi$ and $\varphi^{\prime}$ coincide.

The variation with respect to $\varphi_{i}$ and $\varphi^{\prime}_{i}$ leads to two equations which can be written in the following form,

[TABLE]

These initial boundary conditions simplify when written in terms of frequency components. In the initial asymptotic region ($t\to-\infty$), where $\varphi$ and $\varphi^{\prime}$ are free fields, we can write

[TABLE]

Then the initial boundary conditions become

[TABLE]

Finally, there are two saddle-point equations which come from the variation of the exponent in (15) with respect to $\vartheta$ and $\mathcal{T}$. These equations determine the saddle-point values of $\vartheta$ and $\mathcal{T}$ as functions of $E$ and $N$. In terms of frequency components $f_{\bm{k}}$ and $g_{\bm{k}}$ they read (after using boundary conditions (19))

[TABLE]

One may recognize the usual expressions for the energy and the number of particles contained in the free classical field, $n^{\vphantom{*}}_{\bm{k}}=f^{\vphantom{*}}_{\bm{k}}g^{*}_{{\bm{k}}}$ being the occupation number in the mode with spatial momentum $\bm{k}$.

The field $\varphi^{\prime}(x)$ originates from the complex conjugate amplitude. This suggests that its saddle point value is complex conjugate to that of $\varphi(x)$. Indeed, the Ansatz

[TABLE]

is compatible with the boundary value problem (16)–(21). Then the saddle point values of $\mathcal{T}$ and $\vartheta$ are pure imaginary

[TABLE]

provided the initial energy (20) and particle number (21) are real. The boundary conditions (17) imply then that the field $\varphi$ is real asymptotically at final time

[TABLE]

while Eq. (19) relates the positive and negative frequency components of the field $\varphi$ in the initial asymptotic region

[TABLE]

where

[TABLE]

Until now, the initial time $T_{i}$ was real. However, it is convenient to reformulate the boundary value problem directly in terms of the fields on the contour ABCD, at which the initial time has imaginary part \mathop{\mathrm{missing}}{Im}\nolimits T_{i}=T^{\prime}/2 (see Fig. 7). The analytical continuation in the initial asymptotic region can be done explicitly by means of Eqs. (18). In Eqs. (19)–(21) this continuation results in the substitution of $\gamma_{{\bm{k}}}$ by

[TABLE]

The simplest boundary conditions are obtained in the case when the contour height in imaginary time $T^{\prime}$ is equal to the parameter $T$, leading to ${\bm{k}}$-independent $\gamma$

[TABLE]

In this case one arrives at the boundary condition (9b) and the contour ABCD with height $T/2$ shown in Fig. 1. This formulation will be used in most cases. Then the boundary value problem (16)–(21) is equivalent333The boundary term in (7) is obtained from the boundary terms in (15) by making use of the $\theta$ boundary conditions (9b). to (9)–(11). This is the boundary value problem we solve numerically in the present paper.

Let us discuss some subtle points of this boundary value problem. First, one notices that the condition of asymptotic reality (9c) does not always coincide with the condition of reality at finite time. Of course, if the solution approaches the vacuum on the part CD of the contour, the asymptotic reality condition (9c) implies that the solution is real at any finite positive $t$. Indeed, at large enough time the system evolves in that case in the linear regime, so the condition (9c) means that all physical modes should be real. Due to the equations of motion the fields are then real on the entire CD–part of the contour. This situation corresponds to the transition directly to the neighboring vacuum. However, the situation can be drastically different if the solution on the final part of the time contour remains in the interaction region, i.e. close to the sphaleron. Since one of the excitations about the sphaleron is unstable, there may exist solutions which approach the sphaleron exponentially along the complexified unstable direction. In that case the solution may be complex at any finite time, and become real only asymptotically, as $t\to+\infty$. Such solution corresponds to tunneling to the sphaleron; afterwards the system rolls down classically to the correct vacuum (with probability of order $1$, inessential for the tunneling exponent $F$). We will see in section V that the situation of this sort indeed takes place at high energies $E\gtrsim E_{\mathrm{sph}}$.

Second, the initial boundary conditions (19) (imposed on the real time axis) mean, that $\varphi$ and $\varphi^{\prime}=\varphi^{*}$ are different at large negative time, while at large positive time they coincide because of the condition (17). For solutions ending in the vacuum at positive time (so that the fields are exactly real at finite $t>0$), this means that there should exist a branch point in the complex time plane: the contour in Fig. 1 winds around this point and cannot be deformed to the real time axis. This argument does not work for solutions ending on the sphaleron at $t\to+\infty$, so branch points between the AB–part of the contour and the real time axis may be absent. We have found that this is indeed the case at high energies (cf. Ref. Bezrukov:2003yf ).

III.2 Reduction to spherically symmetric configurations

Here we consider spherically symmetric configurations Ratra:1988dp of the $SU(2)$–Higgs theory. The reason is that one can entertain the expectation that the most important tunneling configurations possess maximum spatial symmetry. On the other hand, without the simplification provided by spherical symmetry the computational cost of the numerical analysis would be prohibitive.

In the spherically symmetric Ansatz the original fields are expressed in terms of six real two-dimensional fields $a_{0}$, $a_{1}$, $\alpha$, $\beta$, $\mu$ and $\nu$ as follows

[TABLE]

where $\bm{n}$ is the unit three-vector in the radial direction and $\xi$ is an arbitrary constant two-component complex unit column. This Ansatz is symmetric under spatial rotations complemented by appropriate rotations in the gauge group and custodial global symmetry transformations. The action (2) expressed in terms of the new fields becomes

[TABLE]

where the indices $\mu$, $\nu$ run from 0 to 1 and

[TABLE]

Note that the overbar on $\phi$, $\chi$ and $D_{\mu}$ denotes changing $i\to-i$ in the definitions (25) above, which is the same as complex conjugation only if the six fields $a_{\mu}$, $\alpha$, $\beta$, $\mu$ and $\nu$ are real. In the boundary value problem (9) these fields become complex and overbar no longer corresponds to normal complex conjugation.

The equations of motion obtained from (24) are

[TABLE]

Equation (26a) is of the first order in time—it is Gauss’ law.

The spherical Ansatz (23) has a residual $U(1)$ gauge invariance

[TABLE]

with gauge function $\Omega(r,t)$. The complex “scalar” fields $\chi$ and $\phi$ have $U(1)$ charges $1$ and $1/2$ respectively. $a_{\mu}$ is the $U(1)$ gauge field, $f_{\mu\nu}$ is the field strength tensor, and $D_{\mu}$ in (25) is the covariant derivative. The residual $U(1)$ gauge invariance must be fixed when solving the equations numerically. We choose the temporal gauge $a_{0}=0$. In this gauge, if Gauss’ law is obeyed at some moment of time, the other five equations guarantee that it is obeyed at any time. This means, in fact, that one of the equations is redundant, and one of the fields is not physical—it can be expressed in terms of the other four fields and their derivatives using Gauss’ law. However, numerically it is easier to solve five second order equations of motion imposing Gauss’ law as one of the boundary conditions. Also, in the $a_{0}=0$ gauge, there remains a gauge freedom with time independent gauge function, and this should also be fixed by boundary conditions.

The trivial space-independent vacuum of the model is

[TABLE]

Other vacua are obtained from the trivial one by the gauge transformations

[TABLE]

By regularity, $\Omega$ should be zero at the origin. Vacua with different winding numbers correspond to $\Omega\to 2n\pi$ as $r\to\infty$. For such values of $\Omega$, the fields of the original four-dimensional model are constant at spatial infinity, and this is the standard choice. It allows for a simple description of the topological properties of vacua: since the sphere $S^{2}$ at spatial infinity is mapped to one point in field space, one can compactify the space to $S^{3}$ and consider mappings $S^{3}\to SU(2)$, corresponding to pure gauge field configurations.

One can also make other choice of gauge transformation function $\Omega(r)$ at spatial infinity (as long as the fields are pure gauge and constant in time there). In our case it is convenient to set $\Omega\to(2n-1)\pi$ at $r\to\infty$. This choice, called “periodic instanton gauge” in this paper, in the original 4-dimensional theory corresponds to mapping of the sphere $S^{2}$ at spatial infinity onto the equatorial sphere $S^{2}$ of the $SU(2)$ gauge group, parameterizing the pure gauge field configuration. This behavior of $\Omega$ is equivalent to the requirement that the fields satisfy the following boundary conditions at $r=0$ and $r=\infty$,

[TABLE]

The conditions for the field $\phi$ at $r\to 0$ make the original field $\Phi$ regular at the origin.

In this gauge no $r$-independent vacuum exists, but transition between vacua with $n=0$ and $n=1$ is described in a very symmetric way. The behavior of the fields $\chi$ and $\phi$ for such transition is shown in Fig. 8. In the original 4-dimensional model this topology changing process corresponds to a transition where the fields wind over the lower hemisphere of $SU(2)$ before the transition and over the upper hemisphere after the transition.

The initial $\theta$–boundary conditions in gauge theory are quite complicated. The basic reason is that there is a redundant field among the five fields $a_{1}$, $\phi$, $\bar{\phi}$, $\chi$, $\bar{\chi}$, while the $\theta$-boundary conditions (9b) are to be imposed on physical fields only. The analytic expressions for the modes $f_{\bm{k}}$, $g_{\bm{k}}$ in terms of the fields $\chi$, $\phi$, $a_{1}$ are cumbersome (see Refs. Rebbi:1996zx ; Farhi:1996aq ) and will not be presented here. It is simpler, and more precise in the lattice case, to perform this expansion numerically in the discretized version of the model. This expansion will be described in the following section.

To complete the boundary value problem, one has to impose Gauss’ law and the equation fixing the time independent gauge invariance. Note that both of these equations are not full complex valued equations (unlike the $\theta$–boundary conditions), otherwise the system would have been overdetermined. The point is that, the reality conditions at final time (9c) forbid gauge transformations with imaginary gauge functions and also guarantee that Gauss’ law does not have imaginary part. So, only the real part of Gauss’ law (26a) and equation fixing only real-valued gauge transformations must be used. Together with four $\theta$ boundary conditions this gives the right number of boundary conditions for the system with five complex valued fields $a_{1}$, $\alpha$, $\beta$, $\mu$, $\nu$. The exact form of the gauge fixing condition will be given in section IV, because it is again most conveniently expressed in lattice terms.

One more complication of the problem is the invariance of the equations under translations along the real time. To solve the equations numerically this should be fixed in a controlled way, to make sure the contour winds around the branching points of the solution, and does not get too close to them. A method of removing this invariance will also be described in Section IV.

IV Computational challenges

IV.1 Discretized action

To obtain a self-consistent system of equations, the discretization of the equations (26b)–(26f) should be done in a gauge invariant way.

First, let us consider the discretized version of the action (24). The spatial axis is discretized by introducing sites $r_{i}$, $i=0,\dots,N$, where $r_{0}=0$, $r_{N}=L$. The time grid consists of sites $t_{j}$, $j=-1,\dots,N_{t}+1$. We are working in the $a_{0}=0$ gauge, and omit the subscript in the spatial component of the gauge field, $a_{1}(r,t)\equiv a(r,t)$. Field variables $\chi_{ij}$, $\bar{\chi}_{ij}$ and $\phi_{ij}$, $\bar{\phi}_{ij}$ correspond to field values on the space-time lattice sites, while $a_{ij}$ are defined on spatial links and temporal sites. We also absorb the $\Delta r_{i}$ factors in the definition of $a_{ij}$. The boundary conditions in the periodic instanton gauge, Eq. (30), are

[TABLE]

for all $j$. In the boundary condition for $\phi_{0j}$, the spatial derivative in Eq. (30) was changed into a covariant one to preserve exact lattice gauge invariance. Thus the complex lattice field variables left are

[TABLE]

The discretized action reads

[TABLE]

where the weights are

[TABLE]

with $\Delta t_{j+1/2}=t_{j+1}-t_{j}$, $\Delta t_{j}=(\Delta t_{j+1}+\Delta t_{j})/2$, and analogous expressions for $\Delta r_{i}$; $h_{j}=1$ for $j=0,N_{t}$ and $1/2$ for $j=-1,N_{t}+1$. Lattice field equations are derived from (32) by variation over the lattice fields, after the boundary conditions (31) have been used to exclude $\chi_{0j}$, $\bar{\chi}_{0j}$, $\phi_{0j}$, $\bar{\phi}_{0j}$ and $\chi_{Nj}$, $\bar{\chi}_{Nj}$, $\phi_{Nj}$, $\bar{\phi}_{Nj}$ from the lattice action.

The action (32) is exactly invariant under time-independent lattice gauge transformations of the form

[TABLE]

This gauge freedom has to be fixed by boundary conditions.

IV.2 Boundary term: normal modes

To obtain lattice version of (9b) one notes that plane waves are no longer eigenfunctions of the Hamiltonian on the lattice. To find their analogue one brings the quadratic part of the action (32), taken in the limit of continuous time, to the canonical form. We expand it near the space-independent vacuum (28)

[TABLE]

(performing in the end a gauge transformation to the vacuum (30) is straightforward). It is also useful to change to the notations (25),

[TABLE]

In these terms, the quadratic part of the action (32) is

[TABLE]

As seen from Eq. (35), the variables $\tilde{\alpha}_{0}$, $\tilde{\beta}_{0}$, $\tilde{\mu}_{0}$ and $\tilde{\nu}_{0}$ do not have kinetic terms. Three of them are fixed by the boundary conditions at $r=0$,

[TABLE]

The fourth one, $\nu_{0}$, is determined from the field equation, which for this variable reads

[TABLE]

After the variables with $i=0$ have been excluded in the manner described above, the quadratic action takes the form

[TABLE]

where the real valued coefficients $d_{I}$ and $S_{IJ}$ are to be read off from Eq. (35), indices $I,J$ label fields and space points, and $\tilde{\varphi}_{I}$ stands for the fields $\{\tilde{\alpha}_{i},\tilde{\beta}_{i},\tilde{\mu}_{i},\tilde{\nu}_{i},a_{i}\}$. The change of variables

[TABLE]

brings the kinetic term to the canonical form,

[TABLE]

where

[TABLE]

The symmetric matrix $\tilde{S}_{IJ}$ is then diagonalized

[TABLE]

where $O_{KJ}$ is an orthogonal matrix. Introducing yet another set of variables $z_{I}$ by the relations

[TABLE]

we finally bring the action to the diagonal canonical form,

[TABLE]

Therefore, vectors

[TABLE]

are normal modes in the lattice formulation of the theory, and should be used instead of the usual spherical waves. The corresponding frequencies are $\omega^{2}_{K}$.

The matrix $O_{KJ}$ and frequencies $\omega_{K}$ are found numerically. Since they depend only on the spatial lattice parameters (size and spacing) and coupling constant $\lambda$, and do not depend on the background vacuum field configuration, it is sufficient to perform this diagonalization once for a given lattice. The first $4N-3$ eigenvectors $\xi_{(K)I}$ and eigenvalues $\omega_{K}$ correspond to physical modes, and the rest $N-1$ of them have $\omega=0$ and thus correspond to the gauge (unphysical) degrees of freedom.

IV.3 Boundary conditions.

$\theta$ boundary conditions.

To derive the lattice version of the boundary conditions, one takes the variation of the exponent for the total probability (15), which can be written in the following form

[TABLE]

where dots denote terms irrelevant in the current context and $\gamma=\mathrm{e}^{-\theta}$. One has to vary the lattice version of (36) with respect to $z_{I,-1}$ (values of $z$ at the first time slice) and set $z^{\prime}=z^{*}$. The variational equation reads

[TABLE]

which leads to

[TABLE]

where the derivatives of the action are equal to the classical momenta of the modes

[TABLE]

Here the index $I=1,\ldots,4N-3$ labels physical degrees of freedom. One can go back to the original notations by means of the relation

[TABLE]

where $I=1,\ldots,4N-3$.

Finally, we use the gauge transformation

[TABLE]

with $c=0.5$ to transform the fields from the periodic instanton gauge (30) to the form $\tilde{\varphi}_{J}$ (34) (of course, other gauge choices with $\Omega(0)=\pi$, $\Omega(L)=0$ are possible). This transformation has the form

[TABLE]

where $\varphi_{J}$ stands for the fields ${\alpha_{j},\beta_{j},\mu_{j},\nu_{j},a_{j}}$ in periodic instanton gauge. The matrix $g_{KJ}$ and vector $\tilde{\varphi}^{\mathrm{vac}}_{J}$ can be easily read off from the expression for the lattice gauge transformation (33) and definition (34).

Zero modes part.

The $\theta$-boundary conditions (37) give only $4N-3$ (complex) equations, while $5N-4$ boundary conditions are required at the initial time. The “left over” $N-1$ conditions correspond to $N-1$ gauge degrees of freedom in the model. As described in section III.2, for these $N-1$ equations one has to use the real part of Gauss’ law (26a),

[TABLE]

where $i=1,\ldots,N-1$. One also makes use of equations that fix the remaining real gauge freedom. The latter equations are

[TABLE]

for all $L=4N-2,\ldots,5N-4$. These modes have zero frequency $\omega_{L}=0$ and correspond to the unphysical degrees of freedom which change under gauge transformation, so Eq. (40) fixes the residual gauge invariance with real gauge functions completely. Gauge transformations with imaginary gauge functions are forbidden by the reality conditions at final time.

Final boundary conditions.

It is straightforward to implement the reality conditions at final time (9c). Supposing that the last two time grid points, $N_{t}$ and $N_{t}+1$ are on the real time axis, they are

[TABLE]

For energies below the bifurcation line $E_{1}(N)$ the time $N_{t}$ can be chosen to coincide with the point C of the time contour (so there are only two lattice points $N_{t}$ and $N_{t}+1$ on the whole CD part). For higher energies, though, the fields are not real along the most part of the real axis, so the part CD of the time contour has to be as long as possible (see Section V).

Fixing time translational invariance.

One more complication is that, in the continuous formulation, the boundary value problem (9) has an invariance under translations along real time (both field equations and boundary conditions are invariant under such a translation). To define properly the boundary value problem, one has to fix the position of the solution in time. In the lattice version this invariance is violated by the discretization and finite volume effects, but this violation does not enable one to control the position of the time contour relative to the branching points of the solution.

The existence of this invariance means that one of the equations is redundant (if discretization and finite volume effects are discarded). Somewhat arbitrarily, we take as redundant one of the real equations entering the $\theta$-boundary conditions (9b)

[TABLE]

for a specific mode. Provided the system linearizes at initial time, this equation is indeed a consequence of the others. The reason is that reality conditions at final time imply that the (conserved) energy is real. Hence the (linearized) energy (10) is real at initial time. Then one of the modes automatically obeys Eq. (42) provided all other modes obey the $\theta$-boundary condition (9b).

This suggests the following modification of the equations. One of the equations (9b) is changed to

[TABLE]

whose lattice version is (cf. Eq. (37))

[TABLE]

Thus, instead of Eq. (42) one imposes another boundary condition, which is not invariant under time translations. The choice of the latter is a matter of convenience. We control the position of the solution in time by imposing the boundary condition that fixes the “center-of-mass” of the field $\chi$ at the initial time to be equal to a given $R$:

[TABLE]

This prescription works if the mode $z_{K}$ in Eq. (43) is reasonably occupied at the initial time, otherwise the equation (42), which is “thrown away”, is nearly degenerate. Aside from this, the results of the calculations do not depend significantly on the mode chosen.

The relative phase between $f_{\bm{k}}$ and $g_{\bm{k}}$ can be used to check the validity of the calculations. In the linear regime it should be equal to zero, so the actual value of this phase indicates how close the system is to the linear regime at the initial time.

To summarize, the lattice boundary value problem consists of the field equations, obtained from action (32) for all inner lattice points ($i=1,\ldots,N-1$, $j=0,\ldots,N_{t}$, a total of $(N-1)(N_{t}+1)$ equations), the final reality boundary conditions (41) ($N-1$ equations), the $\theta$ boundary conditions (37) for all modes except one mode $z_{K}$ ($N-2$ equations), and a pair of real equations (43), (44) (one complex valued equation). This makes $(N-1)(N_{t}+3)$ complex equations for the same number of variables.

IV.4 Search for solutions

The equations to be solved make a set of discretized partial differential equations which change their signature from hyperbolic on the Minkowskian parts of the time contour to elliptic on the Euclidean part. The problem at hand is a boundary value problem which cannot be transformed into an initial value one. This means that the equations can be solved only globally, as a set of nonlinear equations at all $r,t$ grid coordinates.

To deal with the non-linear system of equations we employ a multidimensional analog of the Newton–Raphson method which approaches the desired solution iteratively. At each iteration, the linearized equations in the background of the current approximation are solved. The next approximation is obtained by adding the solution to the background, and the procedure is repeated. The advantage of the algorithm is that it does not require positive-definiteness of the matrix of second derivatives. It is, however, sensitive to zero modes. In the absence of zero modes, the algorithm converges quadratically; the accuracy of $10^{-9}$ is typically reached in 3-5 iterations. The convergence slows down in the presence of very soft modes, as typically happens near bifurcation points.

IV.5 Elimination algorithm

The discrete version of the equations derived from (32) is

[TABLE]

(here $\varphi_{jI}=\{a,\alpha,\beta,\mu,\nu\}(t_{j},r_{i})$ and $I$ runs from [math] to $5N-4$, $j=0\dots N_{t}$). The Newton–Raphson iteration is

[TABLE]

(all other second derivatives are zero) and before the next Newton–Raphson step the fields are changed according to

[TABLE]

Equations (45) can be rewritten in matrix form

[TABLE]

where $u_{j}$ and $b_{j}=\frac{\partial S}{\partial\varphi_{jI}}$ are $(5N-4)$–dimensional vectors, $\tilde{D}^{(-)}_{j}=\frac{\partial^{2}S}{\partial\varphi_{jI}\partial\varphi_{j-1,K}}$, $\tilde{D}_{j}=\frac{\partial^{2}S}{\partial\varphi_{jI}\partial\varphi_{jK}}$, $\tilde{D}^{(+)}_{j}=\frac{\partial^{2}S}{\partial\varphi_{jI}\partial\varphi_{j+1,K}}$ are $(5N-4)\times(5N-4)$ matrices. By multiplying Eq. (46) by $\tilde{D}^{-1}_{j}$ we get

[TABLE]

with $D^{(\pm)}_{j}=-\tilde{D}^{-1}_{j}\cdot\tilde{D}^{(\pm)}_{j}$, $b_{j}=-\tilde{D}_{j}^{-1}\cdot\tilde{b}_{j}$. This system of linear equations was solved by the following version of “divide–and–conquer” elimination algorithm. Excluding $u_{j}$ for some $j$ gives

[TABLE]

[TABLE]

Since the elimination of an equation changes only adjacent equations, it is possible to eliminate all equations with odd $j$ in parallel, and arrive to a system of the type (47) again, but with two times less variables and equations. This is the second level of elimination. After a series of eliminations we arrive at a system of only two equations for $j=0$ and $j=N_{t}$:

[TABLE]

where $\hat{D}^{(\pm)}$ and $\hat{b}$ have the values resulting from the elimination of all intermediate equations. Solving them together with the boundary conditions444Unlike the field equations at the intermediate points (45), which are analytic, the boundary conditions involve complex conjugation. So, equations (48) and boundary conditions are to be viewed as eight real matrix equations. All the elimination calculations (and reconstruction of field values afterwards) can be done, however, with complex algebra, which is two times more efficient. (also linearized), which involve $u_{-1}$ and $u_{0}$ for initial boundary condition and $u_{N_{t}}$ and $u_{N_{t}+1}$ for the final one, we determine the corrections $u_{-1}$, $u_{0}$, $u_{N_{t}}$ and $u_{N_{t}+1}$. Then it is straightforward to reconstruct $u$ at all intermediate points, using the equations (47) for each elimination level.

IV.6 Solutions below the sphaleron energy

The Newton–Raphson method requires a good initial approximation for the solution. This favors the following general strategy. We first find the periodic instanton solution, which corresponds to $\theta=0$ Bonini:2000mb and can be obtained via a minimization procedure. After the periodic instanton is found, we change the parameters $T$ and $\theta$ in small steps, using the solution from the previous step as a starting configuration. At each step we then calculate the energy $E$, number of particles $N$ and the suppression exponent $F(E,N)$ for the solution obtained.

This procedure is illustrated in Fig. 9, where each dot represents one solution of the boundary value problem. Initial periodic instanton configurations correspond to the points on the upper left line in the figure. Starting from these points, the value of $\theta$ was increased, and lines with constant values of $T$ were obtained until the bifurcation line was met. Data obtained in this way make almost straight lines in the left part of Fig. 9.

The boundary value problem (9)–(11) does not explicitly refer to the topological properties. Hence, it is not guaranteed that its every solution describes a transition between topologically distinct vacua. This is not a problem at $\theta=0$, because of the proper topological structure of the periodic instanton solutions. But at non-zero $\theta$ one should check that the solution indeed has correct topology.

The topological properties of a given solution are associated with the behavior of the phases of the fields, see Fig. 8. A very useful tool to control the properties of the solution is visualization of the field behavior. The visualization of a representative field configuration is presented in Fig. 10. It describes the field $\chi(r,t)$, with the phase of the field encoded in color. The Euclidean part of the time contour is inclined to make it distinct from the Minkowskian parts. In the initial state (left part of the surface) the field is close to its vacuum value, with excitation in the form of the incoming spherical wave moving towards to $r=0$. The final state (right part of the surface) contains the outgoing wave. The phase of the field clearly behaves differently in the initial and final states. This confirms that the topological transition indeed has occurred (compare to upper three images in Fig. 8). Several other important properties of the solution may also be seen immediately. These are: the moment when the field goes through $\chi=0$ in the middle of the Euclidean evolution, which of course should happen with the field evolving between neighboring vacua; the wide outgoing wave, suggesting that a large number of low energy particles is created after the transition; the small and relatively sharp incoming wave, meaning that higher energy modes are occupied and the number of particle in the incoming state is smaller.

V Going over the sphaleron energy

The procedure described above works as it is for relatively low energies $E\lesssim E_{\mathrm{sph}}$ only. With growing energy, the solutions on the CD part of the contour tend to stay for a long time close to the sphaleron. As the energy approaches some $N$-dependent value $E_{1}(N)$ this time tends to infinity, and if one continues to search for solutions to the boundary value problem (9) with reality condition imposed at finite positive time, the solutions above this energy have wrong topological properties, i.e. they end up in the same topological vacuum as the initial one (see Fig. 11). This situation is not specific to the $SU(2)$ gauge model studied here, but appears quite generally in quantum mechanical tunneling with multiple degrees of freedom. It was observed also in the study of the false vacuum decay in scalar field theory  Kuznetsov:1997az and in quantum mechanics with two degrees of freedom Bonini:1999cn ; Bonini:1999kj . The phenomenon was studied in detail in Ref. Bezrukov:2003yf in the case of quantum mechanics of two degree of freedom, and a general method of dealing with this difficulty was proposed there and checked against the exact solution of the Schrödinger equation. We describe here its gauge field version.

As suggested in Ref. Bezrukov:2003yf , the line $E_{1}(N)$ is the bifurcation line at which two types of solutions to the boundary value problem (9) meet. These are (i) solutions which end up close to the same vacuum as the initial one and (ii) solutions that arrive at the sphaleron with excited positive modes (in the case of field theory these excitations fly away quickly in the form of spherical waves in the sphaleron background). The former solutions are unphysical, while solutions of the latter form determine the tunneling exponent. For the interesting solutions of type (ii), the condition (9c) is satisfied only asymptotically, so it is very hard to find them numerically. A way out is to introduce a small regularization parameter into the equations of motion, which would not allow a solution to stay close to the sphaleron for infinite time. The final result is then obtained in the limit of zero regularization parameter.

To implement these ideas we start with the regularized expression for the cross section

[TABLE]

where $\epsilon$ is a small parameter and $T_{\mathrm{int}}$ is a functional proportional to the time the system spends in the interaction region. In case of gauge–Higgs theory we use the functional

[TABLE]

The path integral for (49) is no longer saturated by classical solutions spending infinite time close to the sphaleron, while the original cross section $\sigma(E,N)$ is obtained in the limit $\epsilon\to 0$. The boundary value problem for (49) coincides with the unregularized one, but with the action modified by adding an imaginary term of the form

[TABLE]

The equations of motion (26) are modified accordingly.

The functional (50) has several important features: i) it is gauge invariant; ii) it is large and positive on configurations close to the sphaleron (where the fields stay for a long time away from their vacuum values); iii) it does not change the free dynamics in the linear region, since it does not produce quadratic terms in the expansion of the action about vacuum (this is important for the boundary conditions to be unaffected).

With this regularization one obtains results for all energies $E$. The procedure is as follows. One introduces small but non-zero $\epsilon$ at energies below $E_{1}(N)$, then obtains solutions with proper topology for any energy. Then one takes the limit $\epsilon\to 0$. Upon taking this limit, the configurations with $E>E_{1}(N)$ stay for longer time close to the sphaleron, which means that in the limit of $\epsilon\to 0$ the solution tunnels “onto” the sphaleron.

Moreover, at the boundary of the classically allowed region, the solutions to the regularized problem merge smoothly with the classical over-barrier topology changing solutions, because the bifurcation on the boundary of the classically allowed domain is regularized exactly in the same way as the bifurcation at $E=E_{1}(N)$. At the boundary of the classically allowed region $F=0$ by definition, so the regularized version of this functional $F_{\epsilon}$ is proportional to $\epsilon$. This means that $T$ and $\theta$ are also proportional to $\epsilon$ there, and as the regularization is turned off, $\epsilon\to 0$, both $T$ and $\theta$ disappear, leading to purely real classical boundary value problem in real time.

There is one more complication in the $SU(2)$ field theory, which is relevant to this procedure. For energies $E<E_{1}(N)$ the amount of time $T$ spent on the Euclidean part of the contour is a growing function of energy (as opposed to the situation in two dimensional quantum mechanics Bezrukov:2003yf ), while it is zero at the boundary of the classically allowed region $E_{0}(N)$. This means that $T$ as function of energy at fixed $\theta$ has a maximum somewhere in between $E_{1}(N)$ and $E_{0}(N)$ (the calculations show that the maximum is actually at $E_{1}(N)$, see Fig. 13). For $T$ close to this maximum value, the Newton-Raphson method fails because of the presence of two nearby solutions with equal values of $T$. This new bifurcation is absent, if one searches for solutions with fixed energy $E$, instead of fixed $T$. To formulate the boundary value problem with fixed $E$ instead of fixed $T$, one simply sets the contour height $T^{\prime}$ to some conveniently chosen value, and leaves $T$ as a free variable. This leads to a trivial modification of the initial boundary conditions,

[TABLE]

An additional equation is then required to find one extra undetermined variable $T$. This is the equation involving energy of the solution, Eq. (10). With this modification of the procedure, the “bifurcation” corresponding to the maximum of $T$ disappears. Note however, that this method had to be applied with great care—on realistic grids it is hard to achieve fully linear regime in the initial state, therefore the difference between $T$ and $T^{\prime}$ must be small, as the dependence of the fields on imaginary time is exponential. In our calculations, the version of the boundary value problem with fixed $E$ ($T^{\prime}$ different from $T$) was used only to cross the maximum of $T$.

In our numerical calculations we introduced small non-zero $\epsilon$ as the energy of the solutions approached $E_{1}(N)$. Simultaneously, the modification (52) was used to get past the maximum of $T$ at $E_{1}(N)$. In Fig. 9 the solutions obtained with this modification are represented by points on the line that crosses the bifurcation line $E_{1}(N)$. At higher energies the modification (52) is no longer needed, and only regularization with non-zero $\epsilon$ was used. In Fig. 9 solutions to the regularized problem in region A.II correspond to points on the curved lines in the right part of the plot (lines of constant $T$). The line with the highest energy has zero suppression exponent and corresponds to the boundary of the classically allowed region. To connect lines of constant $T$, we obtained a set of solutions represented by the irregular line in right part of Fig. 9.

VI Numerical results

There are several factors affecting the choice of the lattice size and shape. The physical spatial size of the lattice $L$ is chosen large enough to make comfortable room for the sphaleron. More importantly, $L$ determines how close to the linear regime the system is in the initial state: the farther away from the origin, the smaller becomes the amplitude of incoming spherical waves. After $L$ is chosen, the lengths of the parts AB and CD of the time contour are determined completely; the length of the AB part $T_{AB}$ is slightly smaller than $L$, so that the incoming wave does not reach the spatial boundary $r=L$. The length of the CD part $T_{CD}$ is zero for energies below the bifurcation energy $E<E_{1}(N)$. For higher energies, $T_{CD}$ is adjusted to be long enough, so that the solution gets close to the vacuum configuration, and the regularization (51) does not contribute significantly to the equations of motion at the final moment of time when the reality boundary conditions are imposed.

The lattice spacing $\Delta r$ constrains the precision of the discretization in two different ways. First, it is chosen to be substantially smaller than the size of the instanton-like part of the configuration, i.e. the characteristic scale of nonlinear dynamics that occurs near $r=0$ during the topological transition itself. Second, $\Delta r$ controls the energy of the hardest mode in the initial state, thus limiting the lowest particle number $N$ that can be reached for given energy $E$. The time spacing $\Delta t$ is chosen to be smaller than $\Delta r$ to guarantee stability of the numerical procedure.

The amount of computer memory required for a lattice of spatial size $N_{r}$ and time length $N_{t}$ is approximately $2\times N_{t}(5N_{r})^{2}\times 16$ bytes (see Section IV.5), while the CPU time of one Newton–Raphson iteration scales roughly as $N_{t}(5N_{r})^{3}$. It was noted in Section IV.5 that the algorithm is suitable for parallel execution, so one divides this time by the number of processors available for the calculations555The parallelization algorithm is effective only if $N_{\mathrm{processors}}<\sqrt{N_{t}}$, so the shortest possible wall clock time in an ideal situation is proportional to $\sqrt{N_{t}}(5N_{r})^{3}$. Overall, the most strongly constrained is the spatial size $N_{r}$: a two times larger spatial grid means eight times longer processor time.

The main results we present in this paper were obtained on a grid with spatial size $L=8$ (i.e. $L=8/\sqrt{2}M_{W}$) and number of spatial grid points $N_{r}=90$. The length of the initial Minkowskian part of the contour $T_{AB}$ was equal to 6. The number of time grid points $N_{t}$ on the part AB of the contour was equal to 200, while on the Euclidean part BC it was equal to 150. The number of points on CD part varied from 2 for energies $E<E_{1}(N)$ to about 400 for higher energies (when the $\epsilon$–regularization was used). On the largest grids the amount of memory used was 4Gb, and it took 3 minutes for one Newton–Raphson iteration on a 16 processor IBM-RS/6000 supercomputer, or about 15 minutes for one full solution.

We obtained the results for the suppression factor in the region of $E$ and $N$ shown in Fig. 9. For the lattice parameters we used, this region is limited mainly by the effects of nonlinearity at the initial time, preventing us from reaching smaller particle numbers. When energy and particle number are small simultaneously (bottom-left part of the plot), effects of the spatial discretization (finite $\Delta r$) are also important.

To check the discretization effects, a limited set of calculations was performed on smaller grids. The results presented here coincide with results obtained with $N_{r}=64$ with precision better than 1% (except for very small energies). With $N_{r}=45$, on the other hand, the results coincide only for sufficiently large initial particle numbers, exactly as one would expect.

The linearization of the system in the initial state can be checked by evaluating the time dependence of the linear energy (10) and particle number (11) on the part AB of the contour. For linearized system, these should be independent of time. For a typical configuration this test is shown in Fig. 12. The linear energy coincides with the exact one in the initial state with precision of order of 1% or better, which confirms that the solution is quite close to the linear regime. Another test of linearity is the amount of the violation of the initial boundary condition (42), which is discarded to impose the time translation invariance fixing relation (Section IV). This amount grows towards smaller $N$, and apparently this is one of the effects preventing us from going to lower $N$ with the current spatial lattice size $L=8$. Larger lattices are needed to achieve better linearization on the initial part of the time contour and thus reach smaller particle numbers.

We made additional checks of the precision of the numerical calculations, including conservation of energy and the inverse Legendre transform

[TABLE]

These checks are satisfied with precision better than $10^{-3}$. This means that the precision of the final results is determined mostly by the quality of the linearization in the initial state (of the order of 1%).

Lines of constant $T$ and constant $\theta$ are shown in Figs. 13 and 14. One observes from Fig. 14 that $\theta$ grows as $N$ decreases, as expected, and $\theta$ is equal to zero on the periodic instanton line and on the boundary of the classically allowed region $E_{0}(N)$. The lines of constant $T$ show that $T$ also equals to zero at the boundary $E_{0}(N)$, and reaches a maximum for given $N$ (and for given $\theta$ also) at the bifurcation line $E_{1}(N)$. Close to this line we made use of the modification of the boundary value problem described at the end of Section V.

Representative solutions are shown in Fig. 15. They correspond to deep tunneling regime ($E<E_{1}(N)$), tunneling onto the sphaleron ($E>E_{1}(N)$) and classical overbarrier transition at $E_{0}(N)$, all for $N=1$. One can see from the color patterns that the field indeed undergoes the topology changing transition of the form illustrated in Fig. 8. The incoming wave is present in the left part of the pictures, becoming sharper and sharper for higher energy (the particle number is the same for all plots). In the first picture the topological transition is seen on the Euclidean part of the contour. In the second and third pictures, a sphaleron like configuration is visible on the right, with “extra” waves (excitations about the sphaleron) flying away, while the sphaleron itself starts to decay quite close to the right end of the plot (with the regularization parameter $\epsilon$ tending to zero, the moment of sphaleron decay moves towards larger times). At large times, the wave reflected from the boundary $r=L$ appears due to the Dirichlet boundary conditions (31) imposed at $r=L$. This wave does not alter the results, as it occurs in the linear regime666To get rid of the reflected wave a much larger spatial grid would be needed..

The fact that for $E>E_{1}(N)$ the solution after tunneling has the form of the sphaleron plus spherical excitations in its background is illustrated by plotting the spatial energy density at different times after tunneling. In Fig. 16 the energy density distribution is shown for the middle solution of Fig. 15. As the time increases, the bump on the right (spherical wave) moves towards larger $r$, while the energy density profile approaches that of the sphaleron.

It is also instructive to see that with the number of incoming particles decreasing, the occupied modes have higher frequencies. This is demonstrated in Fig. 17 for energy slightly smaller than the sphaleron energy.

Finally, let us discuss the extrapolation of the results to zero number of particles, which we performed to obtain predictions for the suppression exponent $F_{HG}(E)$ of the two-particle cross section.

Two ways of obtaining the lower bounds on $F_{HG}(E)$ were explained in Section II. One is to continue $F(N)$ at each energy linearly in $N$ to $N=0$ (this is justified by recalling that $\partial F/\partial N\propto-\theta$ increases as $N\to 0$), while the other is to continue lines of constant $F$ (see Fig. 3) linearly to $N=0$ (this gives a lower bound since the lines of constant $F$ have positive curvature). Both these extrapolations are straightforward to make, insofar as the required derivatives of $F(E,N)$ are given for each configuration by the values of $T$ and $\theta$ through the relations (53) and (54). In this way we obtained the bound shown in Fig. 5.

We now elaborate on our estimate of the function $F_{HG}(E)$ itself, Fig 6. Perturbative calculation at low energies Bezrukov:2003zn shows that while the exponent $F(E,N)$ has singular behavior of the form $N\log(N)$, the function $T(E,N)$ is regular in $N$ and close to a linear function. At large $N$ the numerical data also demonstrate that the behavior of $T$ is close to linear. This is shown in Fig. 18. For $E>E_{\mathrm{sph}}$ numerical results suggest that $T(N)$ is almost constant at small $N$. These properties justify a linear extrapolation of $T(N)$ to $N=0$ with energy kept constant. After obtaining $T(E)$ at zero particle number, $F_{HG}$ is readily found by integrating Eq. (54) starting from the instanton value $F_{HG}(E=0)=1$. The resulting estimate is presented in Fig. 6, solid line.

VII Conclusions

Our study shows that the semiclassical procedure coupled to suitable computational techniques is capable of producing quantitative results for semi-inclusive, weakly coupled non-perturbative processes, and, in particular, for particle collisions.

In this paper we applied this technique to study the suppression factor for topology changing transitions, and accompanying baryon and lepton number violation, in the SU(2) sector of the electroweak theory up to energies well above the sphaleron energy. We imposed spatial spherical symmetry, so our results are valid, strictly speaking, for s-wave scattering.

Our results show that the known analytic expression for the suppression exponent, which contains three terms of low-energy expansion, works well up to the sphaleron energy, but underestimates the suppression at higher energies.

By numerical analysis we have found that baryon and lepton number violation, accompanying topology changing s-wave particle collisions in the electroweak theory, remains highly suppressed up to energies of at least $\sim 250$ TeV (and likely much higher).

Acknowledgements.

The authors are indebted to A. Kuznetsov for helpful discussions. We wish to thank Boston University’s Center for Computational Science and Office of Information Technology for generous allocations of supercomputer time. Part of this work was completed during visits by C.R. at the Institute for Nuclear Research of the Russian Academy of Sciences and by F.B., D.L., V.R. and P.T. at Boston University, and we all gratefully acknowledge the warm hospitality extended to us by the hosting institutions. The work of PT was supported in part by the SNSF grant 21-58947.99. This research was supported by Russian Foundation for Basic Research grant 02-02-17398, U.S. Civilian Research and Development Foundation for Independent States of FSU (CRDF) award RP1-2364-MO-02, and DOE grant US DE-FG02-91ER40676.