text:
# Drag Reduction by Polymers in Turbulent Channel Flows: Energy Redistribution Between Invariant Empirical Modes.

## I Introduction

“Drag reduction” refers to the interesting observation that the addition of a few tens of parts per million (by weight) of long-chain polymers to turbulent fluids can bring about a reduction of the friction drag by up to 80% 00SW . Obviously, the phenomenon has far reaching practical implications besides being challenging from the fundamental point of view. In spite of intense interest for an extended period of time 69Lu ; 75Virk ; 90Ge , Sreenivasan and White 00SW  recently concluded that “it is fair to say that the extensive - and continuing - activity has not produced a firm grasp of the mechanisms of drag reduction”. In this paper we want to advance on the basis of recent Direct Numerical Simulation of model viscoelastic hydrodynamic equations 97THKN ; 98DSB ; 00ACP ; 02ACP . Such DNS show unequivocally that drag reduction is reproduced by model equations like the FENE-P model. From the theoretical point of view this is significant, since it indicates that the phenomenon is included in the solutions of the model equations. Understanding drag reduction then becomes a usual challenge of theoretical physics.

Our thinking here is motivated in part by a recent analysis of the stability of laminar channel flows subject to space dependent effective viscosity 01GLP ; 02GLP . It turned out the even small viscosity gradients can lead to a giant stabilization of the most unstable modes, both for primary and secondary instabilities. In these cases one can understand the phenomenon completely by examining the energy budget of the putative unstable modes and their interaction with the mean flow; the most important observation had been that it is the the existence of viscosity gradients positioned at a strategic distance from the wall which is crucial for the existence of a large effect. It seems desirable to do something similar for the viscoelastic turbulent flows as well (in which the space dependent effective viscosity arises due to differential stretching of the polymers). But alas, in distinction from primary and secondary instabilities, where it is obvious which are the relevant modes, for the turbulent flow these are not known apriori. We therefore decided to first initiate a systematic study of the empirical modes that are sustained in the turbulent flow, and then to discuss their interaction with the mean flow and with the polymeric additive, their stabilization or destabilization when we compare the viscoelastic to the Newtonian flow, and their energy budget. In this paper we present the first results of this study.

We will demonstrate that we can determine with reasonable accuracy at least the first thirty most energetic modes that are sustained in the turbulent flow, for both the FENE-P and the Navier-Stokes equations (run at the same friction Reynolds number, and see below for details). These modes can be arranged in descending order according to their relative energy. Unexpectedly we find that the nature of the most relevant modes is unchanged in the two cases. On the other hand the energies associated with the modes and their relative ordering are changed; some modes that are energetic in one flow are strongly suppressed (their energy decreases by a factor of 4) in the other flow, and vice versa. Most importantly, the few most energetic modes of the viscoelastic flow contain a lot more energy that the same number of most energetic Newtonian modes. We propose therefore that drag reduction should be understood by examining the dynamics and relative stability of the energy containing modes rather than focusing only on the dissipative end of the spectrum, as proposed for example in 90Ge .

In Sects. II A and B we summarize the FENE-P equations and the numerical approach. In Sect. II.3 we present the essential results regarding the observation of drag reduction. In Sect. III we review the Karhunen-Loéve method for determining the best empirical modes, and apply it to the problem at hand. In Sect. IV we discuss the results, demonstrate the invariance of the modes, and present the relative ordering. Sect. V is devoted to a discussion of the findings and of the road ahead.

## II Equations of Motion and Direct Numerical Simulations

### II.1 The FENE-P model for dilute polymers

The addition of a dilute polymer to a Newtonian fluid gives rise to an extra stress tensor $`𝓣(𝒓,t)`$ which affects the Navier-Stokes equations:

$`{\displaystyle \frac{𝒖}{t}}+(𝒖\mathbf{})𝒖`$ $`=`$ $`\mathbf{}p+\nu _\mathrm{s}^2𝒖+\mathbf{}𝓣,`$
$`\mathbf{}𝒖`$ $`=`$ $`0.`$ (1)
Here $`𝒖(𝒓,t)`$ is the solenoidal velocity field, $`p(𝒓,t)`$ is the pressure and $`\nu _s`$ is the viscosity of the neat fluid. The additional stress tensor $`𝓣`$ is not known exactly, and needs to be modeled. There are a number of competing derivations, see 87BCAH ; 94BE . In our work we adopt the FENE-P model that is derived on the basis of a dumbbell model for the polymers 87BCAH ; 94BE  and is known to reproduce the phenomenon of drag reduction. In this model $`𝓣`$ is determined by the “polymer conformation tensor” $`𝑹`$ according to

$$𝓣(𝒓,t)=\frac{\nu _\mathrm{p}}{\tau _\mathrm{p}}\left[\frac{f(𝒓,t)}{\rho _0^2}𝑹(𝒓,t)\mathrm{𝟏}\right].$$
(2)
Here $`\mathrm{𝟏}`$ is the unit tensor, $`\nu _\mathrm{p}`$ is a viscosity parameter, $`\tau _\mathrm{p}`$ is a relaxation time for the polymer conformation tensor and $`\rho _0`$ is a parameter which in the derivation of the model stands for the rms extension of the polymers in equilibrium. The function $`f(𝒓,t)`$ limits the growth of the trace of $`𝑹`$ to a maximum value $`\rho _\mathrm{m}`$:

$$f(𝒓,t)\frac{\rho _\mathrm{m}^2\rho _0^2}{\rho _\mathrm{m}^2R_{\gamma \gamma }(𝒓,t)}.$$
(3)
The model is closed by the equation of motion for the conformation tensor which reads

$`{\displaystyle \frac{R_{\alpha \beta }}{t}}`$ $`+`$ $`(𝒖\mathbf{})R_{\alpha \beta }={\displaystyle \frac{u_\alpha }{r_\gamma }}R_{\gamma \beta }+R_{\alpha \gamma }{\displaystyle \frac{u_\gamma }{r_\beta }}`$ (4)
$``$ $`{\displaystyle \frac{1}{\tau _\mathrm{p}}}\left[f(𝒓,t)R_{\alpha \beta }\rho _0^2\delta _{\alpha \beta }\right].`$
The model is derived by assuming that the polymer can be characterized completely by an end-to-end vector distance. Nevertheless the resulting equations could be written on the basis of plausible arguments including up to quadratic terms in gradients and available tensors. For our purposes the accuracy of the model in reproducing quantitatively all the phenomena of turbulence in viscoelastic fluids is not an issue of concern. We are mainly interested in the fact that these equations were simulated on the computer in a channel geometry and exhibited the phenomenon of drag reduction as discussed below. Our aim is to understand drag reduction within the FENE-P model.

### II.2 Direct Numerical Simulations

A simple flow geometry which exhibits the phenomenon of drag reduction is channel flow between two parallel planes, separated by the distance $`2H`$ in the $`y`$-direction, see Fig. 1. The computational domain is periodic in the two directions parallel to the wall (streamwise $`x`$, spanwise $`z`$).

The Navier-Stokes equation (1) is written in terms of wall-normal component of velocity $`u_y(𝒓,t)`$ and the vorticity $`\omega _y(𝒓,t)`$, where $`𝝎=\mathbf{}\times 𝒖`$,

$`{\displaystyle \frac{u_y}{t}}+(𝒖\mathbf{})u_y`$ $`=`$ $`{\displaystyle \frac{p}{y}}+\nu _s^2u_y+_\beta 𝒯_{\beta y}`$
$`{\displaystyle \frac{\omega _y}{t}}+\left(𝒖\mathbf{}\right)\omega _y`$ $`=`$ $`\left(𝝎\mathbf{}\right)u_y+\nu _\mathrm{s}^2\omega _y`$ (5)
$`+ϵ_{y\beta \gamma }_\beta _\delta 𝒯_{\delta \gamma }`$
with $`ϵ_{\alpha \beta \gamma }`$ the fully antisymmetric Levi-Civita tensor. Given $`u_y`$ and $`\omega _y`$, the components of the velocity field in the two directions parallel to the $`xz`$ plane follow from the continuity equation and the definition of $`\omega _y`$,

$`\mathbf{}_{}𝒖_{}={\displaystyle \frac{u_y}{y}}`$ $`,`$ $`\mathbf{}_{}\times 𝒖_{}=\omega _y,`$ (6)
where the subscript $``$ denotes projection of vectors on the $`xz`$ plane. The system of Eqs. (II.2,6) is one of the standard formulations used for the direct numerical simulation of turbulent channel flows, since the procedure yields a solenoidal velocity field to machine accuracy.

In light of the periodicity in the $`xz`$ plane, it is natural to expand the planar components of the velocity field in Fourier modes. For the wall-normal direction one uses Chebyshev polynomials. The time stepping is carried out by a mixed Crank-Nicolson/Runge-Kutta scheme for the viscous and the nonlinear terms, respectively. The integration in the normal direction is done by Chebyshev tau-method and a standard de-aliasing technique is adopted for the nonlinear terms. The typical simulations have been performed on a computational grid of $`96\times 129\times 96`$ nodes in a domain of dimensions $`2\pi H\times 2H\times 1.2\pi H`$. Turbulence is maintained by enforcing the same constant pressure gradient for the corresponding Newtonian and viscoelastic simulations.

In discussing the simulations, one notes that the channel half-height $`H`$ is only one of the parameters which sets up the external lengthscale. An additional important control parameter is the enforced pressure gradient at the wall, which determines the friction. Denoting by pointed brackets an average over time the friction parameter is defined as

$$\overline{\tau }_\mathrm{w}=H\frac{p}{x}.$$
(7)
This is the basic control parameter of the flow. By considering also the overall kinematic viscosity, $`\nu _\mathrm{f}=\nu _\mathrm{s}+\nu _\mathrm{p}`$, the traditional friction Reynolds number is defined as

$$\mathrm{Re}_\tau =\frac{u_\tau H}{\nu _\mathrm{f}},u_\tau =\sqrt{\overline{\tau }_\mathrm{w}},$$
(8)
where $`u_\tau `$ is the friction velocity. As customary in wall turbulence, one may introduce the inner, or viscous, length scale

$$\mathrm{}=\frac{\nu _\mathrm{f}}{u_\tau }.$$
(9)
The inner velocity scale $`u_\tau `$ has its counterpart bulk velocity

$$U_0=\frac{1}{2H}𝑑yu_x.$$
(10)
Correspondingly we have the outer and inner time scales $`T_0=H/U_0`$ and $`\nu _\mathrm{f}/u_\tau ^2`$ respectively.

In the sequel all quantities are made dimensionless unless stated otherwise. Those made dimensionless with respect to the appropriate inner scale will be denoted by the superscript $`+`$; no special symbol is used to denote normalization by an outer scale.

Finally there are parameters associated with the polymer. Foremost is the Deborah number which is defined as the ratio of the relaxation time $`\tau _\mathrm{p}`$ and a typical time scale of the fluid motion, i.e.

$`\mathrm{De}={\displaystyle \frac{\tau _\mathrm{p}}{T_0}}.`$ (11)
The parameter $`r_\mathrm{p}=\nu _\mathrm{p}/\nu _\mathrm{s}`$ measures the relative viscosity of the polymers with respect to the Newtonian solvent, and the non-linear characteristics of the spring is defined in terms of the ratio $`\rho _\mathrm{m}^2/\rho _0^2`$.

In all our comparisons of Newtonian and viscoelastic simulations, the friction Reynolds number is kept constant, at a typical value $`\mathrm{Re}_\tau =125`$. The correspondence between the two flows is obtained by fixing the computational domain, and choosing for the Newtonian case a viscosity equal to the overall viscosity of the solution. The parameters chosen for all the viscoelastic simulations are $`\mathrm{De}=25`$, $`\eta _\mathrm{p}=0.1`$ and $`\rho _\mathrm{m}^2/\rho _0^2=1000`$.

### II.3 Overview of drag reduction

In this subsection we review briefly the main results of the present simulations which demonstrate the phenomenon of drag reduction. For more detailed description, see 02ACP . The analysis presented in the following sections employs the very same DNS.

In comparing the viscoelastic and the Newtonian flows we maintain the friction Reynolds number (8) fixed. We reiterate that to achieve this we need to choose the viscosity of the Newtonian flow properly, since the viscoelastic wall stress contains a small viscoelastic contribution,

$$\overline{\tau }_\mathrm{w}=\nu _\mathrm{s}\frac{d\overline{U}}{dy}+\overline{𝒯}_{yx}.$$
(12)
The component of the extra-stress $`𝒯_{yx}`$ does not vanish on the average, contributing in our simulation about $`10\%`$ of the total drag.

The main observations regarding drag reductions are as follows:

(i) For a fixed mean pressure gradient at the wall, the viscoelastic flow exhibits an increased flow rate through the channel, see Fig. 2.

Shown there are the mean profiles in the streamwise direction, the other means vanish by symmetry:

$`\overline{U}_x(y)`$ $`=`$ $`\overline{U}(y)u_x(𝒓,t),`$ (13)
$`\overline{U}_y(y)`$ $`=`$ $`\overline{U}_z(y)=0.`$
The increase in the throughput entails an increased bulk Reynolds number,

$$\mathrm{Re}_\mathrm{b}=\frac{U_0H}{\nu _\mathrm{f}}.$$
(14)
(ii) Some typical length scales increase in the viscoelastic flow. The data in Fig. 2 can be recast in terms of the log-law of the wall,

$$\overline{U}^+=\frac{1}{k}\mathrm{log}y^++A,$$
(15)
where $`k0.4`$ is the Von Karman constant. The constant $`A`$ depends on the thickness of the buffer plus viscous sub-layers, defined as the distance between the wall and the beginning of the log-region. The log-law exists equally well for the viscoelastic as for the Newtonian flow with the same $`k`$, but the numerical value of the constant $`A`$ is substantially increased in the former. This thickening of the buffer layer is directly related to the increased flow rate. The effect of thickening buffer layer can be also measured by the increase of the span-wise scale of the streamwise velocity fluctuations. Decomposing the velocity field into a mean and a fluctuating part

$$u_\alpha (𝒓,t)=\overline{U}_\alpha (y)+\stackrel{~}{u}_\alpha (𝒓,t),$$
(16)
one introduces the correlation tensor

$$K_{\alpha \beta }(𝒓,𝒓^{})\stackrel{~}{u}_\alpha (𝒓,t)\stackrel{~}{u}_\beta (𝒓^{},t).$$
(17)
In Fig. 3 we show $`K_{xx}(𝒓,𝒓^{})`$ for $`𝒓=(x,y,z)`$, $`𝒓^{}=(x,y,z+Z)`$. \[Note that by homogeneity $`K_{xx}=K_{xx}(y,Z)`$\]. Fig. 3 demonstrated the increase in the span-wise correlation length of the streamwise velocity fluctuations which can be defined as

$$L_{xz}(y)\frac{1}{K_{xx}(y,0)}𝑑ZK_{xx}(y,Z).$$
(18)
(iii) Comparison of the time signals between the viscoelastic and Newtonian cases shows an alteration of the characteristic frequencies for the viscoelastic model. This corresponds quite well to the experimental data, for example of Luchik et al. Tiederman , in which a decrease of the bursting frequency is observed in drag reducing solutions. We do not reproduce these results here.

(iv) Finally, the root mean square fluctuations change significantly as a function of $`y`$. Denoting

$$U_\alpha (y)\sqrt{|\stackrel{~}{u}_\alpha (𝒓,t)|^2},$$
(19)
we display in Fig. 4 the $`y`$ dependence of the three components $`\alpha =x,y,z`$. The streamwise fluctuations are shown to increase with respect to the corresponding Newtonian flow, while both the span-wise and the wall-normal fluctuations decrease, in qualitative agreement with available experimental results Tiederman .

In conclusion, the FENE-P model is shown to exhibit the phenomenon of drag reduction in close similarity to experimental observations. In addition, DNS of this model provides complete information about the velocity field and the covariance tensor field $`𝑹(𝒓,t)`$ as a function of space and time. We thus feel confident that a sufficiently savvy analysis of this model and its turbulent flow pattern should provide insight into the mechanism of drag reduction.

## III Analysis Using Empirical Modes

In this section we present the analysis of the difference between viscoelastic and Newtonian flows in term of the empirical modes that are sustained in the turbulent flow. In choosing a method to extract the modes we are led by the desire to find the “best” modes, and “best” in our case will mean those that are most energetic. In fluid mechanics the energy is quadratic in the field, so “best” is related to closeness in the $`L^2`$ norm. As is well known, the standard method to find best representation in $`L^2`$ norms is the Karhunen-Loéve method. The aim of the Karhunen Loéve method is to provide a set of modes that optimally decompose the field of interest, in our case the velocity field averaged over time. The approach is guaranteed to yield the best set of truncated modes, meaning that the field cannot be approximated better (in $`L^2`$ norm) by any other set of the same number of modes. The method had been applied to fluid mechanics in a number of contexts, see for example Siro1 ; Siro2 ; Holmes  for details and relevant references. We first adapt the Karhunen-Loéve method to the present context, and then discuss the results of the analysis.

The core object for the analysis is the simultaneous correlation function of the velocity fluctuations which was introduced already in Eq.(17). Due to stationarity this object is time independent, and due to homogeneity in the $`x`$ and $`z`$ direction, we can write

$$K_{\alpha \beta }=K_{\alpha \beta }(x^{}x,z^{}z;y,y^{}).$$
(20)
In the translationally invariant directions $`xz`$ one cannot do better than Fourier decomposition. Accordingly we consider the partly decomposed object $`\widehat{Q}(𝒒|y,y^{})`$ defined as

$`\widehat{Q}(𝒒|y,y^{})`$ $``$ $`{\displaystyle \frac{1}{N_{}}}{\displaystyle \underset{𝒓_{}}{}}K_{\alpha \beta }(x^{}x,z^{}z;y,y^{})`$
$`\times \mathrm{exp}\{ı[q_x(x^{}x)+q_z(z^{}z)]\}.`$
We denoted the discrete two dimensional wave-vector in the $`xz`$ plane as $`𝒒=(q_x,q_z)=(\stackrel{~}{q}_x2\pi /\mathrm{\Lambda }_x,\stackrel{~}{q}_z2\pi /\mathrm{\Lambda }_z)`$ where $`\stackrel{~}{q}_x,\stackrel{~}{q}_z`$ are integers. $`N_{}=N_xN_z`$ and the sum is taken over the discrete set of $`xz`$ points with $`𝒓_{}=(j_x\mathrm{\Lambda }_x/N_x,j_z\mathrm{\Lambda }_z/N_z)`$ the $`xz`$ projection of $`𝒓`$.

The non-trivial empirical modes are obtained from the remaining dependence on $`y,y^{}`$, for each given planar $`𝒒`$. The Karhunen-Loéve method consists of finding the eigenfunctions $`\mathrm{\Psi }_{\alpha \beta }(𝒒,p|y)`$ which solve the eigenvalue equation

$$_1^1\widehat{Q}_{\alpha \beta }(𝒒|y,y^{})\mathrm{\Psi }_\beta (𝒒,\mathrm{p}|y^{})𝑑y^{}=E(𝒒,\mathrm{p})\mathrm{\Psi }_\alpha (𝒒,\mathrm{p}|y).$$
(22)
Here $`E(𝒒,\mathrm{p})`$ is the energy associated with the mode $`(𝒒,\mathrm{p})`$, ordered in decreasing energy order, with a given $`𝒒`$ wave-vector in the wall parallel directions, and p mode number in the $`y`$-direction, orthogonal to the plane walls. Denoting by $``$ the total energy in the system, we have

$$d^3𝒓\frac{\stackrel{~}{u}^2}{2}=\underset{𝒒,\mathrm{p}}{}E(𝒒,\mathrm{p}).$$
(23)
The modes labeled by $`(𝒒,\mathrm{p})`$ can be relabeled in order of decreasing energy by introducing an index $`n`$. In this notation the mode whose associated energy is $`E_n=E(𝒒,\mathrm{p})`$ are denoted by

$$\mathrm{\Psi }_\alpha (n|y)\mathrm{\Psi }_\alpha (𝒒_n,\mathrm{p}_\mathrm{n}|y),$$
(24)
where $`(𝒒_n,\mathrm{p}_n)`$ is the label of the $`n^{\mathrm{th}}`$ mode. It is useful to introduce also a set of fields indexed by $`n`$,

$$\mathrm{\Phi }_\alpha (n|𝒓)=\mathrm{\Psi }_\alpha (n|y)\mathrm{exp}\left[ı(𝒒_n𝒓_{})\right].$$
(25)
These function are orthogonal in the sense

$$\frac{1}{2N_{}}\underset{𝒓_{}}{}_1^1𝑑y\mathrm{\Phi }_\alpha (n|y,𝒓_{})\mathrm{\Phi }_\alpha ^{}(n^{}|y,𝒓_{})=\delta _{nn^{}}.$$
(26)
The fluctuation velocity field can be now expanded in terms of these fields,

$$\stackrel{~}{u}_\alpha (t,𝒓)=\underset{n}{}a_n(t)\mathrm{\Phi }_\alpha (n|𝒓).$$
(27)
The main advantage of the Karhunen-Loéve method is that any finite truncation of this expansion can be shown to yield a best approximation for the velocity field (in the $`L^2`$ norm). Moreover, due to the symmetry of the correlation matrix the modes are orthonormal (after suitable normalization), and correlation functions in the basis of these modes are diagonal:

$$a_na_n^{}^{}=E_n\delta _{nn^{}}.$$
(28)
Similarly one can have an optimal representation of the correlation matrix itself as

$$K_{\alpha \beta }=\underset{n}{}E_n\mathrm{\Psi }_\alpha (n|y)\mathrm{\Psi }_\beta (n|y^{})\mathrm{exp}\left[ı𝒒_n(𝒓_{}𝒓_{}^{})\right].$$
(29)
In analyzing the DNS the modes are determined by using an expansion in terms of Chebyshev polynomials $`T_k(y)`$ in the $`y`$ direction. Denote $`y_j`$ the j<sup>th</sup> node in the list of $`N_y`$ Chebyshev nodes,

$$y_j=\mathrm{cos}\left(j\frac{\pi }{N_y}\right)j=1,\mathrm{},N_y,$$
(30)
and express

$$\mathrm{\Psi }_\alpha (𝒒,\mathrm{p}|y)=\underset{k=0}{\overset{N_y1}{}}\widehat{\mathrm{\Psi }}_\alpha (𝒒,\mathrm{p}|k)T_k(y).$$
(31)
Then Eq. (22) is discretized as

$`{\displaystyle \underset{k=0}{\overset{N_y1}{}}}`$ $`A_{\alpha \beta }(𝒒|j,k)\widehat{\mathrm{\Psi }}_\beta (𝒒,\mathrm{p}|k)`$ (32)
$`=E(𝒒,\mathrm{p}){\displaystyle \underset{k^{}=0}{\overset{N_y1}{}}}\widehat{\mathrm{\Psi }}_\alpha (𝒒,\mathrm{p}|k^{})T_k^{}(y_j).`$
where

$$A_{\alpha \beta }(𝒒|j,k)\underset{j^{}=1}{\overset{N_y}{}}\widehat{Q}_{\alpha \beta }(𝒒|y_j,y_j^{})T_k(y_j^{})w_j^{}.$$
(33)
with $`w_j`$ the integration weights.

For a given $`𝒒`$, Eq. (32) results in a linear algebraic eigenvalue problem of order $`(N_y\times N_y)`$, with the p<sup>th</sup> eigenvector given by $`\widehat{\mathrm{\Psi }}_\alpha (𝒒,\mathrm{p}|k)`$, $`k=0,\mathrm{},N_y`$ and the corresponding eigenvalue $`E(𝒒,\mathrm{p})`$. Since we need to solve for every discrete wave-vector $`𝒒`$, there exist in total $`N_xN_z`$ eigenvalue problems to be solved in order to determine the full set of eigenfunctions. From the discrete eigenvectors, the corresponding discrete modes are constructed according to the discrete analog of Eq. (25),

$$\mathrm{\Phi }_\alpha (𝒒,\mathrm{p}|𝒓_{},y_j)=\mathrm{\Psi }_\alpha (𝒒,\mathrm{p}|y_j)\mathrm{exp}\left[ı(𝒒_n𝒓_{})\right],$$
(34)
where Eq. (31) is used to evaluate $`\mathrm{\Psi }_\alpha `$ at the Chebyshev nodes. The discrete modes are normalized according to the discrete version of Eq. (26),

$$\frac{1}{2N_{}}\underset{j=1}{\overset{N_y}{}}\underset{𝒓_{}}{}\mathrm{\Phi }_\alpha (𝒒,\mathrm{p}|𝒓_{},y_j)\mathrm{\Phi }_\alpha ^{}(𝒒,\mathrm{p}|𝒓_{},y_j)w_j=1.$$
(35)
In practice we ran DNS with and without polymers for $`5000`$ large-eddy turnover times $`T_0`$, see Sec. II.2, in statistically stationary conditions. The time step used for time-advancement, in terms of viscous units $`\nu _\mathrm{f}/u_\tau ^2`$, Sec. II.2, is $`Dt^+=0.05`$. Since the same pressure gradient is enforced in the two cases, the friction velocity is identical, while the bulk velocity increases by $`24\%`$ in the viscoelastic simulation. This corresponds to a reduction of the large-eddy turnover time, $`T_0`$.

We have collected $`N_{\mathrm{fields}}`$ fields, $`N_{\mathrm{fields}}100`$, displaced in time by $`50T_0`$. These fields were used to construct the discrete correlation function, Eq. (17), at the nodes of the computational grid. Eigenvalues and eigenvectors of its discrete Fourier transform in the wall parallel directions has then been evaluated, according to Eq. (32). The modes corresponding to the computed eigenvectors have been arranged in order of decreasing eigenvalue. From Eq. (35), the eigenvalues correspond to the average energies in the considered modes, so that the chosen arrangement is in fact an ordering in terms of decreasing energy content. When needed, the index of the energetic ordering of the Newtonian and viscoelastic simulations will be denoted by $`n__\mathrm{N}`$ and $`n_{_{\mathrm{VE}}}`$ respectively, dropping the subscript when $`n`$ refers to both cases.

## IV Analysis and results

### IV.1 The dominant modes are approximately invariant

The first discovery in the study of the empirical modes was admittedly surprising for the present authors, and in hindsight very serendipitous for the discussion of drag reduction. We found that the dominant modes are approximately invariant. In other words, the same modes that carry a sizable fraction of the energy of the viscoelastic flow appear essentially unchanged in the Newtonian flow, having practically the same spatial $`y`$-dependence. We will first demonstrate this surprising finding, and later focus on the difference between the two flows.

To discuss meaningfully the correspondence between the empirical modes in the two cases we seek a criterion of matching the viscoelastic modes to the corresponding Newtonian modes. This can be done easily, since both sets of modes form a complete orthonormal basis for solenoidal fields in the same domain. Each viscoelastic mode can then be expanded in terms of the Newtonian set,

$$\mathrm{\Psi }_\alpha ^{\mathrm{VE}}(𝒒,\mathrm{p}_{_{\mathrm{VE}}}|y)=\underset{\mathrm{p}__\mathrm{N}}{}A(𝒒,\mathrm{p}_{_{\mathrm{VE}}}|\mathrm{p}_\mathrm{N})\mathrm{\Psi }_\alpha ^\mathrm{N}(𝒒,\mathrm{p}__\mathrm{N}|y).$$
(36)
The complex amplitude $`A(𝒒,\mathrm{p}_{_{\mathrm{VE}}}|\mathrm{p}__\mathrm{N})`$ is given by the projection of the viscoelastic mode $`\mathrm{p}_{_{\mathrm{VE}}}`$ on Newtonian mode $`\mathrm{p}__\mathrm{N}`$, with identical wall-parallel wave-vector, $`𝒒`$:

$`A(𝒒,\mathrm{p}_{_{\mathrm{VE}}}|\mathrm{p}__\mathrm{N})=`$ (37)
$`{\displaystyle \frac{1}{2N_y}}{\displaystyle \underset{j=1}{\overset{N_y}{}}}\mathrm{\Psi }_\alpha ^{\mathrm{VE}}(𝒒,\mathrm{p}_{_{\mathrm{VE}}}|y_j)\mathrm{\Psi }_\alpha ^\mathrm{N}(𝒒,\mathrm{p}__\mathrm{N}|y_j)w_j.`$
We find that all the most energetic modes of the viscoelastic flow have one amplitude whose magnitude is close to unity. For example we plot the absolute magnitude $`|A(𝒒,\mathrm{p}_{_{\mathrm{VE}}}|\mathrm{p}__\mathrm{N})|`$ as a function of $`\mathrm{p}__\mathrm{N}`$ in Fig. 5 for the first three most energetic viscoelastic modes. Each of these modes displays an amplitude maximum very close to unity, implying that it matches well a single Newtonian mode.

This procedure furnishes a correspondence between viscoelastic and Newtonian modes: a given viscoelastic mode corresponds to the Newtonian mode for which the amplitude $`|A(𝒒,\mathrm{p}_{_{\mathrm{VE}}}|\mathrm{p}__\mathrm{N})|`$ is maximal. This unambiguously associates a single Newtonian mode to each viscoelastic mode. The value of the maximal amplitude, hereafter called matching parameter, gives a quantitative estimate of the difference in shape between two corresponding modes: A matching parameter equal to one implies absolute identity between the two modes in the energy norm. The correspondence in the spatial structure can then be verified by direct inspection. For instance, matching modes for the cases discussed in Fig. 5, where the matching parameter is well above $`0.9`$, are almost indistinguishable. Physically, the matching parameter represents the fraction of the energy in the viscoelastic mode that is ascribed to the matching Newtonian mode.

Fig. 6 shows the matching parameter for the first thirty most energetic viscoelastic modes ($`n_{_{\mathrm{VE}}}30`$). All the VE-modes except the 23rd and 29th have a matching parameter above $`0.9`$. We conclude that at least as far as the most energetic modes are concerned, the spatial structure of the modes is almost not altered by the polymer. Loosely speaking this can be expressed by saying that the modes are approximately invariant, i.e. fixed in shape. For practical purposes the difference between the Newtonian and viscoelastic empirical modes can be safely disregarded.

Next we discuss what changes from Newtonian to viscoelastic turbulence.

### IV.2 The energy and the relative ordering change

In light of the first discovery the second may be already anticipated: although the dominant modes hardly change, their energies and relative energy ordering change very significantly. We propose that understanding the change of energies and relative ordering of the modes will take us a long way in understanding drag reduction. We begin by comparing the most dominant modes in the two respective flows. Each empirical mode is identified by three numbers, $`(\stackrel{~}{q}_x,\stackrel{~}{q}_z)`$, and p, where as explained above, $`\mathrm{p}`$ corresponds to the index of the energetic ordering for fixed wave-vector $`𝒒`$.

#### IV.2.1 The most dominant mode

In Fig. 7 we present pictorially the mode (0,3,1) which is the most energetic for both the Newtonian and viscoelastic simulation. Its amplitude $`a_1(t)`$,

$$a_1(t)=\frac{1}{2N_y}\underset{j=1}{\overset{N_y}{}}\underset{𝒓_{}}{}\stackrel{~}{𝒖}_\alpha (t,y_j,𝒓_{})\mathrm{\Phi }_\alpha ^{}(n=1|y_j,𝒓_{})w_j,$$
(38)
is real in both flows, implying that the mode does not propagate, neither in the $`x`$ nor in the $`z`$ direction. Moreover, for this particular mode $`q_x=0`$, i.e. the field $`\mathrm{\Phi }_\alpha (n=1|𝒓)`$ is constant in $`x`$. It belongs to the class of non-propagating roll-modes discussed in Siro2 . The figure presents an iso-surface of the streamwise component of the velocity field associated with the mode, $`\mathrm{\Phi }_x(n=1|𝒓)=\mathrm{const}`$. The span-wise wavelength of the mode is $`\lambda _z=\mathrm{\Lambda }_z/\stackrel{~}{q}_z=2\pi H/3`$, and it appears to be confined relatively close to the channel walls.

A further pictorial presentation of the Newtonian mode (0,3,1) is shown in Fig. 8, where the real part of the mode (top panel) and imaginary part (bottom panel) are plotted separately. From the plots it is apparent that the mode is more or less localized in the vicinity of the walls, as already commented. The phase difference between the various components is also worth mentioning. Comparing the top and the middle panels, the stream-wise and span-wise components, $`\mathrm{\Psi }_x`$ and $`\mathrm{\Psi }_z`$ are real for all $`y`$, while the wall-normal component, $`\mathrm{\Psi }_y`$, has a constant phase lag of $`\pi /2`$ with respect to the other two, i.e. it is purely imaginary.

We note that although the first, most dominant mode, is the same for the two flows, the actual energy associated with this mode is twice larger in the viscoelastic mode. This is typical for all the leading modes in the viscoelastic flow as compared with the Newtonian flow; this central point of distinction between the two flows will be addressed in the next sub-section. The agreement between the leading modes ends with the first one; the majority of higher modes change in relative position in the energy descending ordering, as explained below. Note that all the seven leading modes in both flows are associated with $`q_x=0`$. This is surely due to the relatively short channel length in our computation. It is known that in longer channels the “roll” modes become oscillatory modes with a finite value of $`q_x`$.

For the sake of completeness we present in Fig. 9 the portrait of a higher order viscoelastic mode, namely $`\stackrel{~}{𝒒}=(0,\mathrm{\hspace{0.17em}3})`$, $`\mathrm{p}=15`$. We see that for this mode $`\mathrm{\Psi }_z`$ is essentially purely imaginary, whereas $`\mathrm{\Psi }_x`$ and $`\mathrm{\Psi }_y`$ are essentially purely real. In fact, this particular mode is very low in energy, and we present it just to demonstrate that the numerics is still not noisy even for rather low energy modes.

#### IV.2.2 Sub-dominant modes and energy redistribution

The full comparison between the first 30 modes of the Newtonian and viscoelastic flows can be seen in Table 1, in which these modes are listed together with their energy (in percent of the total sum of energies). Also included in the table is the cumulative sum of energies up to the mode listed.

To understand the Table we recall that for each value of $`𝒒`$ the modes are ordered by their energy and labeled by p. The over-all energy label is $`n__\mathrm{N}`$ or $`n_{_{\mathrm{VE}}}`$. In the columns of the viscoelastic modes, instead of writing $`n_{_{\mathrm{VE}}}`$ in obvious increasing order, we provided the value of the Newtonian index $`n__\mathrm{N}`$ for which the viscoelastic mode has the highest matching parameter. Thus for example the second leading viscoelastic mode with $`n_{_{\mathrm{VE}}}=2`$ matches almost perfectly the sixth Newtonian mode $`n__\mathrm{N}=6`$, etc.

One glaring difference between the two flows that stands out from Table 1 is the energy concentration in the few most dominant modes of the viscoelastic flow compared to the flat distribution of energy in the Newtonian flow. For example, the first, second and third dominant viscoelastic modes have all twice the energy of the corresponding first, second and third dominant Newtonian modes. The sum of the first four viscoelastic energies contain as much as as the first ten Newtonian modes. The first 15 viscoelastic modes already contain more than 50% of the energy, whereas one needs to collect as many as 80 Newtonian modes to reach the same fraction of the total energy. In Fig. 10 we display the sum $`(n)=_{j=1}^nE(j)`$ as a function of $`n`$ for both flows, where the sum is computed in the respective energy ordering. We note that the difference between the two curves is established within the first ten modes; for larger values of $`n`$ the lines are almost parallel, indicating a similar energy distribution between the less dominant modes.

### IV.3 The energy spectrum

An interesting and illuminating way to discuss the difference between the two flows is provided by the energy spectrum. It was proposed in 90Ge  that the main difference between the spectra of the two flows should appear in the position of the dissipative scale that separates a spectral power-law from exponential decay. The increase of the dissipative scale was indeed observed in recent DNS of homogoneous isotropic turbulence in the FENE-P model 02ACBP . However, the present results indicate that important changes involve the energy containing scales.

In Fig. 11 the energy is plotted for the two flows as a function of $`n`$. The Newtonian case displays a spectral plateau for the most dominant modes which crosses over to a power law for $`n__\mathrm{N}8`$. In the case of the viscoelastic flow the plateau is dramatically higher, but also the power law changes its slope. It appears that the whole spectral curve is tilted in favor of the energy containing modes and on the expense of the lower energy modes and the dissipative modes.

The difference between the power laws can be made clearer by comparing with what can be obtained from the Kolmogorov law for high $`k`$-vectors. For relatively large $`k`$ we can expect the flow to isotropize 99ALP ; 02LPT , and Fourier modes would again become “best”. In this asymptotic situation (neglecting intermittency corrections) the spectrum $`E(k)`$ is expected to be $`E(k)k^{11/3}`$. In an isotropic environment we also expect that $`n`$ is proportional to the volume of the sphere of radius $`k`$, i.e. $`nk^3`$. Thus we can expect for large values of $`n`$

$$E_nn^{11/9},n\text{ large}.$$
(39)
The law proposed in Eq. (39) is displayed as the dashed-dotted line in Fig. 11. We see that it is in rough agreement with the power law section of the Newtonian spectrum, but it is certainly not in agreement with the viscoelastic spectrum which, as said before, is becoming steeper on the whole. This is a clear demonstration of the increase in energy of the energy containing modes on the expense of the others. We propose that even the power law section of the viscoelastic spectrum will not have a “universal” slope, but rather a slope that depends on the concentration of the polymer and the degree of drag reduction.

What emerges from this study is that understanding drag reduction lies in the reordering and energy redistribution of the energy containing modes. To examine these phenomena further we consider now the energy contents of the viscoelastic modes ordered by $`n__\mathrm{N}`$ instead of their own ordering. This plot is a very vivid graphic representation of the data of Table 1, stressing very strongly the concentration of energy in the dominant viscoelastic modes, which are however rearranged in dominance compared to the Newtonian case.

Next we want to reiterate that one should focus on the most energetic modes. Noticing from Table 1 that all the dominant modes in both flows are associated with $`\stackrel{~}{q}_x=0`$, we consider next these modes as a function of $`\stackrel{~}{q}_z`$ for $`p=1,2`$. In Fig. 13 we show the energy of these modes. We note that the dramatic redistribution of energy occurs only for the most energetic modes with $`p=1`$; The less energetic modes are not affected much. Already the modes with $`p=2`$ are seen in the figure to be affected in a negligible way. In our opinion this is a clear message that to understand drag reduction we need to understand the rearrangement of the energy containing modes. For our flow configuration the most relevant are the modes which are space homogeneous in the spanwise direction. It should be stressed however that different geometries, and even channel geometry with different aspect ratios, may bring forth other modes as the most relevant ones. Nevertheless we expect that drag reduction would always be associated with a substantial increase in the energy containing modes, whatever these are for a given flow configuration.

## V conclusions

In this paper we initiated a systematic study of drag reduction on the basis of DNS of the FENE-P model. We investigated simulations of Newtonian and viscoelastic flows in channel configuration at the same friction Reynolds number. Our main aim is to understand the mechanism of drag reduction. Since drag reduction involved modifications of the mean flow and of the large scale gradients in the flow, we are motivated to understand more the energy containing modes, rather than focusing only on phenomena of small scales. To this aim we have found first the list of empirical modes that represent the velocity field in an optimal way in an energy decreasing ordering. The first important discovery was that this list contains the very same modes for the Newtonian and viscoelastic flows. We propose that this finding will offer a huge simplification in any future theory of drag reduction. While we cannot offer a clean explanation of this finding, we can proceed at this time taking the approximate invariance of the modes as an empirical fact. What needs to be understood is just the energy distribution and reordering of the modes in the viscoelastic case.

We should stress that the point of view proposed here differs in a fundamental way from the approach presented for example in Refs. 69Lu ; 90Ge ; 01BFL . The thinking there focuses on the energy cascade in the turbulent flow, and on the modification of the small dissipative scales. Roughly speaking, the maximal gradients of the velocity are estimated from balancing the RHS of Eq. (4). Neglecting the statistical correlation between the conformation tensor and the velocity gradient, one can estimate the maximal velocity gradient by $`u/r1/\tau _\mathrm{p}`$. Thus the matching of the polymer relaxation time with an eddy turn over time is used to predict a decrease in the maximal velocity gradient which is interpreted as drag reduction. We take an exception to this approach. First, a careful analysis of the space dependence of the conformation tensor shows that it is highly correlated with the velocity gradients (see for example 98DSB ; 02ACBP ; 02BP . Thus the estimate taken above is questionable at best. But moreover, we have shown that the main changes between the Newtonian and the viscoelastic modes occur in the energy containing modes. The energy containing modes are highly anisotropic, they are not Fourier modes, and their connection to the small scales where the isotropised Kolmogorov picture is tenable is very unclear. A theory that assumes a K41 spectrum down to a modified viscous scale does not appear tenable for the FENE-P flow, as seen in Fig. 11. We propose that a theory of drag reduction entails an understanding of the relative energy of the modes that characterize the largest scales of the flow.

To understand the relative ordering and the relative energy of the most dominant modes one needs to study the energy intake by these modes from the mean flow 01GLP ; 02GLP , the energy exchange between the modes, and the energy exchange with the viscoelastic subsystem represented by the conformation tensor field $`𝑹(𝒓,t)`$ 02ACBP ; 02BP . Such an investigation calls for measuring additional statistical objects like 3rd order correlation functions. The necessary simulation are in progress and the results will be published elsewhere.

###### Acknowledgements.

We are indebted to Anna Pomyalov for helpful suggestions in the course of this work. This work has been supported in part by the European Commission under a TMR grant, the German Israeli Foundation, the Minerva Foundation, Munich, Germany, the Israeli Science Foundation and the Naftali and Anna Backenroth-Bronicki Fund for Research in Chaos and Complexity.


format: text

text:
# Factorizing one-loop contributions to two-loop Bhabha scattering and automatization of Feynman diagram calculations

## 1 DIANA

We are moving in the direction of two-loop Bhabha scattering, which is extremely important, in particular at higher energies, for the luminosity determination of the coming accellerators. By means of DIANA  we produced the following one-loop diagrams:

Factorizing one-loop contributions are obtained e.g. in mass-renormalization: in any Feynman diagram we replace

$$\frac{1}{k^2(m^2+\delta (m^2))}\frac{1}{k^2m^2}(1+\frac{\delta (m^2)}{k^2m^2}),$$

thus obtaining ‘dotted’ diagrams.

## 2 Calculation of 1-loop integrals

The evaluation of one-loop integrals as used here was introduced in detail in . The method exists in first reducing tensor integrals to scalar ones in higher dimension and higher indices (i.e. higher powers of scalar propagators). Before further reducing these scalar integrals arready the on-shell conditions can be implemented and the amplitude decomposition performed. The reduction to ‘master integrals’ is performed by the use of recursion relations. One of these, reducing the space time dimension, is given by:

$`(d{\displaystyle \underset{i=1}{\overset{n}{}}}\nu _i+1)\left(\right)_nI_n^{(d+2)}=`$
$`\left[\left({\displaystyle \genfrac{}{}{0pt}{}{0}{0}}\right)_n{\displaystyle \underset{k=1}{\overset{n}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{0}{k}}\right)_n𝐤^{}\right]I_n^{(d)},`$
where the objects $`\left(\genfrac{}{}{0pt}{}{i}{j}\right)_n`$ are Gram determinants. As a particular problem occurred in this connection zero Gram determinants , i.e. $`\left(\genfrac{}{}{0pt}{}{0}{0}\right)_n=0:`$

$$(d\underset{i=1}{\overset{n}{}}\nu _i+1)\left(\right)_nI_n^{(d+2)}=\underset{k=1}{\overset{n}{}}\left(\genfrac{}{}{0pt}{}{0}{k}\right)_n𝐤^{}I_n^{(d)}.$$

To increase again the dimension $`d`$, one then uses the relation

$$\underset{j=1}{\overset{n}{}}\nu _j𝐣^+I_n^{(d+2)}=I_n^{(d)}.$$

## 3 Decomposition of the diagrams into Amplitudes

In order to avoid $`\gamma _5`$, we decompose the amplitdes in terms of the following 12 structures:

$$\begin{array}{cccccccc}O_1=II,O_2=\widehat{p}_4I,O_3=I\widehat{p}_2,\hfill & & & & & & & \\ O_4=\widehat{p}_4\widehat{p}_2,O_5=\gamma _\mu \gamma _\mu ,O_6=\gamma _\mu \widehat{p}_4\gamma _\mu ,\hfill & & & & & & & \\ O_7=\gamma _\mu \gamma _\mu \widehat{p}_2,O_8=\gamma _\mu \gamma _\nu \gamma _\nu \gamma _\mu ,\hfill & & & & & & & \\ O_9=\gamma _\mu \gamma _\nu \widehat{p}_4\gamma _\nu \gamma _\mu ,O_{10}=\gamma _\mu \gamma _\nu \gamma _\nu \gamma _\mu \widehat{p}_2,\hfill & & & & & & & \\ O_{11}=\gamma _\mu \gamma _\nu \widehat{p}_4\gamma _\nu \gamma _\mu \widehat{p}_2,O_{12}=\gamma _\mu \gamma _\nu \gamma _\rho \gamma _\rho \gamma _\nu \gamma _\mu .\hfill & & & & & & & \end{array}$$

Thus the on-shell diagrams read

$$\mathrm{𝐃𝐢𝐚𝐠𝐫𝐚𝐦}=\underset{𝐣=\mathrm{𝟏}}{\overset{\mathrm{𝟏𝟐}}{}}𝐀_𝐣𝐎_𝐣.$$

with amplitudes $`𝐀_𝐣`$. Finally we obtained ‘crossing relations’ like e.g. crossing diagrams 8 and 9 into each other:

$`A_1^{(9)}=A_1^{(8)}2m_eA_7^{(8)}4m_eA_9^{(8)}+2dA_8^{(8)}`$
$`A_2^{(9)}=A_2^{(8)}+4A_7^{(8)}2(d4)A_9^{(8)}8m_eA_{11}^{(8)}`$
$`A_3^{(9)}=A_3^{(8)}+2A_6^{(8)}2(d2)A_{10}^{(8)}4m_eA_{11}^{(8)}`$
$`A_4^{(9)}=A_4^{(8)}2(d2)A_{11}^{(8)}`$
$`A_5^{(9)}=A_5^{(8)}+2m_eA_6^{(8)}+4m_eA_7^{(8)}+8m_eA_9^{(8)}`$
$`+4m_eA_{10}^{(8)}12m_e^2A_{11}^{(8)}(6d4)A_{12}^{(8)}`$
$`A_6^{(9)}=4m_eA_{11}^{(8)}+A_6^{(8)}`$
$`A_7^{(9)}=4m_eA_{11}^{(8)}+A_7^{(8)}`$
$`A_8^{(9)}=A_8^{(8)}`$
$`A_9^{(9)}=A_9^{(8)}`$
$`A_{10}^{(9)}=A_{10}^{(8)}`$
$`A_{11}^{(9)}=A_{11}^{(8)}`$
$`A_{12}^{(9)}=A_{12}^{(8)}`$
and exchanging t and u. These relations also hold for crossing diagrams 7 to 10 (s and u exchanged). For the crossing of 7 and 8 (s and t exchanged) only a change of sign occurs in the amplitudes. Thus:  we need to calculate only one selfenergy, vertex and box.

Further relations between amplitudes of one diagram (obtained by solving a system of equations) are:

$`A_3=A_2`$
$`A_7=A_6`$
$`A_{10}=A_9`$
$`A_8={\displaystyle \frac{1}{4}}A_1{\displaystyle \frac{m_e}{2}}A_2+(A_{12})`$
$`A_9={\displaystyle \frac{1}{4m_e}}A_1{\displaystyle \frac{1}{2}}A_2{\displaystyle \frac{1}{2}}A_6`$
$`A_{11}={\displaystyle \frac{1}{4m_e^2}}A_1{\displaystyle \frac{1}{2m_e}}A_2`$
Thus: there are only 6 indpendent amplitudes.

The amplitudes for any of the diagrams, even for non-zero electron mass, are relatively short. They were presented in the talk and can be obtained from ACAT’s webpage . What remains to do is to expand the master integrals, given in arbitrary $`d`$, in $`\epsilon ,(d=42\epsilon )`$ in order to obtain properly all finite contributions. This work is presently being done.


text:
# 1 Introduction

## 1 Introduction

Even if nobody will ever see a magnetic monopole, there is surely much to be gained by studying the theory of monopoles. In the seminal paper of Julia and Zee  the existence of a profound connection between the angular momentum and the electric charge in a Yang-Mills-Higgs (YMH) theory has been suggested. A well know result in classical electrodynamics implies that in the presence of a magnetic monopole an electrically charged particle acquires an additional angular momentum oriented toward the position of the monopole and with a magnitude equal to the product of the electric and magnetic charges . In the nonabelian case the electric charge distribution and the magnetic charge distribution are not separated. To compute the angular momentum for this type of configuration we have to solve extremely complicated partial differential equations.

Surprisingly, a few years ago it has been shown that Julia-Zee dyons do not admit slowly rotating excitations . However the absence or presence of rotating perturbative solutions, though indicative, is not in general conclusive to establish the absence or presence of exact solutions. For instance in boson star models slowly rotating perturbative solutions are absent, but solutions with a quantized angular momentum exist  (see also the discussion in ). It is therefore important to identify non-perturbative criteria to settle such questions.

In a recent paper a general relation for the total angular momentum in an Einstein-Yang-Mills-Higgs theory as a surface integral in terms of Yang-Mills (YM) fields has been derived . Within an ansatz which satisfies some extra symmetries it has been found that there are no asymptotically flat rotating dyon solutions, while a monopole-antimonopole configuration with zero net magnetic charge possesses an angular momentum proportional to the electric charge. However Ref. leaves unsolved the question of possible existence of finite energy solutions with a nonvanishing total angular momentum and net magnetic charge for a more general ansatz.

In this paper we argue that, for the most general axially symmetric YMH monopole configuration, the total angular momentum is zero. We adress also the same question for a different asymptotic structure of spacetime and present arguments for the existence of axially symmetric dyon solutions in an anti-de Sitter (AdS) spacetime. Qualitatively, the behavior of AdS solutions is very similar to that corresponding to Minkowski spacetime configurations. However, we argue that these solutions present also a vanishing total angular momentum.

## 2 General ansatz and total angular momentum

To simplify the general picture, in this paper we ignore the effects of gravity. We follow also the notations and conventions used by Hartmann, Kleihaus and Kunz in . The action for a non-Abelian $`SU(2)`$ gauge field coupled to a triplet Higgs field with the usual potential $`V(\mathrm{\Phi })=\frac{\lambda }{8}Tr(\mathrm{\Phi }^2\eta ^2)^2`$ is

$$S=d^4x\sqrt{g}[Tr\{\frac{1}{2}F_{\mu \nu }F^{\mu \nu }\}+Tr\{\frac{1}{2}D_\mu \mathrm{\Phi }D^\mu \mathrm{\Phi }\}+V(\mathrm{\Phi })],$$
(1)
where the field strength tensor is

$`F_{\mu \nu }`$ $`=`$ $`_\mu A_\nu _\nu A_\mu +i[A_\mu ,A_\nu ],`$
and the covariant derivative

$`D_\mu `$ $`=`$ $`_\mu +i[A_\mu ,].`$
Varying the action (1) with respect to $`A_\mu `$ and $`\mathrm{\Phi }`$ we have the field equations

$`{\displaystyle \frac{1}{\sqrt{g}}}D_\mu (\sqrt{g}F^{\mu \nu })`$ $`=`$ $`{\displaystyle \frac{1}{4}}i[\mathrm{\Phi },D^\nu \mathrm{\Phi }],`$ (2)
$`{\displaystyle \frac{1}{\sqrt{g}}}D_\mu (\sqrt{g}D^\mu \mathrm{\Phi })+\lambda (\mathrm{\Phi }^2\eta ^2)\mathrm{\Phi }`$ $`=`$ $`0.`$ (3)
Varying the Lagrangian with respect to the metric $`g_{\mu \nu }`$ yields the energy-momentum tensor

$`T_{\mu \nu }=2Tr\{F_{\mu \alpha }F_{\nu \beta }g^{\alpha \beta }{\displaystyle \frac{1}{4}}g_{\mu \nu }F_{\alpha \beta }F^{\alpha \beta }\}+Tr\{{\displaystyle \frac{1}{2}}D_\mu \mathrm{\Phi }D_\nu \mathrm{\Phi }{\displaystyle \frac{1}{4}}g_{\mu \nu }(D_\alpha \mathrm{\Phi }D^\alpha \mathrm{\Phi }\}+V(\mathrm{\Phi })g_{\mu \nu }.`$ (4)
The energy density of a solution of the YMH equations is given by the $`tt`$component of the energy-momentum tensor. Integration over all space yields the total energy

$`E`$ $`=`$ $`{\displaystyle }T_t^t\sqrt{g}d^3x={\displaystyle _0^{\mathrm{}}}dr{\displaystyle _0^\pi }d\theta {\displaystyle _0^{2\pi }}d\phi \sqrt{g}Tr\{F_{r\theta }F^{r\theta }+F_{r\phi }F^{r\phi }+F_{\theta \phi }F^{\theta \phi }`$
$``$ $`F_{rt}F^{rt}F_{\theta t}F^{\theta t}F_{\phi t}F^{\phi t}+{\displaystyle \frac{1}{4}}(D_r\mathrm{\Phi }D^r\mathrm{\Phi }+D_\theta \mathrm{\Phi }D^\theta \mathrm{\Phi }+D_\phi \mathrm{\Phi }D^\phi \mathrm{\Phi }D_t\mathrm{\Phi }D^t\mathrm{\Phi })+{\displaystyle \frac{\lambda }{8}}(\mathrm{\Phi }^2\eta ^2)^2\},`$
where $`r,\theta ,\phi `$ are the usual spherical coordinates. For the energy integral to converge, each of the terms in (2) must vanish at large $`r`$. The finite energy condition imposed on this expression will give us the acceptable asymptotic behavior of the gauge functions and Higgs field as $`r\mathrm{}`$. The additional requirements to have a finite, locally integrable energy density impose boundary conditions at the origin and on the $`z`$axis.

The total angular momentum is

$`J`$ $`=`$ $`{\displaystyle T_\phi ^t\sqrt{g}d^3x}={\displaystyle 2Tr\{F_{r\phi }F^{rt}+F_{\theta \phi }F^{\theta t}+D_\phi \mathrm{\Phi }D^t\mathrm{\Phi }\}\sqrt{g}d^3x}.`$ (6)
We mention also the expression for the electric and magnetic charges derived by using the ’t Hooft field strength tensor

$`Q_e`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle _{\mathrm{}}}𝑑S_\mu Tr\{\widehat{\mathrm{\Phi }}F_{\mu t}\},`$ (7)
$`Q_m`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle _{\mathrm{}}}𝑑S_\mu {\displaystyle \frac{1}{2}}ϵ_{\mu \nu \alpha }Tr\{\widehat{\mathrm{\Phi }}F_{\nu \alpha }\},`$ (8)
where the integration is over a surface at spatial infinity and $`\widehat{\mathrm{\Phi }}=\mathrm{\Phi }/|\mathrm{\Phi }|`$.

We are interested in static, axially symmetric finite energy solutions of the equations (2)-(3). As proven in , for static axially symmetric configurations, the volume integral (6) can be converted into a surface integral in terms of Yang-Mills potentials. The symmetry of the gauge field under a spacetime symmetry means that the action of an isometry can be compensated by a suitable gauge transformation . For the time translational symmetry, we choose a natural gauge such that $`A/t`$=0. However, a rotation around the $`z`$axis can be compensated by a gauge rotation

$`_\phi A=D\mathrm{\Psi },`$ (9)
and therefore

$`F_{\mu \phi }=`$ $`D_\mu W,`$ (10)
$`D_\phi \mathrm{\Phi }=`$ $`i[W,\mathrm{\Phi }],`$
where $`W=A_\phi \mathrm{\Psi }`$.

Using the potential $`W`$ we find

$`T_\phi ^t`$ $`=`$ $`2Tr\left\{(D_rW)F^{rt}+(D_\theta W)F^{\theta t}+{\displaystyle \frac{i}{4}}[W,\mathrm{\Phi }]D^t\mathrm{\Phi }\right\}`$ (11)
$`=`$ $`2Tr\{{\displaystyle \frac{1}{\sqrt{g}}}D_r(WF^{rt}\sqrt{g})+{\displaystyle \frac{1}{\sqrt{g}}}D_\theta (WF^{\theta t}\sqrt{g})`$
$`W({\displaystyle \frac{1}{\sqrt{g}}}D_r(\sqrt{g}F^{rt})+{\displaystyle \frac{1}{\sqrt{g}}}D_\theta (\sqrt{g}F^{\theta t}))+{\displaystyle \frac{i}{4}}[W,\mathrm{\Phi }]D^t\mathrm{\Phi }\}.`$
As a consequence of the YM equations (2) and making use of the fact that the trace of a commutator vanishes we obtain

$`T_\phi ^t=2Tr\left\{{\displaystyle \frac{1}{\sqrt{g}}}_\mu (WF^{\mu t}\sqrt{g})\right\}.`$ (12)
Thus, ignoring possible singularities the expression of the total angular momentum is

$`J`$ $`=`$ $`{\displaystyle _{\mathrm{}}}2Tr\{WF^{\mu t}\}𝑑S_\mu `$ (13)
$`=`$ $`2\pi \underset{r\mathrm{}}{lim}{\displaystyle _0^\pi }𝑑\theta \mathrm{sin}\theta ^{}r^2[W^{(r)}F^{rt(r)}+W^{(\theta )}F^{rt(\theta )}+W^{(\phi )}F^{rt(\phi )}].`$
Note that this relation has been derived without fixing the gauge. It greatly simplifies the evaluation of the angular momentum of a regular YMH configuration, since we need only the asymptotics of some gauge functions. In what follows, we use (13) to prove that the total angular momentum of a finite energy dyon solution in flat space is zero.

The method of proof is, in principle, quite straightforward. From the assumptions of finite energy density and finite total energy, it follows that every term in the energy integral (2) must be finite. The most crucial requirement is that $`D_\mu \mathrm{\Phi }`$ vanish sufficiently rapidly as $`r\mathrm{}`$. When evaluating (13) we find a vanishing $`J`$ for this behavior.

The construction of an axially symmetric YMH ansatz has been discussed by many authors starting with the pioneering papers by Manton  and Rebbi and Rossi . The most general axially symmetric Yang-Mills ansatz contains 12 functions: 9 magnetic and 3 electric potentials and can be easily obtained in cylindrical coordinates $`x^\mu =(\rho ,\phi ,z,t`$) (with $`\rho =r\mathrm{sin}\theta ,z=r\mathrm{cos}\theta )`$

$`A_\mu ={\displaystyle \frac{1}{2}}A_\mu ^\rho (\rho ,z)\tau _\rho ^n+{\displaystyle \frac{1}{2}}A_\mu ^\phi (\rho ,z)\tau _\phi ^n+{\displaystyle \frac{1}{2}}A_\mu ^z(\rho ,z)\tau _z,`$ (14)
where the only $`\phi `$-dependent terms are the $`SU(2)`$ matrices (composed of the standard $`(\tau _x,\tau _y,\tau _z)`$ Pauli matrices)

$`\tau _\rho ^n`$ $`=`$ $`\mathrm{cos}n\phi \tau _x+\mathrm{sin}n\phi \tau _y,`$
$`\tau _\phi ^n`$ $`=`$ $`\mathrm{sin}n\phi \tau _x+\mathrm{cos}n\phi \tau _y.`$ (15)
We can reduce this general ansatz by imposing some extra symmetries on the fields. To our knowledge there are no analytical or numerical results for the most general ansatz.

Transforming to spherical coordinates, it is convenient to introduce, without any loss of generality, a new $`SU(2)`$ basis $`(\tau _r^n,\tau _\theta ^n,\tau _\phi ^n)`$, with

$`\tau _r^n`$ $`=`$ $`\mathrm{sin}\theta \tau _\rho ^n+\mathrm{cos}\theta \tau _z,`$
$`\tau _\theta ^n`$ $`=`$ $`\mathrm{cos}\theta \tau _\rho ^n\mathrm{sin}\theta \tau _z.`$ (16)
Thus, similar to (14), we find the general expression

$`A_\mu ={\displaystyle \frac{1}{2}}A_\mu ^r(r,\theta )\tau _r^n+{\displaystyle \frac{1}{2}}A_\mu ^\theta (r,\theta )\tau _\theta ^n+{\displaystyle \frac{1}{2}}A_\mu ^\phi (r,\theta )\tau _\phi ^n,`$ (17)
where $`A_\mu ^adx^\mu =A_r^adr+A_\theta ^ad\theta +A_\phi ^ad\phi +A_t^adt`$. For the Higgs field we have a similar relation $`\mathrm{\Phi }=\mathrm{\Phi }^a(r,\theta )\tau _a^n`$. The gauge invariant quantities expressed in terms of these functions will be independent of the azimutal angle $`\phi `$ and hence axially symmetric. The winding number $`n`$ corresponds to the topological charge of the solutions and is a constant of motion.

For this parametrization

$$\mathrm{\Psi }=n\frac{\tau _z}{2}=n\mathrm{cos}\theta \frac{\tau _r^n}{2}n\mathrm{sin}\theta \frac{\tau _\theta ^n}{2}.$$
(18)
The Higgs field must be nonvanishing at spatial infinity in order that the potential energy be zero there. Thus, the spontaneous symmetry-breaking mechanism requires

$$\underset{r\mathrm{}}{lim}\mathrm{\Phi }^a\mathrm{\Phi }^a=\eta ^2.$$
(19)
As proven in , any regular axially symmetric magnetic charge distribution can be located only at isolated points situated on the axis of symmetry, with equal and opposite values of the charge at alternate points. In particular, if only one sign of the charge is allowed, all the charge must be concentrated at a single point.

For a multimonopole solution containing only magnetic charges with the same sign, the suitable boundary conditions for the Higgs field at infinity are $`\mathrm{\Phi }^r=\eta ,\mathrm{\Phi }^\theta =\mathrm{\Phi }^\phi =0`$. Note that a more general asymptotic behavior of the Higgs field is allowed by (19), but this would correspond to composed configurations containing both monopoles and antimonopoles.

The requirement of finite total energy implies the decay of the $`F_{rt}^a`$ as $`r\mathrm{}`$ is faster than $`1/r^{1.5}`$. Also, the contribution $`(D_\phi \mathrm{\Phi })^a(D_\phi \mathrm{\Phi })^ag^{\phi \phi }`$ term to the energy is finite if $`(D_\phi \mathrm{\Phi })^a`$ approaches zero to large $`r`$ sufficiently rapid. This implies the fall-off conditions $`(W_\theta ,W_\phi )1/r^{0.5+ϵ}`$ for large $`r`$ (with $`ϵ>0`$). Therefore the last two terms in (13) give null contribution to the total angular momentum. The contribution of the $`A_\phi ^r`$ to the total angular momentum is also zero, since it vanishes asymptotically for a regular configuration.

By using these results we find

$$\underset{r\mathrm{}}{lim}Tr(r^2WF_{rt})=\frac{nQ_e}{2}\mathrm{cos}\theta $$
(20)
where $`Q_e`$ is the electric charge of the configuration given by (7). Therefore the total angular momentum of the solution is clearly zero. Thus, an axially symmetric regular dyon has a vanishing total angular momentum.

## 3 Axially symmetric dyons in AdS background

This nonexistence result can be partially attributed to the asymptotic structure of spacetime. One may ask whether it remains valid in the presence of a cosmological constant $`\mathrm{\Lambda }`$.

However, up to date, all known axially symmetric dyon regular solutions have been obtained in a fixed Minkowski background. In what follows we consider dyon solutions in an AdS spacetime. Further motivation for this attempt comes from the surprising behavior found for a $`SU(2)`$gauge field in the presence of a negative cosmological constant. The existence of stable particle-like and black hole solutions in an Einstein-Yang-Mills theory with $`\mathrm{\Lambda }<0`$ is the most exciting result .

The global existence of a solution of the Cauchy problem for the YMH equations in AdS spacetime is considered in . Some analytical properties of spherically symmetric monopole and dyon solutions for YMH theory with a negative cosmological constant have been discussed in . Unlike the flat spacetime case, in an AdS spacetime there are no analytic solutions that might be used as a guiding line. Moreover, when a cosmological constant is included (no matter how small this constant is) no solution close to the BPS configuration can be found. As proven in , it is the change in the asymptotic behavior of the Higgs field that prevents such a solution. Numerical monopole solutions have been exhibited in , where the effects of gravity are also included. A distinctive feature of AdS solutions concerns the asymptotic behavior of the fields. When $`\mathrm{\Lambda }<0`$, the Higgs field approaches its vacuum expectation value faster than in the flat space case. The radius of the monopole core decreases, as the magnetic field concentrates near the origin.

In what follows, we present numerical arguments for the existence of spherically- and axially symmetric dyon solutions in an AdS spacetime with line-element

$$ds^2=\frac{dr^2}{1\frac{\mathrm{\Lambda }}{3}r^2}+r^2(d\theta ^2+r^2\mathrm{sin}^2\theta d\phi ^2)(1\frac{\mathrm{\Lambda }}{3}r^2)dt^2.$$
(21)
Searching for general axially symmetric solutions within the most general ansatz is a difficult task. In our approach we consider a reduced ansatz used also to obtain flat space - axially symmetric dyons

$`A_r^r`$ $`=`$ $`A_r^\theta =A_\theta ^r=A_\theta ^\theta =A_\phi ^\phi =\mathrm{\Phi }^\phi =0,`$
$`A_r^r`$ $`=`$ $`{\displaystyle \frac{H_1}{r}},A_\theta ^r=1H_2,A_\phi ^r=n\mathrm{sin}\theta H_3,A_\phi ^\theta =n\mathrm{sin}\theta (1H_4),`$ (22)
$`A_t^r`$ $`=`$ $`\eta H_5,A_t^\theta =\eta H_6,\mathrm{\Phi }^r=\eta \mathrm{\Phi }_1,\mathrm{\Phi }^\theta =\eta \mathrm{\Phi }_2.`$
This ansatz satisfies some discrete symmetries  and was used by various authors when discussing numerically solutions . The boundary conditions at infinity consistent with the requirements of regularity, finite energy and symmetry, are

$$H_i=0,i=1,2,3,4,6;H_5=\alpha ;\mathrm{\Phi }_1=1;\mathrm{\Phi }_2=0$$
(23)
at infinity, and

$$H_i=0,i=1,3,5,6;H_2=H_4=1,\mathrm{\Phi }_1=\mathrm{\Phi }_2=0,$$
(24)
at the origin. Given the parity reflection symmetry, we need to consider solutions only in the region $`0\theta \pi /2`$; on the $`z`$\- and $`\rho `$-axis the functions $`H_1,H_3,H_6,\mathrm{\Phi }_2`$ and the derivatives $`_\theta H_2,_\theta H_4,_\theta H_5`$ and $`_\theta \mathrm{\Phi }_1`$ are to vanish. To fix the residual abelian gauge invariance we choose the usual gauge condition $`r_rH_1_\theta H_2=0`$ .

Within this ansatz, we have found dyon solutions for any value of the cosmological constant. We change to dimensionless coordinates and Higgs field by rescaling $`rr/\eta `$, $`\mathrm{\Lambda }\mathrm{\Lambda }\eta ^2`$ and $`\mathrm{\Phi }\eta \mathrm{\Phi }`$, respectively. The numerical calculations for $`n>1`$ were performed with the software package CADSOL/FIDISOL, based on the Newton-Raphson method . We use the general relations (2)-(8) to compute the energy, angular momentum and charges of these configurations. The solutions with $`n=1`$ corresponds to spherically symmetric dyons with $`H_1=H_3=\mathrm{\Phi }_2=0,H_2=H_4=w(r),\mathrm{\Phi }_1=H(r),H_5=J(r)`$, and generalize the Julia-Zee solution  for a negative cosmological constant. Typical spherically symmetric solutions are presented in Fig. 1.

Although the results we find are broadly similar to those valid for the $`\mathrm{\Lambda }=0`$ case, there are some differences. The most interesting difference, not discussed in , is that the magnitude of the electric potential $`A_t`$ at infinity is no longer restricted. In an asymptotically Minkowski spacetime, the constant $`\alpha `$ in eq. (23) is restricted to $`\alpha 1`$. For $`\alpha >1`$ some gauge field functions become oscillating instead of asymptotically decaying. The $`A_t^a`$ components of the gauge filed act like an isotriplet Higgs field with negative metric, and by themselves would cause the other components of the gauge field to oscillate rather than decrease exponentially as $`r\mathrm{}`$ . However, in a AdS spacetime, finite energy solutions with arbitrary $`\alpha `$ are allowed. The Higgs field forces $`H_1=H_2=H_4=H_6=0`$ asymptotically, but does not restrict the value of $`H_5`$. In the limit of vanishing Higgs field ($`\eta 0`$) we recover the (pure-) YM solutions discussed in  (in the absence of the Higgs field the value at infinity of the magnetic potentials $`H_2,H_4`$ is arbitrary). Another distinctive feature of the AdS solutions is a smaller dyon radius as compared to the $`\mathrm{\Lambda }=0`$ case.

Apart from these properties, qualitatively, the Higgs field and Yang-Mills field behavior is very similar to that corresponding to Minkowski spacetime dyons. We notice a similar shape for the functions $`H_i`$ and $`\mathrm{\Phi }_i`$ and also for the energy density. Increasing in size at the origin for increasing $`|\mathrm{\Lambda }|`$, the energy density is localized in a decreasing region of space. The energy of the multidyons is on order of $`n`$ times the corresponding one-dyon energy. In Fig. 2 the energy and the electric charge of the solutions is plotted as a function of the cosmological constant for various winding numbers. The toroidal shape familiar from the energy density of the dyon solutions in flat space is retained for solutions in AdS space, as is illustrated in Fig. 3.

Most of our results have been obtained in the Prasad-Sommerfeld limit, but similar conclusions hold for a nonvanishing Higgs self-coupling. Further details on these solutions will be presented elsewhere.

Returning to the question of angular momentum, by using the ansatz (3) with the boundary conditions (23) we find again

$$J=2\pi \underset{r\mathrm{}}{lim}_0^\pi 𝑑\theta \mathrm{sin}\theta r^2W^{(r)}F_{rt}^{(r)}=2\pi nQ_e_0^\pi 𝑑\theta \mathrm{sin}\theta \mathrm{cos}\theta =0.$$
(25)
It is worthwhile to notice that these solutions carry a nonvanishing angular momentum density, even in the Prasad-Sommerfeld limit, as ilustrated in Fig. 3c (see also Fig. 4). The fact that $`J=0`$ can be atributed again to the extra symmetries of the reduced ansatz, since as $`\theta \pi \theta `$ we have $`H_1H_1,H_2H_2,H_3H_3,H_4H_4,H_5H_5,H_6H_6`$ which implies the vanishing of the integral (6). We find also that the YM dyon solutions discussed in  have a vanishing angular momentum.

Evidently, dyon solutions within the most general ansatz may also exist. The construction of such solutions represents a difficult challenge. However, we argue that even for the most general ansatz, an AdS finite energy dyon solution presents a vanishing total angular momentum.

The asymptotic behavior of the functions $`A_\mu ^a,\mathrm{\Phi }^a`$ is complicated by the presence of a negative cosmological constant. Similar to the $`n=1`$ case , the transversal components of the YM potentials present a complicated power decay, instead of the expected exponential law. Nevertheless, following the analysis in the $`\mathrm{\Lambda }=0`$ case, we find that the required fall-off conditions for a finite total energy are again $`F_{rt}^a1/r^{1.5+ϵ}`$, $`W_\theta `$ and $`W_\phi 1/r^{0.5+ϵ}`$, while $`A_\phi ^r`$ shold also vanish for a regular configuration. Thus, after replacing in (12) we find again $`J=0`$.

## 4 Further discussion

We presented arguments that the results obtained within a perturbative approach in  remain valid in the general case, for the most general YMH axially symmetric ansatz. A flat space axially symmetric dyon configuration necessarily possesses a vanishing total angular momentum.

We emphasize the role played by the Higgs field boundary condition at infinity in obtaining this nonexistence result. We can of course consider a different behavior at infinity, consistent also with (19) and finite energy conditions, corresponding, however, to composed configurations. A sufficiently general asymptotic behavior of the Higgs field at infinity will read

$$\mathrm{\Phi }^r=\eta \mathrm{cos}f(\theta ),\mathrm{\Phi }^\theta =\eta \mathrm{sin}f(\theta ),\mathrm{\Phi }^\phi =0,$$
(26)
where $`f(\theta )`$ is an arbitrary function of $`\theta `$ reflecting a gauge freedom. The regularity along the $`z`$axis imposes that $`\mathrm{sin}f(0)=\mathrm{sin}f(\pi )=0`$. Therefore we find two possible sets of boundary conditions: $`\mathrm{cos}f(0)=\mathrm{cos}f(\pi )=\pm 1`$ or $`\mathrm{cos}f(0)=\mathrm{cos}f(\pi )=\pm 1`$ (see also the discussion in ). We call the first set even boundary conditions, while the second set (for $`f(\pi )=f(0)+(2k+1)\pi )`$ will be refered to as odd boundary conditions.

For the generic Higgs field boundary conditions (26) and within the reduced ansatz (3), we find the asymptotic behavior in a flat space background

$`H_1`$ $`=`$ $`0,H_2=f^{}(\theta ),H_5=\alpha \mathrm{cos}f(\theta ),H_6=\alpha \mathrm{sin}f(\theta ),`$
$`F_{rt}^{(r)}`$ $`=`$ $`{\displaystyle \frac{Q_e\mathrm{cos}f(\theta )}{r^2}},F_{rt}^{(\theta )}={\displaystyle \frac{Q_e\mathrm{sin}f(\theta )}{r^2}}.`$ (27)
For even boundary conditions, the regularity and finite energy conditions impose the asymptotic expression of the gauge potentials $`H_3,H_4`$

$`H_3={\displaystyle \frac{\mathrm{cos}\theta \mathrm{cos}f(0)}{\mathrm{sin}\theta }}(\mathrm{cos}f(\theta )\mathrm{cos}f(0)),H_4={\displaystyle \frac{\mathrm{cos}\theta }{\mathrm{sin}\theta }}\mathrm{sin}f(\theta )\mathrm{cos}f(0),`$ (28)
while for the second type of boundary conditions we find

$`H_3={\displaystyle \frac{\mathrm{cos}f(0)}{\mathrm{sin}\theta }}(\mathrm{cos}f(\theta )\mathrm{cos}\theta \mathrm{cos}f(0)),H_4={\displaystyle \frac{\mathrm{cos}f(0)\mathrm{sin}f(\theta )}{\mathrm{sin}\theta }}.`$ (29)
By using these expressions and the relations (8), (13) we find that the even boundary configurations have a magnetic charge $`Q_m=\pm n`$ and a vanishing total angular momentum. The configurations satisfying odd-boundary conditions present a complementary picture: they posses a zero net magnetic charge and a nonzero angular momentum proportional to the electric charge $`J=\pm 4\pi nQ_e`$. Thus, as conjectured in , it seems that the vanishing of the total angular momentum is related to the presence of a net magnetic charge.

At this stage, we remark that the ansatz (3) is invariant under a gauge transformation $`U=\mathrm{exp}\{i\mathrm{\Gamma }(r,\theta )\tau _\phi ^n/2\}`$ , and the Higgs field transforms as

$`\mathrm{\Phi }^r`$ $``$ $`\mathrm{cos}\mathrm{\Gamma }\mathrm{\Phi }^r+\mathrm{sin}\mathrm{\Gamma }\mathrm{\Phi }^\theta ,`$
$`\mathrm{\Phi }^\theta `$ $``$ $`\mathrm{sin}\mathrm{\Gamma }\mathrm{\Phi }^r+\mathrm{cos}\mathrm{\Gamma }\mathrm{\Phi }^\theta ,`$ (30)
which implies $`ff\mathrm{\Gamma }(\theta )`$ in (26). We may use this relation to set $`f(\theta )=k\theta `$ in the Higgs field boundary conditions (with $`k`$ an integer number).

A number of well-known configurations (in a flat space background) are obtained for this choice of $`f`$. The dyon configurations discussed above correspond to $`k=0`$. For $`k=1`$, we have as $`r\mathrm{}`$

$`\mathrm{\Phi }^r=\eta \mathrm{cos}\theta ,\mathrm{\Phi }^\theta =\eta \mathrm{sin}\theta ,\mathrm{\Phi }^\phi =0,`$ (31)
which implies the asymptotic behavior, within the reduced ansatz (3)

$`H_1`$ $`=`$ $`H_3=0,H_2=H_4=1,H_5=\alpha \mathrm{cos}\theta ,H_6=\alpha \mathrm{sin}\theta ,`$
$`W^{(r)}`$ $`=`$ $`n\mathrm{cos}\theta ,W^{(\theta )}=n\mathrm{sin}\theta ,W^{(\phi )}=0,`$ (32)
$`F_{rt}^{(r)}`$ $`=`$ $`{\displaystyle \frac{Q_e\mathrm{cos}\theta }{r^2}},F_{rt}^{(\theta )}={\displaystyle \frac{Q_e\mathrm{sin}\theta }{r^2}}.`$
This corresponds to a monopole-antimonopole configuration, with a vanishing net magnetic charge but a nonzero electric charge. This regular solution is discussed in  for $`n=1,H_5=H_6=0`$ and a different basis $`\tau _a^n`$ in spherical coordinates. The authors of these papers use a different combination of $`\tau _\rho ^n`$ and $`\tau _z`$, more suitable for these boundary conditions ($`\tau _r^n=\mathrm{sin}2\theta \tau _\rho ^n+\mathrm{cos}2\theta \tau _z,\tau _\theta ^n=\mathrm{cos}2\theta \tau _\rho ^n\mathrm{sin}2\theta \tau _z`$). In  it was found that the total angular momentum of this solution is proportional to the electric charge.

For $`k=2`$ we find a configuration consisting in two monopoles and one antimonopole (we consider only the case $`n=1`$). This configuration is sustained by the boundary conditions at infinity

$`\mathrm{\Phi }^r`$ $`=`$ $`\eta \mathrm{cos}2\theta ,\mathrm{\Phi }^\theta =\eta \mathrm{sin}2\theta ,\mathrm{\Phi }^\phi =0,`$ (33)
which imply asymptotically

$`W^{(r)}`$ $`=`$ $`\mathrm{cos}\theta \mathrm{cos}2\theta ,W^{(\theta )}=\mathrm{cos}\theta \mathrm{sin}2\theta ,W^{(\phi )}=0,`$
$`F_{rt}^{(r)}`$ $`=`$ $`{\displaystyle \frac{Q_e\mathrm{cos}2\theta }{r^2}},F_{rt}^{(\theta )}={\displaystyle \frac{Q_e\mathrm{sin}2\theta }{r^2}},`$ (34)
(corresponding to $`H_1=0,H_2=2,H_3=\mathrm{sin}2\theta ,H_4=\mathrm{cos}2\theta 1,H_5=\alpha \mathrm{cos}2\theta ,H_6=\alpha \mathrm{sin}2\theta `$ for the ansatz (3)). One can verify that this configuration has a unit magnetic charge and a zero total angular momentum.

One may speculate on the existence of solutions for an arbitrary integer $`k`$, representing composed configurations, containing both monopoles and anti-monopoles. It appears that further analysis of such problems will require at least some use of numerical techniques. The role played by the Higgs field topology and the extra-symmetries of the ansatz needs further study also.

The inclusion of a negative cosmological constant in a YMH theory leads to the same general picture. Numerical results have been presented showing the existence of dyon solutions in fixed AdS spacetime. However, similar to the flat space counterparts, these dyon solutions present also a vanishing total angular momentum. We can consider of course more general configurations containing both monopoles and antimonopoles, by using the set of boundary conditions at infinity derived above. In this case, it is possible to write a similar relation between the electric charge and the angular momentum.

We expect that the inclusion of gravity in this analysis will not change these conclusions, since the gravitational field does not affect the behavior at infinity of the YMH system.

Although the obtained results are in agreement with the perturbative approach, we still are confronted by the problem of the physical mechanism preventing a configuration with net magnetic charge from rotating. A more general proof, beyond a specific ansatz (possibly based on topological arguments), is clearly necessary.
Acknowledgement
One of the authors (E.R.) acknowledges an important discussion with Erick Wöhnert.

This work was performed in the context of the Graduiertenkolleg of the Deutsche Forschungsgemeinschaft (DFG): Nichtlineare Differentialgleichungen: Modellierung,Theorie, Numerik, Visualisierung.

Figure Captions
Figure 1:
Typical spherically symmetric dyon solutions for the same value of Higgs self-coupling constant $`\lambda =1`$ and different values of $`\mathrm{\Lambda },\alpha `$. The total energy $`E`$ (in units $`4\pi \eta `$) and the electric charge $`Q_e`$ are also marked.
Figure 2:
The energy per topological charge $`E/n`$ in units $`4\pi \eta `$ (Figure 2a) and the electric charge per topological charge $`Q_e/n`$ (Figure 2b) are plotted as a function of the cosmological constant for dyon solutions with $`\alpha =0.2`$ and $`\lambda =0.5`$.
Figure 3:
The energy density $`T_t^t`$ (Figure 3a), surfaces of constant energy density (Figure 3b), and $`T_\phi ^t`$ component of the energy-momentum tensor (in units of $`4\pi \eta `$) (Figure 3c) for a dyon solution with topological charge $`n=2`$, total energy $`E=5.522`$, electric charge $`Q_e=0.086`$ and a vanishing Higgs self-coupling, as a function of the dimensionless compactified coordinates $`\rho =\overline{r}\mathrm{sin}\theta `$ and $`z=\overline{r}\mathrm{cos}\theta `$. Here $`\overline{x}=r/(1+r)`$ is the coordinate used in the numerical calculations. In Figure 3b, the inner surface correspond to a value of the energy density $`ϵ=4`$, while the value for the outer surface is 9.5. The cosmological constant is $`\mathrm{\Lambda }=6`$.
Figure 4:
The $`T_\phi ^t`$ component of the energy-momentum tensor (in units $`4\pi \eta `$) at $`\theta =\pi /4`$ for dyon solutions with $`\alpha =0.2,\lambda =0.5`$, is shown for several values of the cosmological constant. The topological charge of these solutions is $`n=3`$.

Figure 1.

Figure 2a.

Figure 2b.

Figure 3a.

Figure 3b.

Figure 3c.

Figure 4.


text:
# Folded solitary waves and Foldons in 2+1 dimensions

## Abstract.

A general type of localized excitations, folded solitary waves and foldons, are defined and studied both analytically and graphically. The folded solitary waves and foldons may be “folded” in quite complicated ways and possess quite rich structures and abundant interaction properties. The folded phenomenon is quite universal in the real natural world. The folded solitary waves and foldons are obtained from a quite universal formula and the universal formula is valid for some quite universal (2+1)-dimensional physical models. The “universal” formula is also extended to a more general form with many more independent arbitrary functions.

*PACS.*05.45.Yv, 02.30.Jr, 02.30.Ik.

In the study of nonlinear science, soliton theory plays a very important role and has been applied in almost all the natural sciences especially in all the physical branches such as the fluid physics, condense matter physics, biophysics, plasma physics, nonlinear optics, quantum field theory and particle physics etc.. Almost all the previous studies of soliton theory especially in high dimensions are restricted in single valued situations. However, the real natural phenomena are very complicated. In various cases, it is even impossible to describe the natural phenomena by single valued functions. For instance, in the real natural world, there exist very complicated folded phenomena such as the folded protein , folded brain and skin surfaces and many many other kinds of folded biologic systems. The simplest multi-valued (folded) waves may be the bubbles on (or under) a fluid surface. Various kinds of ocean waves are really folded waves also.

To study the complicated folded natural phenomena is very difficult. Similar to the single valued cases, the first important question we should ask is: Are there any stable multi-valued (folded) localized excitations? For convenience later, we define the multi-valued localized excitations folded solitary waves. Furthermore, if the interactions among the folded solitary waves are completely elastic, we call them foldons. In (1+1)-dimensional case, the simplest foldons are so-called loop solitons which can be found in many (1+1)-dimensional integrable systems and have been applied in some possible physical fields like the string interaction with external field, quantum field theory and particle physics. However, in our knowledge, there is no study at all for the possible higher dimensional foldons. That is the main topic of this letter.

In (1+1)-dimensions, it has been proven that when a physical model can be expressed by a partial differential equation (PDE), then under some suitable approximations, one can always find nonlinear Schördinger type equations that is why the (1+1)-dimensional soliton theory can be successively applied in almost all the physical branches. Similarly, if (2+1)-dimensional physical model can be expressed by a PDE, then under some suitable approximations, one can find the Davey-Stewartson (DS) type equations.

Recently, it is found that a quite “universal” formula

$`U{\displaystyle \frac{2\mathrm{\Delta }q_yp_x}{(a_0+a_1p+a_2q+a_3pq)^2}},\mathrm{\Delta }a_0a_3a_1a_2,`$ (1)
is valid for suitable fields or potentials of various (2+1)-dimensional physically interesting integrable models including the DS equation, the dispersive long wave equation (DLWE), the Broer-Kaup-Kupershmidt (BKK), and so on. In Eq. (1), $`pp(x,t)`$ is an arbitrary function of $`\{x,t\}`$, $`qq(y,t)`$ may be either an arbitrary function for some kinds of models such as the DS system or an arbitrary solution of a Riccati equation for some others, say, the DLWE while $`a_0,a_1,a_2`$ and $`a_3`$ are taken as constants. One of the most important results obtained from (1) is that for all the models mentioned above there are quite rich localized excitations such as the solitoffs, dromions, lumps, breathers, instantons, ring solitons, peakons, compactons, localized chaotic and fractal patterns and so on.

Now the natural and important question is: Can we find some types of (2+1)-dimensional foldons from the “universal” formula? Fortunately, the answer is obviously positive because of the arbitrariness of the functions $`p`$ and $`q`$ included in the universal formula. To find more general types of foldons, we extend the universal formula to a more general form for the DLWE before to give out concrete foldons. For the (2+1)-dimensional DLWE

$$u_{yt}+\eta _{xx}+u_xu_y+uu_{xy}=0,\eta _t+u_x+\eta u_x+u\eta _x+u_{xxy}=0,$$
(2)
(1) is valid for the field $`v\eta 1`$ and the related exact solution reads

$$v\eta 1=U,u=\pm \frac{2p_x(a_1+a_3q)}{a_0+a_1p+a_2q+a_3pq}+u_0,$$
(3)
with $`p`$ being an arbitrary functions of $`\{x,t\}`$, $`q=q(y,t)`$ being an arbitrary solution of the Riccati equation, $`q_ta_0c_0(a_1c_1+a_2c_0a_0c_2)q(a_3c_1a_2c_2)q^2=0`$, and $`u_0=p_x^1[p_t\pm p_{xx}a_0c_1(a_1c_1+a_2c_0+a_0c_2)p(a_1c_2+a_3c_0)p^2].`$ The (1+1)-dimensional DLWE ($`y=x`$ of (2)) is also called the classical Boussinesq equation. There exist a large number of papers to discuss the possible applications and exact solutions of the (1+1)-dimensional DLWE. Various interesting properties of the (2+1)-dimensional DLWE have been studied by many authors.

To extend the universal formula to a more general form, we use the multi-linear variable separation approach again. For the DLWE (2), the transformation $`\{u=u_0\pm 2\frac{f_x}{f},\eta =2\frac{f_{xy}}{f}2\frac{f_xf_y}{f^2}1\}`$ with $`u_0u_0(x,t)`$ degenerates two equations of (2) to a trilinear form

$`[f_{xxxy}\pm (f_{xyt}+u_{0x}f_{xy}+u_0f_{xxy})]f^2[(f_{xxx}\pm u_0f_{xx}\pm f_{xt})f_y\pm f_xf_yu_{0x}`$
$`+(f_{xxy}\pm 2u_0f_{xy}\pm f_{yt})f_x+f_{xy}f_{xx}\pm f_tf_{xy}]f+2f_xf_y(f_{xx}\pm u_0f_x\pm f_t)=0.`$ (4)
To solve the trilinear equation (Folded solitary waves and Foldons in 2+1 dimensions), one has to use some prior ansatz. The variable separation solution (3) is resulted from the ansatz $`f=a_0+a_1p+a_2q+a_3pq`$.

For the DLWE (2), there are two sets of infinitely many symmetries and every symmetry possesses an arbitrary function of $`t`$ or $`y`$. That means infinitely many arbitrary functions of $`y`$ and $`t`$ can be entered into the solutions of (2). So it is possible to extend the solution (3) to a more general one with more arbitrary functions. After finishing the detailed calculations, we find that the following variable separation ansatz,

$`f=q_0+{\displaystyle \underset{i=1}{\overset{N}{}}}p_iq_i,`$ (5)
where $`\{q_i,i=0,1,2,\mathrm{},N\}`$, and $`\{p_i,i=1,2,\mathrm{},N\}`$ are functions of $`\{y,t\}`$ and $`\{x,t\}`$ respectively, solves the trilinear equation (Folded solitary waves and Foldons in 2+1 dimensions) under the conditions

$`q_{it}={\displaystyle \underset{j=0}{\overset{N}{}}}(c_{i,j}+q_iC_j)q_j,p_{it}=(c_{00}u_0_x_x^2)p_ic_{0i}+{\displaystyle \underset{j=1}{\overset{N}{}}}(c_{j0}p_ic_{ji})p_j,`$ (6)
for all $`i=0,1,2,\mathrm{},N`$ and $`\{c_{ij},i,j=0,1,2,\mathrm{},N,C_j,j=1,2,\mathrm{},N\}`$ being arbitrary functions of $`t`$. Obviously, the general ansatz (5) will return back to the known one for $`N=2,q_0=a_0,q_1=a_1,q_2=a_2q,q_3=a_3q,p_1=p_3=p,p_2=1`$.

The corresponding solution for the field $`v=\eta 1`$ now reads

$$v=\frac{2_{i=1}^Np_{ix}q_{iy}}{q_0+_{i=1}^Np_iq_i}+\frac{2_{i=1}^Np_{ix}q_i\left(q_{0y}+_{j=1}^Np_jq_{jy}\right)}{(q_0+_{i=1}^Np_iq_i)^2}U_E$$
(7)
while the quantity $`u`$ is given by $`u=\pm (2_{i=1}^Np_{ix}q_i)/(q_0+_{i=1}^Np_iq_i)+u_0.`$ It is clear that in addition to a (1+1)-dimensional arbitrary function of $`\{x,t\}`$ (one of $`u_0`$ and $`p_i`$), $`(N+1)(N+2)1`$ arbitrary functions of $`t`$, $`c_{ij},C_j`$ have been included in the general solution. Furthermore, various arbitrary functions of $`y`$ and $`\{x,t\}`$ will also be included in (7) after the coupled systems of (6) being solved. Because of the complexity, we have to leave these problems for our future studies. Here we just write down the simplest nontrivial case for the later uses. It is quite trivial that, when we fix $`N=1`$, $`c_{ij}=C_i=0,p_1=p`$ and $`q_0a_0+q_0`$, the formula (7) is simplified to

$$v=\frac{2p_x\left(q_1q_{0y}(a_0+q_0)q_{1y}\right)}{(a_0+q_0+pq_1)^2}$$
(8)
with $`q_0`$ and $`q_1`$ being arbitrary functions of $`y`$ and $`p`$ being an arbitrary function of $`\{x,t\}`$. We can prove that the simplified quantity expressed by the right hand side of (8) does work for all the known models that allow the universal formula (1). It is interesting also that though we have not yet proven the validity of (7) for all the models mentioned in , we do have proven that the extended form (7) is really valid at least for some of them such as the BKK system, the (2+1)-dimensional Burgers system and the (1+1)-component AKNS system. For instance, the quantities $`g2G=U_E`$ and $`H=(\mathrm{ln}f)_x+u_0/2`$ with the same $`f,p_i`$ and $`q_i`$ determined by (5) and (6) solve the BKK system $`\{H_{ty}H_{xxy}+2(HH_x)_y+2G_{xx}=0,G_t+G_{xx}+2(HG)_x=0\}`$.

As pointed out in , because of the arbitrariness of the functions appeared in the universal formula and its extended forms, various class of localized excitations can be constructed. Now we concentrate on how to find some types of folded solitary waves and foldons from the universal formula and its extended forms. First of all, we write a localized function, $`px`$, in the form

$`px{\displaystyle \underset{j=1}{\overset{M}{}}}f_j(\xi +v_jt),x=\xi +{\displaystyle \underset{i=1}{\overset{M}{}}}g_j(\xi +v_jt),`$ (9)
where $`v_1<v_2<\mathrm{}v_M`$ are all arbitrary constants and $`\{f_j,g_j\},j`$ are all localized functions with the properties $`f_j(\pm \mathrm{})=0,g_i(\pm \mathrm{})=G_i^\pm =consts.`$ From the second equation of (9), we know that $`\xi `$ may be a multi-valued function in some possible regions of $`x`$ by selecting the functions $`g_j`$ suitably. So the function $`px`$ may be a multi-valued function of $`x`$ in these regions though it is a single valued function of $`\xi `$. It is also clear that $`px`$ is an interaction solution of $`M`$ localized excitations because of the property $`\xi |_x\mathrm{}\mathrm{}`$. Actually, most of the known multi-loop solutions are the special cases of (9). Now if we take all the arbitrary functions appeared in the universal formula (1) and/or its slightly general one (8) (and/or even its more general extended one (7)) possess the forms similar to (9), then we can get various types of folded solitary waves and/or foldons.

In Fig. 1, four typical special folded solitary waves are plotted for the quantity $`v`$ shown by (8) with $`p_x=px`$, $`p=^\xi p_xx_\xi d\xi `$ and the functions $`q_0`$ and $`q_1`$ being given in a similar way

$`q_{\delta y}={\displaystyle \underset{j=1}{\overset{M_\delta }{}}}Q_{\delta j}(\theta ),y=\theta +R(\theta ),q_\delta ={\displaystyle ^\theta }q_{\delta y}y_\theta d\theta ,(\delta =0,1).`$ (10)
In (10), $`Q_{\delta j}(\theta ),j`$ and $`R(\theta )`$ are localized functions of $`\theta `$. The more detailed function selections of the figures are directly given in the figure caption.

To construct the possible localized excitations with completely elastic interaction behaviors, especially, foldons, we should discuss the asymptotic behaviors of (8) and (9) when $`t\pm \mathrm{}`$. Assuming further that $`F_j(\omega _j)F_j(\xi +v_jt)f_j\mathrm{dx}|_{\omega _j\pm \mathrm{}}F_j^\pm `$, we can gaze at the $`j`$th excitation to see the interaction properties among the localized excitations expressed by (8) with (9) and (10). In other words, we can take $`\omega _j`$ invariant and then take $`t\mathrm{}`$ because $`q_0`$ and $`q_1`$ have been fixed as $`t`$-independent. The results read

$`v|_t\mathrm{}{\displaystyle \frac{2f_j(\omega _j)[q_1q_{0y}(a_0+q_0)q_y]}{[a_0+q_0+q_1(F_j(\omega _j)+\mathrm{\Omega }_j^{})]^2}},x|_t\mathrm{}\xi +\mathrm{\Gamma }_j^{}+g_j(\xi +v_jt),`$ (11)
where $`\mathrm{\Omega }_j^{}=_{i>j}F_i^{}+_{i<j}F_i^\pm `$ and $`\mathrm{\Gamma }_j^{}=_{i>j}G_i^{}+_{i<j}G_i^\pm `$. From the asymptotic results (11), we know that, (i) the $`j`$th localized excitation given by (8) with (9) and (10) is a travelling wave moving in the velocity $`v_j`$ along the negative ($`v_j>0`$) or positive ($`v_j<0`$) $`x`$-direction; (ii) the multi-valued properties (i.e., the structures) of the $`j`$th localized excitation is only determined by $`g_j`$ of (9); (iii) the shape of the $`j`$th excitation will be changed if $`\mathrm{\Omega }_j^+\mathrm{\Omega }_j^{}`$, however, it will preserve its shape if $`\mathrm{\Omega }_j^+=\mathrm{\Omega }_j^{}`$; (iv) the total phase shift for the $`j`$th excitation is $`\mathrm{\Gamma }_j^+\mathrm{\Gamma }_j^{}`$.

From the above discussions, it is seen that to construct some (2+1)-dimensional localized soliton type excitations with completely elastic interaction behaviors becomes an easy task. The only thing one has to do is to select suitable localized functions $`F_j`$ (and then $`f_j`$) with $`\mathrm{\Omega }_j^+=\mathrm{\Omega }_j^{}`$ in (9)! The localized soliton like solutions can still possess quite rich structures. Especially, various types of known single-valued localized excitations are also in the catalog due to the arbitrariness of the functions included in (9) and (10). For instance, if we select $`f_i`$ and $`Q_i`$ of (9) and (10) possesses the peakon shape solutions and/or the compacton shape solutions with the property $`\mathrm{\Omega }_i^+=\mathrm{\Omega }_i^{}`$, then we can find some (2+1)-dimensional peakon type and/or compacton type solutions with completely elastic interaction properties.

Figs. 2a–2c are plotted to show the possible existence of foldons which are given by (8) with $`q_1=1`$, (9) and (10) for $`\delta =0`$ and the concrete function selections are also given in the figure caption. From Figs. 2a and 2c, we can see that the interaction of two foldons is completely elastic. Because one of the velocities of foldons has been selected as zero, it can also be seen that there are phase shifts for two foldons. Especially, before the interaction, the static foldon (the large one) is located at $`x=1.5`$ and after the interaction, the large foldon is shifted to $`x=1.5`$.

In summary, starting from the “universal” formula and/or its extended form, various kinds of localized excitations with and without completely elastic interaction behaviors can be constructed easily by selecting the arbitrary functions suitably according to the asymptotic results (11). Especially, a new kinds of localized excitations, folded solitary waves and foldons, are investigated both analytically and graphically. The foldons may be folded quite freely and complicatedly and then possess quite rich structures and interaction behaviors. The explicit phase shifts for all the localized excitations given by (8) with (9) and (10) have been given.

On the one hand, there are a large number of complicated “folded” and/or the multi-valued phenomena in the real natural world. On the other hand, there is no good analytical way to treat these kinds of complicated phenomena. This letter is only a beginning attempt, to find some types of possible stable multi-valued localized excitations, folded solitary waves and foldons, for some real physical models including some quite “universal” systems such as the well known DS equation and the DLWE.

Because the formula (8) is valid for various (2+1)-dimensional interesting models, which are divergently applied in many physical fields, we believe that the foldons may be useful in the future studies on the complicated “folded” natural world. The more about both the (extended) “universal” formula and the general (or special) foldons especially its real possible applications should be studied further.

The authors are in debt to thank the helpful discussions with Prof. C. Rogers, Dr. W. Schief, Prof. C-z Qu and Dr. H-c Hu. The work was supported by the National Outstanding Youth Foundation of China (No. 19925522), the National Natural Science Foundation of China (No. 90203001) and the Research Fund for the Doctoral Program of Higher Education of China (Grant. No. 2000024832) and ARC Large Grant A00103139 of the University of New South Wales University, Australia.