text:
# 1 Introduction: The “spinless-Salpeter” Hamiltonian

## 1 Introduction: The “spinless-Salpeter” Hamiltonian

We investigate the discrete spectrum—that is, the set of eigenvalues $`E_n\mathrm{}`$ below the onset of the essential spectrum, corresponding to bound states characterized by the radial quantum number $`n=1,2,3,\mathrm{}`$ and the angular-momentum quantum number $`\mathrm{}=0,1,2,\mathrm{}`$—of the semirelativistic “spinless-Salpeter” Hamiltonian for particles of mass $`m`$ and momentum $`𝐩,`$

$$H=\beta \sqrt{m^2+𝐩^2}+V(r),r|𝐱|,\beta >0,$$
(1)
where $`V(r)`$ is an attractive central potential in three spatial dimensions. As is evident from the presence of the relativistic kinetic-energy operator $`\beta \sqrt{m^2+𝐩^2}`$ (with $`\beta >0`$ allowing, e.g., for more than one particle), this Hamiltonian may be regarded as the generalization of the nonrelativistic Schrödinger Hamiltonian towards relativistic kinematics. The eigenvalue equation of the operator (1), called the “spinless Salpeter equation,” arises as a well-defined approximation to the Bethe–Salpeter formalism for the description of bound states within a relativistic quantum field theory. Our main goal is to discuss bounds on the energy levels $`E.`$

## 2 Envelope theory for energy bounds

The envelope theory, developed in Refs. , constructs bounds on the discrete eigenvalues $`E_n\mathrm{}`$ of a given operator $`H`$ by comparing the spectrum of $`H`$ with the spectrum of a suitably formulated “tangential problem,” for which some spectral information is known (or may be obtained with considerable more ease than for the original problem). More specifically, when applied to the (semirelativistic) Hamiltonian $`H`$ defined in Eq. (1) , the envelope theory compares $`H`$ with a corresponding “tangential Hamiltonian”

$$\stackrel{~}{H}=\beta \sqrt{m^2+𝐩^2}+vh(r),v>0,$$

involving some “basis potential” $`h(r).`$ Under the assumption that the interaction potential under consideration, $`V(r),`$ is a smooth transformation $`V(r)=g(h(r))`$ of the basis potential $`h(r),`$ where the function $`g(h)`$ has definite convexity, this comparison yields rigorous bounds on the discrete eigenvalues $`E_n\mathrm{}`$ of $`H.`$ More precisely, if the function $`g(h)`$ is convex, that is, $`g^{\prime \prime }(h)>0,`$ we obtain lower energy bounds; if the function $`g(h)`$ is concave, that is, $`g^{\prime \prime }(h)<0,`$ we obtain upper energy bounds. If the eigenvalues $`\stackrel{~}{E}`$ of $`\stackrel{~}{H}`$ are not known (analytically), the envelope theory may take advantage of the knowledge of suitable bounds on the eigenvalues, namely, of lower bounds on $`\stackrel{~}{E}_n\mathrm{}`$ for convex $`g(h)`$ or of upper bounds on $`\stackrel{~}{E}_n\mathrm{}`$ for concave $`g(h).`$

Suppressing the quantum numbers $`n,`$ $`\mathrm{}`$ identifying the bound state under consideration, we denote, for a given value of the coupling $`v`$ entering into the “tangential Hamiltonian” $`\stackrel{~}{H},`$ the eigenvalues of $`\stackrel{~}{H}`$ or the appropriate—in the above sense—bounds to the latter by $`e(v).`$ The strict bounds on the discrete eigenvalues $`E`$ of $`H`$ predicted within the framework of the envelope technique may then be summarized in form of the “principal envelope formula”

$$E\underset{v>0}{\mathrm{min}}[e(v)ve^{}(v)+g(e^{}(v))],$$
(2)
where the sign of approximate equality is used to recall that, for a definite convexity of $`g(h),`$ the envelope theory yields lower bounds for convex $`g(h)`$ and upper bounds for concave $`g(h).`$

For the tangential problem posed by the (semirelativistic) “tangential Hamiltonian” $`\stackrel{~}{H},`$ spectral information is, at present, available if the basis potential $`h(r)`$ is either the Coulomb potential, that is, $`h(r)=1/r,`$ or the harmonic-oscillator potential, that is, $`h(r)=r^2.`$ The energy bounds resulting in these two cases are discussed, in turn, in Sec. 4 and Sec. 5. The case of the interaction potential $`V(r)`$ in Eq. (1) being the harmonic-oscillator potential may be traced back to a nonrelativistic Schrödinger problem; this is studied for its own in Sec. 3.

## 3 Harmonic-oscillator potential

The “relativistic harmonic-oscillator problem,” defined by the special case of $`V(r)`$ being the harmonic-oscillator potential $`V(r)=vr^2,`$ may be analyzed by observing  that in momentum-space representation the semirelativistic Hamiltonian $`H=\beta \sqrt{m^2+𝐩^2}+vr^2`$ of Eq. (1) simplifies to the (nonrelativistic) Schrödinger operator $``$ given by

$$=v\mathrm{\Delta }+𝒱(r),𝒱(r)=\beta \sqrt{m^2+r^2}.$$
(3)
Thus the eigenvalue equation for $`H`$ reduces to an easier-to-treat nonrelativistic Schrödinger problem with an effective interaction potential $`𝒱(r)`$ which is reminiscent of the square root.

### 3.1 Upper bounds on harmonic-oscillator energy levels

The potential $`𝒱`$ is a concave transform of a harmonic-oscillator potential. Thus we may find upper bounds on the eigenvalues $`E_n\mathrm{}(v)`$ of $``$ (cf. Eq. (2.2) of Ref.  or Eq. (12) of Ref. ):

$$E_n\mathrm{}(v)\underset{r>0}{\mathrm{min}}\left[\beta \sqrt{m^2+\frac{P_n\mathrm{}^2(2)}{r^2}}+vr^2\right],P_n\mathrm{}(2)=2n+\mathrm{}\frac{1}{2}.$$
(4)
The (dimensionless) numbers $`P_n\mathrm{}(2)`$ introduced here are related to the eigenvalues $`_n\mathrm{}`$ of the Hamiltonian $`\stackrel{~}{H}=𝐩^2+vr^2`$ by $`_n\mathrm{}=2\sqrt{v}P_n\mathrm{}(2).`$

Upper bounds on the eigenvalues of self-adjoint operators bounded from below may also be found by combining the minimum–maximum principle  with appropriate operator inequalities. (For details, see, e.g., Sec. 3.4 of Ref. , Sec. 3.5 of Ref. , Sec. 2.4 of Ref. , or Appendix A of Ref. .) Exploiting the obvious positivity of $`(\sqrt{m^2+𝐩^2}\mu )^2,`$ where $`\mu `$ is an arbitrary real parameter (with the dimension of mass), we obtain a set of inequalities for the square-root operator $`\sqrt{m^2+𝐩^2}`$ entering into the kinetic energy of the Hamiltonian (1); cf. also Ref. :

$$\sqrt{m^2+𝐩^2}\frac{𝐩^2+m^2+\mu ^2}{2\mu }\mu >0.$$

(This operator inequality may also be found by constructing the tangent line $`a(𝐩_1^2)𝐩^2+b(𝐩_1^2)`$ to the square root $`\sqrt{m^2+𝐩^2}`$ at the point of contact $`𝐩_1^2=\mu ^2m^2;`$ see Sec. 4 of Ref. .) Amending this result by the interaction potential $`V(r),`$ one finds that the spinless-Salpeter Hamiltonian (1) satisfies

$$HH_{\mathrm{NR}}(\mu )\beta \frac{𝐩^2+m^2+\mu ^2}{2\mu }+V(r)\mu >0.$$

The min–max principle then tells us that the discrete eigenvalues $`E_n\mathrm{}`$ of the semirelativistic Hamiltonian $`H`$ of Eq. (1) are bounded from above by the corresponding discrete eigenvalues $`E_{\mathrm{NR},n\mathrm{}}(\mu )`$ of the nonrelativistic Hamiltonian $`H_{\mathrm{NR}}(\mu )`$ and consequently also by the minimum $`\overline{E}_n\mathrm{}`$ with respect to $`\mu `$ of all these upper bounds :

$$E_n\mathrm{}\overline{E}_n\mathrm{}\underset{\mu >0}{\mathrm{min}}E_{\mathrm{NR},n\mathrm{}}(\mu ).$$

Applying these general findings to the case of the harmonic-oscillator potential, one ends up with the following upper bounds for all energy levels of the relativistic harmonic-oscillator problem (see Appendix A of Ref. ; for the ease of comparison, we express here the bounds of Ref.  in terms of $`P_n\mathrm{}(2)`$ rather than in terms of $`E_{\mathrm{NR},n\mathrm{}}(\mu )`$):

$$E_n\mathrm{}(v)\overline{E}_n\mathrm{}\underset{\mu >0}{\mathrm{min}}\left[\beta \frac{m^2+\mu ^2}{2\mu }+\sqrt{\frac{2\beta v}{\mu }}P_n\mathrm{}(2)\right],P_n\mathrm{}(2)=2n+\mathrm{}\frac{1}{2}.$$
(5)
Despite their seemingly different form, the upper bounds in Eqs. (4) and (5) are, in fact, identical, in the sense that the two functions on the right-hand side of these inequalities have their minimum at the same critical point and take the same miminal value at this point.

In order to demonstrate this equivalence, let us start from the “min–max” expression (5) for these upper bounds. For simplicity of notation, we suppress for the moment the quantum numbers $`n,`$ $`\mathrm{}`$ which identify the bound states under consideration. The upper bound $`\overline{E}`$ then reads

$$\overline{E}=\underset{\mu >0}{\mathrm{min}}\left[\beta \frac{m^2+\mu ^2}{2\mu }+\sqrt{\frac{\beta v}{2\mu }}2P\right].$$

We introduce the abbreviation $`q=p^2,`$ define the function $`f(q)=\sqrt{m^2+q},`$ which satisfies

$$f^{}(q)=\frac{1}{2f(q)},f^2(q)=m^2+q,$$

and identify the (positive but otherwise arbitrary) parameter $`\mu `$ according to $`\mu =f(q).`$ Next, we introduce a function $`F`$ by the definition $`F(v)=2P\sqrt{v}.`$ A couple of algebraic steps then recasts the upper bound $`\overline{E}`$ into the form

$$\overline{E}=\underset{q>0}{\mathrm{min}}[F(\beta vf^{})+\beta (fqf^{})].$$

The critical points of this latter expression are obtained by differentiation with respect to $`q.`$ Since $`f^{\prime \prime }(q)0,`$ this yields

$$\widehat{q}=vF^{}(\beta vf^{}(\widehat{q}))$$

and

$$\overline{E}=F(\beta vf^{}(\widehat{q}))\beta vf^{}(\widehat{q})F^{}(\beta vf^{}(\widehat{q}))+\beta f(\widehat{q}).$$

Let us define a new problem, related to the above one by a Legendre transformation , which has the same minimum and the same critical point, but the function to be minimized is different. Let $`u(q)=\beta vf^{}(q)>0`$ and consider the energy expression

$$\widehat{E}=\underset{u>0}{\mathrm{min}}[F(u)uF^{}(u)+\beta f(vF^{}(u))].$$

Since $`F^{\prime \prime }(u)0,`$ we find that the critical point is given by $`\widehat{u}=\beta vf^{}(vF^{}(\widehat{u})).`$ But, given the definition of $`u,`$ this equation is equivalent to the earlier critical equation $`\widehat{q}=vF^{}(\beta vf^{}(\widehat{q})).`$ Meanwhile, the critical value $`\widehat{E}`$ is also the same as the critical value $`\overline{E}`$ because we have

$$\widehat{E}=F(\beta vf^{}(\widehat{q}))\beta vf^{}(\widehat{q})F^{}(\beta vf^{}(\widehat{q}))+\beta f(\widehat{q})=\overline{E}.$$

Having shown that $`\overline{E}`$ yields the same results as $`\widehat{E},`$ we now transform $`\widehat{E}`$ into the form (4). First of all, we note that $`F(u)uF^{}(u)=P\sqrt{u}.`$ Next, we define a new dummy minimization variable $`r>0`$ by

$$\sqrt{u}=\frac{vr^2}{P}.$$

From this, it follows that

$$F^{}(u)=\frac{P}{\sqrt{u}}=\frac{P^2}{vr^2}.$$

Finally, remembering the definition $`f(q)=\sqrt{m^2+q}`$ of $`f(q),`$ we may transform $`\widehat{E}`$ into the expression on the right-hand side of the inequality (4); this establishes the equivalence of the upper bounds (4) and (5).

### 3.2 Lower bounds on harmonic-oscillator energy levels

The potential $`𝒱`$ in Eq. (3) is a convex transform of a linear potential. Therefore we may find lower bounds on the eigenvalues $`E_n\mathrm{}(v)`$ of $``$ (cf. Eq. (2.2) of Ref.  or Eq. (12) of Ref. ):

$$E_n\mathrm{}(v)\underset{r>0}{\mathrm{min}}\left[\beta \sqrt{m^2+\frac{P_n\mathrm{}^2(1)}{r^2}}+vr^2\right],$$
(6)
where the (dimensionless) numbers $`P_n\mathrm{}(1)`$ introduced here are related to the eigenvalues $`_n\mathrm{}`$ of the Hamiltonian $`\stackrel{~}{H}=𝐩^2+vr`$ by

$$_n\mathrm{}=3\left[\frac{vP_n\mathrm{}(1)}{2}\right]^{2/3};$$

for the lowest states ($`n=1,2,\mathrm{},5,`$ $`\mathrm{}=0,1,\mathrm{},5`$) these numbers may be found in Table 1.

Lower bounds on the eigenvalues of (self-adjoint) Hamiltonians involving the relativistic kinetic energy $`\sqrt{m^2+𝐩^2}`$ may also be found by taking into account the obvious positivity of $`[𝐩^2\xi ^2m^2(1\xi ^2)]^2`$ for an arbitrary real parameter $`\xi ,`$ which may be converted into the set of operator inequalities

$$\sqrt{m^2+𝐩^2}|𝐩|\sqrt{1\xi ^2}+m\xi ,0\xi 1.$$
(7)
Applying these general findings to the semirelativistic Hamiltonian $`H`$ of Eq. (1), with $`V(r)`$ being the harmonic-oscillator potential, we obtain lower bounds on all discrete energy levels $`E_n\mathrm{}`$ of the relativistic harmonic-oscillator problem (for details, see Appendix A of Ref. ):

$$E_n\mathrm{}(v)\underset{¯}{E}_n\mathrm{}\underset{0\xi 1}{\mathrm{max}}[\beta \xi m+\beta ^{2/3}(1\xi ^2)^{1/3}E_n\mathrm{}(v;m=0)],$$
(8)
where $`E_n\mathrm{}(v;m=0)`$ denote the corresponding discrete eigenvalues of the operator $`|𝐩|+vr^2,`$ i.e., of the Hamiltonian (1) for $`\beta =1`$ and a vanishing particle mass $`m.`$ For states with $`\mathrm{}=0,`$ for instance, these “zero-mass” energy levels are given, in terms of the zeros $`z_n`$ of the Airy function $`\mathrm{Ai}(z)`$  ($`z_n=2.33810,\mathrm{\hspace{0.17em}4.08794},\mathrm{\hspace{0.17em}5.52055},\mathrm{}`$), by $`E_{n0}(v;m=0)=v^{1/3}z_n`$ .

As in Sec. 3.1 we observe that the lower bounds in Eqs. (6) and (8) appear, at first sight, to be rather different. A closer inspection, however, reveals that these bounds are identical.

In order to present an analytical proof of the equivalence of the lower bounds (6) and (8) on the energy eigenvalues of the relativistic harmonic-oscillator problem, we recall that the envelope lower bounds in (6) are based (in the momentum-space representation of $`H`$) on the observation that, since its second derivative with respect to $`|𝐩|`$ is positive, the square-root operator $`\sqrt{m^2+𝐩^2}`$ of the relativistic kinetic energy is bounded from below by tangent lines $`a(t)+b(t)|𝐩|;`$ this may be expressed by the geometrical inequality $`\sqrt{m^2+𝐩^2}a(t)+b(t)|𝐩|,`$ where, by elementary calculus,

$$a(t)=\frac{m^2}{\sqrt{m^2+t^2}},b(t)=\frac{t}{\sqrt{m^2+t^2}},$$

and $`|𝐩|=t`$ is the point of contact between the square root of the relativistic kinetic energy and its straight-line lower approximation. The equivalence of the above two expressions for the harmonic-oscillator lower bound is then shown by a simple change of variable from $`t`$ to $`\xi `$ given explicitly by

$$\xi =\frac{a(t)}{m},$$

which, after a little algebra, leads to $`b(t)=\sqrt{1\xi ^2}.`$ For these coefficient functions $`a(t)`$ and $`b(t)`$ the above geometrical inequality becomes exactly the lower bound (7) on the relativistic kinetic energy. Consequently, with the basic inequalities being identical, the resulting lower energy bounds have to be identical too.

## 4 Convex transform of Coulomb potential

For the “relativistic Coulomb problem” posed by the spinless-Salpeter Hamiltonian (1) with $`V(r)`$ being the Coulomb potential $`V(r)=v/r,`$ for coupling constants $`v`$ smaller than $`\beta v_\mathrm{c},`$ that is, $`v<\beta v_\mathrm{c},`$ where the critical value $`v_\mathrm{c}`$ of the Coulomb coupling constant $`v`$ is given by

$$v_\mathrm{c}=\frac{2}{\pi },$$

a lower bound to the ground-state energy eigenvalue $`E_0`$ of $`H`$ (the bottom of the spectrum of $`H`$) has been derived by Herbst :

$$E_0\beta m\sqrt{1\left(\frac{\sigma v}{\beta }\right)^2},\sigma \frac{\pi }{2}.$$
(9)
For some part of the allowed range of the coupling $`v,`$ this lower bound has been improved by Martin and Roy :

$$E_0\beta m\sqrt{\frac{1}{2}\left[1+\sqrt{1\left(\frac{2v}{\beta }\right)^2}\right]},v<\frac{\beta }{2}.$$
(10)
In Ref. , we recast this Martin–Roy bound into a class of lower bounds on the spectrum of $`H`$ all of which are of the form of Herbst’s bound (9) \[but, of course, slightly weaker than the Martin–Roy bound (10)\]; the new lower bounds are parametrized by a real parameter $`\sigma 1`$:

$$E_0\beta m\sqrt{1\left(\frac{\sigma v}{\beta }\right)^2},v\beta \frac{\sqrt{\sigma ^21}}{\sigma ^2}<\frac{\beta }{2}.$$
(11)
The constraint on the Coulomb coupling $`v`$ in Eq. (11) automatically guarantees $`v\beta /\sigma ;`$ this latter constraint is a necessary and sufficient condition for the reality of the lower energy bounds in Eq. (11).

From the lower bounds (9) or (11) on the spectrum of the relativistic Coulomb problem, for the class of potentials $`V(r)`$ which are convex transforms $`V(r)=g(h(r)),`$ $`g^{\prime \prime }>0,`$ of the Coulomb potential $`h(r)=1/r`$ an (“absolute”) envelope lower bound on the ground-state eigenvalue $`E_0`$ (and, hence, on the entire spectrum) of the spinless-Salpeter Hamiltonian (1) may be found (cf. Eq. (16) of Ref. ):

$$E_0\underset{r>0}{\mathrm{min}}\left[\beta \sqrt{m^2+\frac{P^2}{r^2}}+V(r)\right],v<\beta v_P.$$
(12)
Here we introduced a parameter $`P`$ by defining $`P=1/\sigma .`$ Accordingly, the boundary value $`v_P`$ of the Coulomb coupling constant $`v`$ is given, when arising from demanding the Hamiltonian (1) to be bounded from below, by the critical coupling $`v_\mathrm{c},`$ i.e., $`v_P=v_\mathrm{c},`$ or, when arising from the region of validity of our new family of Coulomb lower bounds (11), by

$$v_P=P\sqrt{1P^2}<\frac{1}{2}.$$
(13)
It is easy to convince oneself that the Coulomb lower energy bound on the right-hand side of the inequality (12) is a monotone increasing function $`\underset{¯}{E}(P)`$ of the parameter $`P.`$ To this end, we rewrite the inequality (12) in the form

$$E_0\underset{¯}{E}(P)=\underset{r>0}{\mathrm{min}}F(r,P),$$

with

$$F(r,P)\beta \sqrt{m^2+\frac{P^2}{r^2}}+V(r).$$

Remembering that $`V(r)`$ is the transform $`V(r)=g(h(r))`$ of $`h(r)=1/r,`$ and assuming that $`g^{}(h)>0,`$ the critical point $`\widehat{r}(P)`$ (corresponding to the minimum of $`F(r,P)`$ for a given $`P`$) may be found, from

$$\frac{F(r,P)}{r}|_{r=\widehat{r}(P)}=0,$$

by solving

$$\frac{\beta P^2}{\sqrt{m^2\widehat{r}^2+P^2}}=g^{}(h(\widehat{r})).$$

The substitution of the solution $`\widehat{r}(P)`$ of the latter equation in $`F(r,P)`$ gives the lower bound $`\underset{¯}{E}(P)=F(\widehat{r}(P),P);`$ here we are interested in the total derivative of $`\underset{¯}{E}(P)`$ with respect to $`P.`$ Formally, this derivative would be given by

$$\frac{\mathrm{d}\underset{¯}{E}(P)}{\mathrm{d}P}=\frac{F(\widehat{r}(P),P)}{P}+\frac{F(r,P)}{r}|_{r=\widehat{r}(P)}\frac{\mathrm{d}\widehat{r}(P)}{\mathrm{d}P}.$$

At the critical point, however, the first factor of the second term vanishes and one is left with

$$\frac{\mathrm{d}\underset{¯}{E}(P)}{\mathrm{d}P}=\frac{F(\widehat{r}(P),P)}{P}=\frac{\beta P}{\widehat{r}\sqrt{m^2\widehat{r}^2+P^2}}>0.$$

From this we conclude that the lower bound $`\underset{¯}{E}(P)`$ is a monotone increasing function of $`P.`$

The existence of the upper bound $`v<\beta v_P<\beta /2`$ on the Coulomb coupling constant $`v`$ is reflected in a corresponding constraint on the parameters of the interaction potential $`V(r).`$ The envelope theory constructs a “tangential potential” to $`V(r)`$ at a point of contact $`r=t.`$ For a potential $`V(r)`$ which is a convex transform $`V(r)=g(h(r))`$ of the basis potential $`h(r),`$ this tangential potential forms a lower bound to $`V(r);`$ for $`h(r)`$ being the Coulomb potential, $`h(r)=1/r,`$ this lower bound is a Coulomb potential with the effective coupling constant $`v(t)=t^2V^{}(t),`$ shifted by the constant $`V(t)+tV^{}(t)`$:

$$V(r)\stackrel{~}{V}(r)\frac{t^2V^{}(t)}{r}+V(t)+tV^{}(t).$$

The relativistic Coulomb problem is posed by the “Coulombic” Hamiltonian

$$H_\mathrm{C}(v)=\beta \sqrt{m^2+𝐩^2}\frac{v}{r},v<\beta v_\mathrm{c},$$

which is nothing else but the tangential Hamiltonian $`\stackrel{~}{H}`$ with $`h(r)=1/r.`$ In the notation of Sec. 2, the eigenvalues $`E_\mathrm{C}(v)`$ of $`H_\mathrm{C}(v)`$ are bounded from below according to $`E_\mathrm{C}(v)e(v),`$ with $`e(v)`$ given by Herbst’s bound (9), our new bound (11), or the Martin–Roy bound (10). Expressed in terms of the Coulombic Hamiltonian $`H_\mathrm{C}(v),`$ the semirelativistic Hamiltonian (1) thus satisfies the inequality $`HH_\mathrm{C}(v(t))+V(t)+tV^{}(t);`$ its eigenvalues $`E`$ are bounded from below according to

$$E_0E_{\mathrm{C},0}(v(t))+V(t)+tV^{}(t)e(v(t))+V(t)+tV^{}(t).$$

The optimized lower bound on the discrete spectrum of $`H`$ is then found by maximizing the expression on the right-hand side of this latter inequality with respect to the point of contact $`r=t.`$ The equation which determines the critical point $`\widehat{t}`$ has to be solved by observing that, at the critical point $`r=\widehat{t},`$ the effective Coulomb coupling $`v(t)`$ still has to satisfy $`v(\widehat{t})<\beta v_P;`$ this nontrivial requirement imposes the announced “Coulomb coupling constant constraint” on the coupling parameters entering in the potential $`V(r).`$

It goes without saying that, when applying our “Herbst-based” lower energy bounds of Eq. (11), care has to be taken in order to assure that the chosen values of the parameters in the potential $`V(r)`$ under consideration respect this Coulomb coupling constant constraint.

Let us illustrate these general considerations by applying them to a Coulomb-plus-linear or (in view of its shape) “funnel” potential

$$V(r)=\frac{c_1}{r}+c_2r,c_10,c_20.$$
(14)
In this case, the effective Coulomb coupling $`v(t)`$ explicitly reads $`v(t)=c_1+c_2t^2`$ and, for the set of “Herbst-like” lower energy bounds $`e(v)`$ given by Eqs. (9) and (11), the critical point $`\widehat{t}`$ is, after some algebra, determined by the relation

$$v(\widehat{t})=\frac{\beta P^2}{\sqrt{m^2\widehat{t}^2+P^2}}.$$

The above Coulomb coupling constant constraint becomes an upper bound on a particular linear combination of the two involved coupling parameters $`c_1`$ and $`c_2`$ (cf. Sec. 6 of Ref. ):

$$c_1+\frac{P^2}{m^2}\left(\frac{P^2}{v_P^2}1\right)c_2<\beta v_P.$$
(15)
The Coulomb lower bounds (12) on the spectrum of the semirelativistic Hamiltonian (1) may be optimized by seeking that value of the parameter $`P=1/\sigma `$ which, for given values of the mass $`m`$ of the bound-state constituent(s) and of the coupling parameters $`c_1`$ and $`c_2`$ in the Coulomb-plus-linear potential (14), solves the Coulomb coupling constant constraint (15). Taking into account the definition (13) of the boundary value $`v_P`$ of the Coulomb coupling $`v,`$ this optimum value of $`P`$ is obtained as the solution of the relation

$$\frac{c_2\mathrm{sin}^4t}{\mathrm{cos}^2t(\beta \mathrm{sin}t\mathrm{cos}tc_1)}=m^2,P\mathrm{sin}t.$$
(16)
The set (12) of envelope lower bounds distinguished from each other by the parameter $`P`$ may be improved somewhat by taking into account the Martin–Roy lower bound (10) on the energy spectrum of the relativistic Coulomb problem. To this end, let us denote the energy bound on the right-hand side of the inequality (10) by $`e(v)`$:

$$e(v)=\beta m\sqrt{\frac{1}{2}\left[1+\sqrt{1\left(\frac{2v}{\beta }\right)^2}\right]}.$$

Performing the change of variables

$$v=\frac{\beta }{2}\mathrm{sin}2t,0t<\frac{\pi }{4},$$

we find

$$e(v)=\beta m\mathrm{cos}t,e(v)ve^{}(v)=\beta m\frac{\mathrm{cos}^3t}{\mathrm{cos}2t},$$

while from $`e^{}(v)=h(r)`$ we get

$$r=\frac{\mathrm{cos}2t}{m\mathrm{sin}t}.$$

Upon inserting these intermediate results in our principal envelope formula (2), we arrive at the “Martin–Roy-based” Coulomb lower bound

$$E_0\underset{0t<\pi /4}{\mathrm{min}}\left[\beta m\frac{\mathrm{cos}^3t}{\mathrm{cos}2t}+V\left(\frac{\mathrm{cos}2t}{m\mathrm{sin}t}\right)\right].$$

## 5 Concave transform of harmonic-oscillator potential

For the class of potentials $`V(r)`$ which are concave transforms $`V(r)=g(h(r)),`$ $`g^{\prime \prime }<0,`$ of the harmonic-oscillator potential $`h(r)=r^2`$ envelope upper bounds on the eigenvalues $`E_n\mathrm{}`$ of the spinless-Salpeter Hamiltonian (1) may be found for all energy levels (cf. Eq. (18) of Ref. ):

$$E_n\mathrm{}\underset{r>0}{\mathrm{min}}\left[\beta \sqrt{m^2+\frac{P^2}{r^2}}+V(r)\right],PP_n\mathrm{}(2)=2n+\mathrm{}\frac{1}{2}.$$
(17)
## 6 Concave transform of linear potential

For the class of potentials $`V(r)`$ which are concave transforms $`V(r)=g(h(r)),`$ $`g^{\prime \prime }<0,`$ of the linear potential $`h(r)=r`$ the “harmonic-oscillator based” envelope upper bounds (17) on the eigenvalues $`E_n\mathrm{}`$ of the spinless-Salpeter Hamiltonian (1) may be (significantly) improved by considering the composition of the envelope approximations for the kinetic-energy operator and the interaction-potential operator along the following lines. In the present investigation we focus our attention to attractive potentials $`V(r),`$ that is, to all potentials which exhibit a monotone increase with increasing radial coordinate $`r`$: $`g^{}>0.`$ Every function $`g(r)`$ that is a monotone increasing \[$`g^{}(r)>0`$\] and concave \[$`g^{\prime \prime }(r)<0`$\] function of $`r`$ is a concave function $`f(r^2)`$ of $`r^2`$: $`g(r)=f(r^2)`$ with $`f^{\prime \prime }(r^2)<0.`$ Consequently, all potentials $`V(r)`$ studied in this section also belong to the class of potentials considered in Sec. 5.

The square-root operator $`k(𝐩^2)\sqrt{m^2+𝐩^2}`$ of the relativistic kinetic energy, regarded as a function of $`𝐩^2,`$ is monotone increasing (i.e., $`k^{}(𝐩^2)>0`$) and concave (i.e., $`k^{\prime \prime }(𝐩^2)<0`$). Let us assume that the interaction potential $`V(r)`$ exhibits a similar behaviour, that is, that the operator $`V(r)=g(r)`$ is monotone increasing, i.e., $`g^{}(r)>0,`$ and concave, i.e., $`g^{\prime \prime }(r)<0.`$ These operators are bounded from above by their respective upper-tangent approximations:

$`k(𝐩^2)a_1(q)+b_1(q)𝐩^2,`$
$`g(r)a_2(t)+b_2(t)r,`$
with

$`b_1(q)=k^{}(q),a_1(q)=k(q)qk^{}(q),`$
$`b_2(t)=g^{}(t),a_2(t)=g(t)tg^{}(t),`$ (18)
where $`𝐩^2=q`$ denotes the point of contact between the kinetic-energy operator $`k(𝐩^2)`$ and its upper-tangent approximation $`a_1(q)+b_1(q)𝐩^2,`$ and $`r=t`$ labels the point of contact between the interaction-potential operator $`g(r)`$ and its upper-tangent approximation $`a_2(t)+b_2(t)r.`$ Consequently, our semirelativistic Hamiltonian $`H`$ in Eq. (1) satisfies the operator inequality

$$H\beta k(𝐩^2)+V(r)\overline{H}(q,t)\beta \left[a_1(q)+b_1(q)𝐩^2\right]+a_2(t)+b_2(t)r.$$

The eigenvalues $`\overline{E}_n\mathrm{}(q,t)`$ of the upper “tangential Hamiltonian” $`\overline{H}(q,t)`$ are given in terms of the numbers $`P_n\mathrm{}(1)`$ listed in Table 1 by the exact formula

$$\overline{E}_n\mathrm{}(q,t)=3\left[\frac{\beta b_1(q)b_2^2(t)P_n\mathrm{}^2(1)}{4}\right]^{1/3}+\beta a_1(q)+a_2(t).$$
(19)
In order to find the best such upper bound, we must minimize with respect to both $`q`$ and $`t.`$ The critical equations immediately yield the necessary conditions

$$\frac{1}{2\beta }\frac{g^{}(\widehat{t})}{k^{}(\widehat{q})}=\frac{\widehat{q}^{3/2}}{P_n\mathrm{}(1)}=\frac{P_n\mathrm{}^2(1)}{\widehat{t}^3}.$$
(20)
From the second of the above equalities we deduce that, for given quantum numbers $`n`$ and $`\mathrm{},`$ the minimizing values of the parameters $`q`$ and $`t`$ are constrained to the curve $`\widehat{q}\widehat{t}^2=P_n\mathrm{}^2(1)`$ in the $`(q,t)`$ parameter plane. On this curve, and with the coefficients (18) of our upper-tangent approximations and the necessary (critical) conditions (20), the minimum of the expression on the right-hand side of the upper energy bound (19) may be cast into a new, simpler form:

$$\underset{\widehat{t}>0}{\mathrm{min}}\left[\beta k\left(\frac{P_n\mathrm{}^2(1)}{\widehat{t}^2}\right)+g(\widehat{t})\right].$$

By remembering the above definitions of $`k(𝐩^2)`$ and $`g(r),`$ this expression eventually becomes

$$E_n\mathrm{}\underset{r>0}{\mathrm{min}}\left[\beta \sqrt{m^2+\frac{P^2}{r^2}}+V(r)\right],PP_n\mathrm{}(1).$$
(21)
These upper bounds on the entire discrete spectrum $`\{E_n\mathrm{}\}`$ are valid for any central potential $`V(r)`$ which is a monotone increasing and concave function of the linear potential $`h(r)=r.`$

## 7 Variational upper bounds

The standard tool for the derivation of rigorous upper bounds on the eigenvalues $`E`$ of some self-adjoint operator $`H`$ is the Rayleigh–Ritz variational technique , which is based on the minimum–maximum principle: If the (discrete) eigenvalues $`E_k,`$ $`k=0,1,2,\mathrm{},`$ of $`H`$ are ordered according to $`E_0E_1E_2\mathrm{},`$ then the first $`d`$ of them are bounded from above by the $`d`$ eigenvalues $`\widehat{E}_k,`$ $`k=0,1,\mathrm{},d1,`$ ordered according to $`\widehat{E}_0\widehat{E}_1\mathrm{}\widehat{E}_{d1},`$ of that operator $`\widehat{H}`$ which is obtained by the restriction of $`H`$ to a $`d`$-dimensional subspace $`D_d`$ of the domain of $`H,`$ that is, $`E_k\widehat{E}_k`$ for all $`k=0,1,\mathrm{},d1.`$ Hence, for the accuracy of the obtained results, an appropriate, “reasonable” definition of this trial subspace $`D_d`$ is crucial. If this $`d`$-dimensional subspace $`D_d`$ is spanned by a chosen set of $`d`$ linearly independent basis vectors $`|\psi _k,`$ $`k=0,1,\mathrm{},d1,`$ the set of eigenvalues $`\widehat{E}`$ may immediately be determined by diagonalizing the $`d\times d`$ matrix $`(\psi _i|\widehat{H}|\psi _j),`$ $`i,j=0,1,\mathrm{},d1,`$ that is, as the $`d`$ roots of the characteristic equation

$$det(\psi _i|\widehat{H}|\psi _j\widehat{E}\psi _i|\psi _j)=0,i,j=0,1,\mathrm{},d1.$$

The relativistic virial theorem derived in Ref. —for a very comprehensive review, see Ref. —allows to define a precise quantitative measure  for the quality of the results found within the framework of variational techniques: according to the analysis presented in Refs. , for some generic Hamiltonian operator $`H`$ consisting of a momentum-dependent kinetic-energy term $`T(𝐩)`$ and a coordinate-dependent interaction-potential term $`V(𝐱),`$ i.e.,

$$H=T(𝐩)+V(𝐱),$$

the accuracy of trial states $`|\phi `$ which approximate the exact bound state under study can be quantitatively estimated by the deviation from zero of the quantity (called $`\nu `$ in Refs. )

$$Q1\frac{\phi |𝐩\frac{}{𝐩}T(𝐩)|\phi }{\phi |𝐱\frac{}{𝐱}V(𝐱)|\phi }(\nu ).$$

The main advantage of this measure for the accuracy of approximate eigenstates $`|\phi `$ is that it does not require any information on the solutions of the investigated eigenvalue problem other than the one provided by the variational approximation technique itself.

In order to get a first idea of the position of the lowest energy level $`E_0,`$ a one-dimensional trial space is sufficient; a basis for the latter might be any of the two-parameter (normalized) trial functions

$$\varphi (𝐱)=\sqrt{\frac{n\alpha ^{3/n}}{4\pi \mathrm{\Gamma }(\frac{3}{n})}}\mathrm{exp}\left(\frac{\alpha }{2}r^n\right).$$

Here the scale parameter $`\alpha >0`$ is varied. The parameter $`n>0`$ allows for later optimization. For the sake of definiteness, let us consider an interaction potential $`V(r)`$ which is the sum of (attractive) pure power-law terms $`a(q)\mathrm{sgn}(q)r^q`$ with $`q0`$ and a logarithmic term $`a(0)\mathrm{ln}r`$:

$$V(r)=\underset{q0}{}a(q)\mathrm{sgn}(q)r^q+a(0)\mathrm{ln}r.$$

We assume, of course, that the potential coefficients $`a(q)0`$ are not all zero. With the help of Jensen’s inequality , and by a tricky redefinition of the scale variable $`\alpha ,`$ the variational upper bound on the ground-state eigenvalue $`E_0`$ of the Hamiltonian (1), resulting from this choice of $`\varphi (𝐱)`$ as basis vector of our one-dimensional trial space, may be cast into the form

$$E_0\underset{r>0}{\mathrm{min}}\left[\beta \sqrt{m^2+\frac{1}{r^2}}+\underset{q0}{}a(q)\mathrm{sgn}(q)(𝒫(n,q)r)^q+a(0)\mathrm{ln}(𝒫(n,0)r)\right],$$
(22)
with the abbreviations

$`𝒫(n,q)`$ $`=`$ $`{\displaystyle \frac{n}{2}}\left({\displaystyle \frac{\mathrm{\Gamma }(2+\frac{1}{n})}{\mathrm{\Gamma }(\frac{3}{n})}}\right)^{1/2}\left({\displaystyle \frac{\mathrm{\Gamma }(\frac{q+3}{n})}{\mathrm{\Gamma }(\frac{3}{n})}}\right)^{1/q},q0,`$
$`𝒫(n,0)`$ $`=`$ $`{\displaystyle \frac{n}{2}}\left({\displaystyle \frac{\mathrm{\Gamma }(2+\frac{1}{n})}{\mathrm{\Gamma }(\frac{3}{n})}}\right)^{1/2}\mathrm{exp}\left({\displaystyle \frac{1}{n}}\psi \left({\displaystyle \frac{3}{n}}\right)\right),`$
where $`\psi `$ denotes the digamma function

$$\psi (z)\frac{1}{\mathrm{\Gamma }(z)}\frac{\mathrm{d}\mathrm{\Gamma }(z)}{\mathrm{d}z}.$$

## 8 The “Laguerre” trial space

Unless the eigenstate corresponding to some eigenvalue $`E_k`$ is already an element of the trial space $`D_d,`$ any upper bounds on the eigenvalues of an operator bounded from below may be improved by enlarging $`D_d`$ to higher $`d,`$ or by spanning $`D_d`$ by more sophisticated basis states.

For the class of spherically symmetric potentials $`V(r)`$ studied here, a very popular choice for the basis states which span the $`d`$-dimensional trial space $`D_d`$ of the variational technique are “Laguerre” trial states, given in their configuration-space representation by

$$\psi _{k,\mathrm{}m}(𝐱)=\sqrt{\frac{(2\mu )^{2\mathrm{}+2\rho +1}k!}{\mathrm{\Gamma }(2\mathrm{}+2\rho +k+1)}}r^{\mathrm{}+\rho 1}\mathrm{exp}(\mu r)L_k^{(2\mathrm{}+2\rho )}(2\mu r)𝒴_\mathrm{}m(\mathrm{\Omega }_𝐱).$$
(23)
Here $`L_k^{(\gamma )}(x)`$ are the generalized Laguerre polynomials (for the parameter $`\gamma `$) , defined by the power series

$$L_k^{(\gamma )}(x)=\underset{t=0}{\overset{k}{}}(1)^t\left(\begin{array}{c}k+\gamma \\ kt\end{array}\right)\frac{x^t}{t!}$$

and orthonormalized, with the weight function $`x^\gamma \mathrm{exp}(x),`$ according to

$$\underset{0}{\overset{\mathrm{}}{}}dxx^\gamma \mathrm{exp}(x)L_k^{(\gamma )}(x)L_k^{}^{(\gamma )}(x)=\frac{\mathrm{\Gamma }(\gamma +k+1)}{k!}\delta _{kk^{}};$$

$`𝒴_\mathrm{}m(\mathrm{\Omega })`$ are the spherical harmonics for angular momentum $`\mathrm{}`$ and projection $`m,`$ depending on the solid angle $`\mathrm{\Omega }`$ and orthonormalized according to

$$d\mathrm{\Omega }𝒴_\mathrm{}m^{}(\mathrm{\Omega })𝒴_\mathrm{}^{}m^{}(\mathrm{\Omega })=\delta _{\mathrm{}\mathrm{}^{}}\delta _{mm^{}}.$$

The trial functions (23) involve two variational parameters, $`\mu `$ (with the dimension of mass) and $`\rho `$ (dimensionless), which, by the requirement of normalizability of these functions, are subject to the constraints $`\mu >0`$ and $`2\rho >1.`$

One of the main advantages of the choice (23) for the variational trial states is the easy availability of an analytic expression for their momentum-space representation, obtained by Fourier transformation of $`\psi _{k,\mathrm{}m}(𝐱)`$:

$`\stackrel{~}{\psi }_{k,\mathrm{}m}(𝐩)`$ $`=`$ $`\sqrt{{\displaystyle \frac{(2\mu )^{2\mathrm{}+2\rho +1}k!}{\mathrm{\Gamma }(2\mathrm{}+2\rho +k+1)}}}{\displaystyle \frac{(\mathrm{i})^{\mathrm{}}|𝐩|^{\mathrm{}}}{2^{\mathrm{}+1/2}\mathrm{\Gamma }\left(\mathrm{}+\frac{3}{2}\right)}}{\displaystyle \underset{t=0}{\overset{k}{}}}{\displaystyle \frac{(1)^t}{t!}}\left(\begin{array}{c}k+2\mathrm{}+2\rho \\ kt\end{array}\right)`$
$`\times `$ $`{\displaystyle \frac{\mathrm{\Gamma }(2\mathrm{}+\rho +t+2)(2\mu )^t}{(𝐩^2+\mu ^2)^{(2\mathrm{}+\rho +t+2)/2}}}F({\displaystyle \frac{2\mathrm{}+\rho +t+2}{2}},{\displaystyle \frac{\rho +t}{2}};\mathrm{}+{\displaystyle \frac{3}{2}};{\displaystyle \frac{𝐩^2}{𝐩^2+\mu ^2}})`$
$`\times `$ $`𝒴_\mathrm{}m(\mathrm{\Omega }_𝐩),`$
with the hypergeometric series $`F,`$ defined by

$$F(u,v;w;z)=\frac{\mathrm{\Gamma }(w)}{\mathrm{\Gamma }(u)\mathrm{\Gamma }(v)}\underset{n=0}{\overset{\mathrm{}}{}}\frac{\mathrm{\Gamma }(u+n)\mathrm{\Gamma }(v+n)}{\mathrm{\Gamma }(w+n)}\frac{z^n}{n!}.$$

The choice of the Laguerre states as the basis of the trial space $`D_d`$ allows to calculate the matrix elements of (any linear combination of) pure power-law potentials $`V(r)=a(q)r^q`$ and under certain circumstances (for instance, for $`\mu =m`$ for either radial or orbital excitations) the matrix elements of the kinetic term on purely analytical grounds, and to arrive therefore at analytical expressions for the $`d\times d`$ Hamiltonian matrix $`(\psi _i|\widehat{H}|\psi _j),`$ $`i,j=0,1,\mathrm{},d1.`$ The explicit algebraic expressions of these matrix elements may be found in Refs. . Up to and including a trial-space dimension $`d=4,`$ the Hamiltonian matrix $`(\psi _i|\widehat{H}|\psi _j)`$ may be diagonalized algebraically. For the funnel potential (14), the case $`d=1`$ yields the bound

$$E_0\left(\frac{64\beta }{15\pi }c_1\right)m+\frac{3c_2}{2m}.$$

## 9 The “local-energy” theorem

Some information on both upper and lower bounds on the isolated eigenvalues $`E_n\mathrm{}`$ of some Hamiltonian $`H`$ may be gained with the help of the so-called “local-energy” theorem . Unfortunately, the proof of this fundamental criterion makes use of the “nodal theorem” for the eigenstates of $`H.`$ As a consequence of this, for Hamiltonians $`H`$ defined in more than one dimension the local-energy theorem can be applied only to the ground state of the spectrum. It has been applied to the (three-dimensional) “relativistic Coulomb problem” in Ref. .

Upper bounds on the energy levels of Hamiltonians bounded from below can be obtained with considerably more efficiency by the variational technique described in Sec. 7. Therefore we employ here the local-energy theorem for the derivation of lower bounds on the spectrum of $`H.`$ In order to formulate the local-energy theorem in momentum-space representation, we introduce the “local energy”

$$(𝐩)\beta \sqrt{m^2+𝐩^2}+\frac{{\displaystyle \mathrm{d}^3q\stackrel{~}{V}(𝐩𝐪)\varphi (𝐪)}}{\varphi (𝐩)}.$$

Here, $`\stackrel{~}{V}(𝐩)`$ is the Fourier transform of the interaction potential $`V(𝐱)`$ under consideration:

$$\stackrel{~}{V}(𝐩)\frac{1}{(2\pi )^3}\mathrm{d}^3x\mathrm{exp}(\mathrm{i}𝐩𝐱)V(𝐱).$$

$`\varphi (𝐩)`$ denotes a suitably chosen, positive trial function: $`\varphi (𝐩)>0.`$ The “lower-bound part” of the local-energy theorem then states that the lowest-lying eigenvalue $`E_0`$ of the Hamiltonian $`H`$ is bounded from below by the infimum of $`(𝐩)`$ with respect to the momentum variable $`𝐩`$:

$$E_0\underset{𝐩}{inf}(𝐩).$$
(25)
The application of the local-energy theorem to confining potentials demands special care for the following reason: the Fourier transform of any confining potential is not well-defined; the potential has to be regarded as a distribution and an appropriate regularization must be applied. For instance, the linear potential $`V_\mathrm{L}(r)=\lambda r,`$ where $`\lambda `$ is a coupling parameter with mass dimension 2, may be regularized by the exponential $`\mathrm{exp}(\epsilon r)`$, $`\epsilon 0,`$ or by writing it as the derivative of this exponential or as the second derivative of a Yukawa-type function:

$$V_\mathrm{L}(r)=\lambda r=\lambda \underset{\epsilon 0}{lim}r\mathrm{exp}(\epsilon r)=\lambda \underset{\epsilon 0}{lim}\frac{}{\epsilon }\mathrm{exp}(\epsilon r)=\lambda \underset{\epsilon 0}{lim}\frac{^2}{\epsilon ^2}\frac{\mathrm{exp}(\epsilon r)}{r}.$$

With these regularizations, the Fourier transform $`\stackrel{~}{V}_\mathrm{L}(𝐩)`$ of the linear potential $`V_\mathrm{L}(r)`$ reads

$$\stackrel{~}{V}_\mathrm{L}(𝐩)=\stackrel{~}{V}_\mathrm{L}(|𝐩|)=\frac{\lambda }{\pi ^2}\underset{\epsilon 0}{lim}\frac{𝐩^23\epsilon ^2}{(𝐩^2+\epsilon ^2)^3}.$$

Note that (as a consistency check) the integral of the Fourier transform $`\stackrel{~}{V}_\mathrm{L}(𝐩)`$ has to vanish: $`\mathrm{d}^3p\stackrel{~}{V}_\mathrm{L}(𝐩)=0.`$ This holds for every power-law potential $`V(r)=a_nr^n`$ with exponent $`n>0.`$ The singularity structure of $`\stackrel{~}{V}_\mathrm{L}(𝐩)`$ becomes manifest when writing the derivative in $`V_\mathrm{L}(r)`$ as a difference quotient,

$$V_\mathrm{L}(r)=\lambda r=\lambda \underset{\epsilon 0}{lim}\frac{\mathrm{exp}(\epsilon r)1}{\epsilon },$$

which yields

$$\stackrel{~}{V}_\mathrm{L}(𝐩)=\stackrel{~}{V}_\mathrm{L}(|𝐩|)=\lambda \underset{\epsilon 0}{lim}\left[\frac{1}{\epsilon }\delta ^{(3)}(𝐩)\frac{1}{\pi ^2(𝐩^2+\epsilon ^2)^2}\right].$$

Thus the Fourier transform $`\stackrel{~}{V}_\mathrm{L}(𝐩)`$ of the linear potential consists of a (negative) regular part and a singular part, which is a distribution localized at the origin $`𝐩=\mathrm{𝟎}.`$ However, a detailed inspection shows that for the case of the funnel potential the momentum-space ground-state eigenfunction is indeed positive. This fact justifies our use of a positive trial function $`\varphi (𝐩).`$

The simplest normalizable trial function which comes to one’s mind is the exponential:

$$\varphi (r)=\sqrt{\frac{\mu ^3}{\pi }}\mathrm{exp}(\mu r).$$

This choice for the ground-state trial function $`\varphi (r)`$ corresponds to the case $`k=\mathrm{}=m=0`$ and $`\rho =1`$ of the Laguerre basis functions $`\psi _{k,\mathrm{}m}(𝐱)`$ introduced in Sec. 8. It is straightforward to calculate the Fourier transform of $`\varphi (r),`$ or to extract it from the general expression given explicitly in Sec. 8:

$$\varphi (𝐩)=\varphi (|𝐩|)=\frac{\sqrt{2\mu ^3}}{\pi |𝐩|(𝐩^2+\mu ^2)}\mathrm{sin}\left(2\mathrm{arctan}\frac{|𝐩|}{\mu }\right)=\frac{2\sqrt{2\mu ^5}}{\pi (𝐩^2+\mu ^2)^2}.$$

In order to illustrate the application of the local-energy theorem, let us consider the example of a potential $`V(r)`$ which is a sum of pure power-law terms:

$$V(r)=\underset{q}{}a(q)r^q.$$

For this power-law potential and the ground-state trial function $`\varphi (r),`$ it is straightforward to obtain the corresponding local energy $`(𝐩)`$ by Fourier transformation of the product $`r^q\varphi (r)`$:

$$(𝐩)=\beta \sqrt{m^2+𝐩^2}+\underset{q}{}\frac{a(q)\mathrm{\Gamma }(q+2)}{2\mu |𝐩|(𝐩^2+\mu ^2)^{q/21}}\mathrm{sin}\left[(q+2)\mathrm{arctan}\frac{|𝐩|}{\mu }\right].$$

The variational parameter $`\mu `$ entering into our ground-state trial function $`\varphi (r)`$ either may be given a chosen fixed value, or it may be adjusted in order to maximize the lower bound (25).

## 10 Sums of distinct potential terms

For the class of potentials $`V(r)`$ which consist of two or more distinct terms $`V^{(i)}(r)V^{(j)}(r)`$ for $`ij,`$ where, for every single component problem defined by the “one-term” Hamiltonian $`H_i\beta \sqrt{m^2+𝐩^2}+V^{(i)}(r),`$ information about the bottom of its spectrum is available, these individual pieces of information may be combined to form a lower bound on the spectrum of the Hamiltonian (1) with such a sum potential $`V(r).`$ The precise requirements for this “sum approximation” to work are the following:

1. The interaction potential $`V(r)`$ entering in the semirelativistic Hamiltonian of Eq. (1) is a sum of (more than one) different components $`V^{(i)}(r)=c_ih^{(i)}(r)`$:

$$V(r)=\underset{i}{}c_ih^{(i)}(r).$$
2. Every component problem $`H_i\beta \sqrt{m^2+𝐩^2}+c_ih^{(i)}(r)`$ supports, for sufficiently large values of its coupling parameter $`c_i,`$ a discrete eigenvalue at the bottom of its spectrum.

The sum approximation has already been studied for the case of Schrödinger Hamiltonians (involving nonrelativistic kinetic energies) in Refs. . The analysis of the corresponding semirelativistic problem has been presented in full generality in Ref. . Therefore, here we summarize only the result for the case in which the potential $`V(r)`$ is a sum of attractive pure power-law terms:

$$V(r)=\underset{q0}{}a(q)\mathrm{sgn}(q)r^q,$$
(26)
where the coupling parameters $`a(q)`$ are non-negative and not all zero. For this potential, a lower bound to the bottom of the spectrum of the semirelativistic Hamiltonian (1), which is, by assumption, the lowest eigenvalue $`E_0`$ of $`H,`$ is given by the expression

$$E_0\underset{r>0}{\mathrm{min}}\left[\beta \sqrt{m^2+\frac{1}{r^2}}+\underset{q0}{}a(q)\mathrm{sgn}(q)(P(q)r)^q\right],$$
(27)
if, for the particular potential $`V(r)`$ under consideration, there exist (or there can be found) $`P`$ numbers $`P(q)`$ such that, whenever there is only one term present in the potential (26), the right-hand side of the inequality (27) yields either the exact ground-state energy eigenvalue of the corresponding “single-term” problem or a lower bound to it.

Consequently, when applying our general sum lower bound (27) to a given sum potential $`V(r),`$ the main task is to derive appropriate $`P`$ numbers for all components in the sum (26). Let us illustrate this by considering the two examples relevant for the funnel potential (14).

* For a Coulomb component potential $`h^{(i)}(r)=1/r,`$ corresponding to the case $`q=1`$ in Eq. (26), we may take advantage from the explicit results presented in Sec. 4. By the change of variables $`rPr,`$ the Coulomb lower bound (12) may be cast into the form

$$E_0\underset{r>0}{\mathrm{min}}\left[\beta \sqrt{m^2+\frac{1}{r^2}}+V(Pr)\right],v<\beta v_P.$$

Specifying this general result to the Coulomb potential $`V(r)=v/r`$ yields the bound

$$E_0\underset{r>0}{\mathrm{min}}\left[\beta \sqrt{m^2+\frac{1}{r^2}}\frac{v}{Pr}\right],v<\beta v_P.$$

This bound is, of course, nothing else but, for $`P=2/\pi ,`$ the Herbst lower bound (9) or, for arbitrary $`P,`$ the Herbst-like lower bound (11) to the relativistic Coulomb problem. Consequently, for the Coulomb component term we may identify the $`P`$ number $`P(1)`$ required for the lower bound (27) with the lower-bound parameter $`P`$ defined in Sec. 4: $`P(1)=P.`$
* For a linear component potential $`h^{(i)}(r)=r,`$ that is, for the case $`q=1`$ in Eq. (26), (a possible value of) the $`P`$ number $`P(1)`$ required for the lower bound (27) turns out to be related to the ground-state eigenvalue $`e_0(v)`$ of the Hamiltonian operator $`\sqrt{𝐩^2}+vr`$ by

$$e_0(v)=2\sqrt{vP(1)};$$

a numerical determination of this lowest energy eigenvalue $`e_0(v)`$ yields $`P(1)=1.2457.`$

## 11 Comparison: The case of the funnel potential

In order to demonstrate the power of the different energy bounds derived so far, let us apply the above results to the very illustrative example of the Coulomb-plus-linear potential (14). By factorizing off an overall coupling constant $`0<v1,`$ we write this potential in the form

$$V(r)=v\left(\frac{a}{r}+br\right).$$
(28)
For the coupling parameters entering in this funnel potential $`V(r),`$ different notations have been used in Eqs. (14), (26), and (28). Evidently, the three sets of potential coefficients have to be identified according to $`a(1)c_1av>0`$ and $`a(1)c_2bv>0.`$

Fig. 1 compares the various energy bounds derived or investigated in Secs. 4 through 10 for the ground-state energy eigenvalue $`E_0`$ of the Hamiltonian (1) with the funnel potential (28) as a function of the overall coupling parameter $`v.`$ These energy bounds necessarily form a bunch of concave curves $`E(v),`$ all starting at the (free-energy) value $`E_0(v=0)=\beta m=1.`$

The Coulomb lower bound (12) may be applied to the funnel potential either with a fixed value of the lower-bound parameter $`P`$ or with the optimized $`P(v)`$ computed from Eq. (16). For “fixed-$`P`$” bounds, the maximum value of $`P`$ allowed by the Coulomb coupling constant constraint discussed in Sec. 4 is found, for given values of the kinetic-term parameters $`\beta ,`$ $`m`$ and potential coefficients $`c_1,`$ $`c_2,`$ as solution of that relation which is obtained when equating the left-hand and right-hand sides of the inequality (15) with the critical coupling $`v_P`$ given by the definition (13), that is, by solving, for the maximum values of $`c_1`$ and $`c_2`$ occuring here, the constraint (16):

$$c_1+\frac{P^4}{1P^2}\frac{c_2}{m^2}=\beta P\sqrt{1P^2}.$$

For a single, unit-mass particle, i.e., $`\beta =1`$ and $`m=1,`$ and the values of the funnel-potential coupling parameters $`a`$ and $`b`$ used in Fig. 1, $`a=0.2`$ and $`b=0.5,`$ this yields $`P=0.728112.`$

We try to approach the state we are interested in as closely as possible by the application of the Rayleigh–Ritz variational technique briefly sketched in Sec. 7, with the trial space $`D_d`$ spanned by the Laguerre basis defined in Sec. 8. By fixing the variational parameter $`\mu `$ to the value $`\mu =m,`$ we are able to obtain the matrix elements of the Hamiltonian (1) analytically, avoiding thereby the necessity to use a numerical integration procedure for the evaluation of these matrix elements; for definiteness, we fix the variational parameter $`\rho `$ to the value $`\rho =1.`$

The accuracy of such an approximation to an eigenstate may be estimated with the help of the “virial-theorem-inspired” quality measure $`Q`$ discussed in Sec. 7. Fig. 2 shows, for the trial-space dimension $`d=25`$ chosen for the variational bound (E) in Fig. 1, the dependence of the accuracy $`Q`$ on the overall coupling $`v.`$ We obtain $`Q=|Q|1.1\times 10^4`$ over the whole parameter range $`0<v1.`$ Such high accuracy should be sufficient for the present purpose.

For the sake of completeness, we plot, in Fig. 3, the behaviour of the ground state of the Hamiltonian (1) with the Coulomb-plus-linear potential (28), for the maximum value $`v=1`$ of the overall coupling parameter $`v,`$ determined by the variational technique of Sec. 7, with the Laguerre basis of Sec. 8, in both configuration and momentum space. As expected, the ground-state momentum-space eigenfunction possesses no node, which is in agreement with our assumption of a positive momentum-space trial function $`\varphi (𝐩)`$ (which is also depicted, as the case $`d=1,`$ in Fig. 3) in the derivation of the local-energy lower energy bound in Sec. 9.

## 12 Summary and conclusions

This review of the discrete spectra of semirelativistic Hamiltonians $`H`$ of the form (1) should be considered as a compilation of rigorous (semi-)analytical results for both upper and lower bounds on the discrete eigenvalues $`E_n\mathrm{}`$ of $`H,`$ presented in a way suitable for immediate use: the algebraic structure and the convexity properties with respect to suitable basis potentials $`h(r)`$ of the interaction potential $`V(r)`$ one wishes to use in one’s semirelativistic Hamiltonian (1) unambiguously determine which of all the energy bounds discussed here actually apply. With the explicit analytical expressions for these bounds on the discrete eigenvalues at one’s disposal, one is left with the (numerical) optimization with respect to a single real variable. Clearly, within the respective limitations, the Rayleigh–Ritz variational technique of Secs. 7 and 8 as well as the local-energy theorem of Sec. 9 apply to arbitrary reasonable potentials.

For the particular example considered in Sec. 11, viz., the Coulomb-plus-linear potential (28), the sum lower bound turns out to be somewhat better than the lower bound provided by the local-energy theorem. Clearly, the sum lower bound cannot always be better than the local-energy bound because one could, in principle, choose in the local-energy theorem (25) the exact eigenfunction for the trial function $`\varphi `$ and thus recover the exact energy eigenvalue.

In conclusion, let us stress that all but two of the various energy bounds we have studied in this work may eventually be expressed by essentially the same general semi-classical form. We can view this in the following way: Since the square of the momentum $`𝐩^2`$ scales like $`1/r^2,`$ we expect that

$$𝐩^2=\frac{P^2}{r^2},$$

where $`P`$ is a suitable coefficient; thus the balance between the kinetic and potential energies required by the minimum–maximum principle might lead to a search for the minimum of an energy expression of the form

$$\sqrt{m^2+\frac{P^2}{r^2}}+V(r).$$

What makes this heuristic argument interesting is the fact that definite values of $`P`$ can be prescribed so that this minimum provides a priori a bound on the (unknown) exact energy. In the case where $`V(r)`$ is a sum of terms, a trivial change of variables is used to move the $`P`$ parameters inside the potential $`V(r),`$ and then optimal distinct parameters can be used for each potential term.

The lower bound by the sum approximation incorporates our best known lower bounds for the sub-problems containing only a single potential term; the final result is the optimized mixture of these contributions. The lower bound by the local-energy theorem (25) could, in principle, achieve arbitrarily high accuracy. The difficulty here, as with variational methods, is to find a suitable trial function. To the extent that we have searched we did not manage to find a better lower bound with the help of the local-energy theorem (25) than that provided immediately by the sum approximation (27).

Our most accurate results overall were upper bounds obtained by using a 25-dimensional trial space spanned by Laguerre functions. We estimate the error for these results to be less than $`0.01\%`$ over the total range of potential parameters studied; such accurate estimates are invaluable for estimating the effectiveness of the large variety of alternative simpler energy bounds which we have reviewed in this article.


format: text

text:
# Stochastic synchronization in globally coupled phase oscillators

## 1 Introduction

Recently, various noise effects for nonlinear systems have been studied. In stochastic resonance, the response of a bistable system or an excitable system to a periodic force is enhanced with the addition of noise. The stochastic resonance improves signal detection by the superposed noise . Noise-enhanced frequency locking is observed in stochastic bistable systems driven by a relatively strong periodic force or a chaotic signal and phase diagrams similar to the Arnold tongues are obtained . Noise induced entrainment among coupled oscillators is found experimentally in Belousov-Zhabotinsky reactions and brain waves . Stochastic resonance has been studied also in globally or locally coupled systems of many bistable elements and it was shown that the coupling can lead to the enhancement of the response . Various types of Fokker-Planck equations have been used to study the stochastic resonance theoretically.

On the other hand, various types of collective dynamics in globally coupled oscillators have been studied . Globally coupled phase oscillators under external noises can be studied with a nonlinear diffusion equation for the number density. The nonlinear diffusion equation is obtained from the mean field approximation of the Fokker-Planck equation for the $`N`$ oscillators, which is considered to be correct in the globally coupled system in the limit $`N\mathrm{}`$ . The nonlinear diffusion equation for the globally coupled phase oscillators can be transformed into coupled nonlinear equations for the Fourier amplitudes of the number density. The numerical simulation of the coupled nonlinear equations for the Fourier modes is relatively easy. Various nonequilibrium phase transitions were found in the globally coupled phase oscillators using the numerical simulation of the coupled equations for the Fourier modes and their bifurcation analysis .

In this paper, we study an extended model of the coupled phase oscillators, in which an external periodic force is added to the model studied in . Each element may exhibit stochastic resonance and stochastic synchronization, however, we study dynamical behaviors of the order parameter. The nonlinear diffusion equation is a deterministic equation, and the stochastic behaviors for each element are averaged out in the description of the order parameter. The resonant response and the synchronization to a periodic force are more clearly shown in this deterministic system.

## 2 Stochastic resonance in globally coupled bistable systems

At first, we consider a globally coupled bistable system. Each phase oscillator evolves according to an equation

$$\frac{d\varphi }{dt}=b\mathrm{sin}2\varphi c\mathrm{sin}\omega _0t\mathrm{sin}\varphi ,$$
(1)
where $`\varphi `$ is the phase of the oscillator, $`b`$ and $`c`$ are positive constants, and $`\omega _0`$ is the frequency of the periodic force. There are two stationary solutions: $`\varphi =0`$ and $`\pi `$. If the coefficient of the second term on the right-hand side is not temporally periodic as $`c\mathrm{sin}\omega _0t`$ but constant $`c_0`$, the two solutions are bistable for $`|c_0|2b`$. In that case, the dynamical system has a potential function $`U(\varphi )=b/2\mathrm{cos}2\varphi c_0\mathrm{cos}\varphi `$, and the solution $`\varphi =0`$ has a lower potential for $`c_0>0`$ and a higher potential for $`c_0<0`$ than the solution $`\varphi =\pi `$. A model of a globally coupled system with a noise term is written as

$$\frac{d\varphi _i}{dt}=b\mathrm{sin}2\varphi _ic\mathrm{sin}\omega _0t\mathrm{sin}\varphi _i\frac{K}{N}\underset{j=1}{\overset{N}{}}\mathrm{sin}(\varphi _i\varphi _j)+\xi _i,$$
(2)
where $`N`$ is the total number of elements, and the second term on the right-hand side represents the periodically forcing term and the third term represents mutual coupling, and the last term represents Gaussian white noise characterized by

$$\xi _i(t)=0,\xi _i(t)\xi _j(t^{})=2D\delta _{ij}\delta (tt^{}),$$

where $`D`$ represents noise strength. The mean field treatment holds exactly in the limit of $`N\mathrm{}`$. The normalized number density defined by

$$n(\varphi ,t)=\frac{1}{N}\underset{m}{}\underset{j}{}\delta (\varphi _j\varphi 2\pi m)$$

obeys a nonlinear diffusion equation

$`{\displaystyle \frac{}{t}}n(\varphi ,t)`$ $`=`$ $`{\displaystyle \frac{}{\varphi }}\left\{b\mathrm{sin}2\varphi c\mathrm{sin}\omega _0t\mathrm{sin}\varphi K{\displaystyle _\pi ^\pi }𝑑\varphi ^{}\mathrm{sin}(\varphi \varphi ^{})n(\varphi ^{},t)\right\}n(\varphi ,t)`$ (3)
$`+D{\displaystyle \frac{^2}{\varphi ^2}}n(\varphi ,t)`$
The $`2\pi `$-periodic function $`n(\varphi ,t)`$ can be expanded as

$$n(\varphi ,t)=\frac{1}{2\pi }\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\rho _m(t)e^{im\varphi }$$

The nonlinear equation (3) is rewritten with coupled ordinary differential equations,

$$\frac{d\rho _m}{dt}=m\{K/2(\rho _1\rho _{m1}\rho _1\rho _{m+1})+b/2(\rho _{m2}\rho _{m+2})+c/2\mathrm{sin}\omega _0t(\rho _{m1}\rho _{m+1})\}Dm^2\rho _m.$$
(4)
We have performed numerical simulation of Eq. (4) retaining the first 50 modes with the Runge-Kutta method of timestep $`\mathrm{\Delta }t=0.0001`$. The order parameter is expressed as

$$\sigma =(1/N)\underset{j=1}{\overset{N}{}}e^{i\varphi _j}=_\pi ^\pi 𝑑\varphi e^{i\varphi }n(\varphi ,t)=\rho _1(t)=\rho _1^{}.$$

In this model, the order parameter takes a real value, that is, $`\mathrm{Im}\rho _1(t)=0`$. Figure 1 displays three time evolutions of $`\sigma (t)`$ at $`D=0.4,0.57`$ and 1 for $`b=0.2,c=0.05,\omega _0=0.3`$ and $`K=1`$. The order parameter is regularly oscillating with the frequency $`\omega _0`$ of the external force. The average value of the order parameter takes a nonzero value at $`D=0.4`$ and zero at $`D=0.57`$ and 1. The amplitude of the periodic oscillation is maximum at $`D=0.57`$. We have calculated the average values $`\sigma `$ of the order parameter and the temporal fluctuations around the average value, $`(\delta \sigma (t))^2^{1/2}`$, where $`\delta \sigma (t)=\sigma (t)\sigma `$. The results are shown in Fig. 2. The average values of order parameters take nonzero values for $`D<D_c0.56`$. This is a kind of symmetry breaking phase transition induced by noises. The amplitude of the periodic oscillation of the order parameter takes a maximum near the phase transition point. The response function to the periodic modulation becomes large in an intermediate range of noise strength. This is interpreted as a kind of stochastic resonance in this globally coupled bistable system. A phase transition in a globally coupled bistable system without a periodic force and the enhancement of the susceptibility to a weak periodic force was discussed in ,.

We have found that a kind of phase transition occurs even for nonzero amplitude of periodic modulation.

## 3 Stochastic synchronization in globally coupled phase oscillators

Next, we consider a globally coupled oscillator system. Each phase oscillator is assumed to obey an equation

$$\frac{d\varphi }{dt}=\omega (b+c\mathrm{sin}\omega _0t)\mathrm{sin}\varphi ,$$
(5)
where $`\varphi `$ is the phase of the oscillator, $`\omega ,b`$ and $`c`$ are positive constants, and $`\omega _0`$ is the frequency of the periodic force. If $`c=b=0`$, this equation describes a simple phase rotator, since the solution is $`\varphi =\omega t`$, where $`\omega `$ is the natural frequency. If $`c=0`$, and $`b0`$ but $`b<\omega `$, $`d\varphi /dt`$ is not constant but always takes a positive value, and this equation describes still a phase rotator. The natural frequency of the oscillator is given by $`\sqrt{\omega ^2b^2}`$. When $`c=0`$ and $`b<\omega `$ but $`b`$ is close to $`\omega `$, the oscillator behaves like the relaxation oscillation (it may be suitable to call it relaxation oscillation of phase), since the phase motion is clearly separated into fast motion and very slow motion near $`\varphi =\pi /2+2n\pi `$ ($`n`$ is an integer). On the other hand, when $`\omega <b`$ and $`c=0`$, the equation has a stable stationary solution given by $`\varphi =\mathrm{sin}^1(\omega /b)`$ and all trajectories are attracted to the stable point. If $`c=0`$, and $`b>\omega `$ but $`b`$ is close to $`\omega `$, the system behaves like an excitable system, since one-cycle rotation is excited by some small perturbation to the stationary state. The parameter $`b`$ is a parameter which determines a transition from an excitable system to an oscillatory system. The threshold value is $`b=b_c=\omega `$. The parametric perturbation $`b+c\mathrm{sin}\omega _0t`$ implies that the threshold value is periodically modulated. It is known in an experiment of light-sensitive Belousov-Zhabotinsky reaction that illumination can control the threshold . The noise entrainment was observed by controlling the illumination . The periodic modulation in our model may be related to such periodic modulation of the threshold value in experimental systems.

A model of a globally coupled system with a noise term is written as

$$\frac{d\varphi _i}{dt}=\omega (b+c\mathrm{sin}\omega _0t)\mathrm{sin}\varphi _i\frac{K}{N}\underset{j=1}{\overset{N}{}}\mathrm{sin}(\varphi _i\varphi _j)+\xi _i,$$
(6)
where $`N`$ is the total number of elements, and the third term represents mutual coupling, and the last term represents Gaussian white noise characterized by

$$\xi _i(t)=0,\xi _i(t)\xi _j(t^{})=2D\delta _{ij}\delta (tt^{}).$$

A coupled phase oscillator model similar to this equation was studied by Kim et al. . In their model, time delay was further assumed. In their system with time delay, there are multistable states, i.e., fast and slow synchronized states and a desynchronized state. They studied noise induced transitions among the multistable states in one oscillator system and a coupled system of 10 oscillators. We do not consider such time delay effect and study the stochastic synchronization using a nonlinear diffusion equation for a coupled system in the limit $`N\mathrm{}`$. The normalized number density obeys

$`{\displaystyle \frac{}{t}}n(\varphi ,t)`$ $`=`$ $`{\displaystyle \frac{}{\varphi }}\left\{\omega (b+c\mathrm{sin}\omega _0t)\mathrm{sin}\varphi K{\displaystyle _\pi ^\pi }𝑑\varphi ^{}\mathrm{sin}(\varphi \varphi ^{})n(\varphi ^{},t)\right\}n(\varphi ,t)`$ (7)
$`+D{\displaystyle \frac{^2}{\varphi ^2}}n(\varphi ,t)`$
The corresponding coupled ordinary differential equations are written as

$$\frac{d\rho _m}{dt}=m\{K/2(\rho _1\rho _{m1}\rho _1\rho _{m+1})+(b+c\mathrm{sin}\omega _0t)/2(\rho _{m1}\rho _{m+1})\}(im\omega +Dm^2)\rho _m.$$
(8)
We have performed numerical simulation of Eq. (8) retaining the first 50 modes with the Runge-Kutta method of timestep $`\mathrm{\Delta }t=0.0005`$. In this model, the order parameter $`\sigma =\rho _1`$ is a complex variable.

Figure 3 displays three time evolutions of $`\mathrm{Re}\sigma (t)`$ at $`D=0.04,\mathrm{\hspace{0.17em}0.08}`$ and 0.12 for $`\omega =1.005,K=1,\omega _0=0.3`$ and $`c=0.02`$. The order parameter is entrained to the external periodic force at $`D=0.08`$. However, the motion of the order parameter is not completely entrained to the periodic force at $`D=0.04`$ and 0.12, and quasi-periodic motions appear.

Figure 4(a) displays the frequency $`\overline{\omega }`$ of the motion of the order parameter as a function of $`D`$ for $`\omega =1.005,K=1,b=1\omega _0=0.3`$ and $`c=0.02`$. The frequency was calculated from the rotation number of $`(\mathrm{Re}\sigma (t),\mathrm{Im}\sigma (t))`$ around the origin $`(0,0)`$. The frequency of the order parameter increases with the noise strength $`D`$. It can be understood from the dynamical behavior of each oscillator. Each oscillator behaves like the relaxation oscillation and the phase motion becomes slow near $`\varphi =\pi /2+2n\pi `$ when $`D=0`$, since $`\omega b`$. An effect of noises is to make the regular slow motion randomly faster, and then the average frequency increases with the noise strength. The periodic modulation is further applied. The frequency entrainment is observed between $`0.058<D<0.094`$. This is a kind of the stochastic synchronization, although the motion of the order parameter is not stochastic. Figure 4(b) displays the Fourier amplitude of the order parameter motion with frequency $`\omega _0`$, that is, $`(1/T)|_0^T𝑑t\sigma (t)\mathrm{exp}(i\omega _0t)|`$. The Fourier amplitude takes a maximum in the parameter range where the stochastic synchronization is observed. Figure 4(c) displays a parameter region in the $`KD`$ parameter space, in which the synchronization is observed. This phase diagram seems to be the Arnold tongues discussed in \[4-6\], however, the ordinate is not the amplitude $`c`$ of the periodic forcing but the coupling constant $`K`$, so the meaning is different. For a fixed value of noise strength, the coupling constant also needs to take an intermediate value for the stochastic synchronization to appear, as shown in this figure. Too large coupling constant makes the distribution of the number density too narrow, and noise effects are effectively reduced.

Even for a parameter range $`\omega <bc`$, in which each elemental system is purely excitable without noises, noises can induce phase slips and the oscillation is observed on an average. When the periodic force is further applied, the stochastic synchronization can occur. We have performed a numerical simulation for $`\omega =0.97,b=1,K=2,\omega _0=0.3`$ and $`c=0.02`$. Figure 5(a) displays the frequency of the motion of the order parameter as a function of $`D`$. The frequency entrainment is observed for $`0.245<D<0.305`$. Figure 5(b) displays the Fourier amplitude of the order parameter motion with frequency $`\omega _0=0.3`$. The Fourier amplitude also has a maximum in an intermediate value of noise strength.

## 4 Summary

We have proposed globally coupled phase oscillator models and performed numerical simulations of the corresponding nonlinear diffusion equations for the number density. The response to the periodic force becomes large in an intermediate range of noise strength. It can be interpreted as a stochastic resonance or a stochastic synchronization, although the motion of the order parameter is regular. We have shown numerical results for parametrically perturbed system, but we have also studied an additively perturbed model:

$$\frac{d\varphi _i}{dt}=\omega b\mathrm{sin}\varphi _i+c\mathrm{sin}\omega _0t\frac{K}{N}\underset{j=1}{\overset{N}{}}\mathrm{sin}(\varphi _i\varphi _j)+\xi _i,$$

and another type of forced model:

$$\frac{d\varphi _i}{dt}=\omega b\mathrm{sin}\varphi _ic\mathrm{sin}(\varphi _i\omega _0t)\frac{K}{N}\underset{j=1}{\overset{N}{}}\mathrm{sin}(\varphi _i\varphi _j)+\xi _i.$$

We have observed similar stochastic synchronization also in these systems. It does not seem to depend on the detailed form of periodic forcing. The stochastic synchronization may be found in many systems driven by noises and periodic forcing, if the parameter for each element is near a transition point between an excitable state and an oscillatory state.


text:
# 1 Introduction

## 1 Introduction

To minimize its fine-tuning problems, such as initial conditions, there have been attempts to unify quintessence with the inflaton field in a single scalar . This way, introducing yet again another unobserved scalar field is avoided. Also, a common theoretical framework can be used to describe both inflation and quintessence. For quintessential inflation one needs a sterile inflaton, with only minimal gravitational coupling to the standard model, because the field should survive until today. The Universe reheats through gravitational particle production . The minimum of the potential (assumed zero) is typically placed at infinity, because it should not have yet been reached. This results in the so-called quintessential tail of the potential. Candidates for the quintessential–inflation scalar field are moduli fields or the radion field.

## 2 Dynamics

The Universe is modeled as a collection of perfect fluids; the background fluid with density $`\rho _B`$ (comprised by matter and radiation) and the scalar field $`\varphi `$ with density $`\rho _\varphi =\rho _{\mathrm{kin}}+V`$ and pressure $`p_\varphi =\rho _{\mathrm{kin}}V`$, where $`\rho _{\mathrm{kin}}\frac{1}{2}\dot{\varphi }^2`$ is the kinetic density of $`\varphi `$ and the dot denotes derivative w.r.t. the cosmic time $`t`$. For every component one defines the barotropic parameter as $`w_ip_i/\rho _i`$. The Universe expansion accelerates when $`\rho _B<\rho _\varphi `$ and $`w_\varphi <\frac{1}{3}`$. Energy-momentum conservation demands $`d(a^3\rho )=pd(a^3)`$, which, for decoupled fluids, gives $`\rho _ia^{3(1+w_i)}`$, where $`a`$ is the scale factor. To study the dynamics of the Universe, one also needs the Friedman equation and the $`\varphi `$ field equation: $`\dot{\rho }_{\mathrm{kin}}+6H\rho _{\mathrm{kin}}+\dot{V}=0`$, where $`H\dot{a}/a`$ is the Hubble parameter. In spatially-flat, FRW cosmology the Friedman equation is $`H^2=\rho /3m_P^2`$, with $`m_P`$ being the reduced Planck mass. This results in the evolution equations:

$`H={\displaystyle \frac{2t^1}{3(1+w)}}`$ $`at^{\frac{2}{3(1+w)}}`$ $`\rho ={\displaystyle \frac{4}{3(1+w)^2}}\left({\displaystyle \frac{m_P}{t}}\right)^2`$ (1)
where $`w`$ corresponds to the dominant fluid component. In brane cosmology the Friedman equation is modified for density above the string tension $`\lambda `$ . Assuming a 5th dimension, the Friedman equation becomes $`H\rho /\sqrt{6\lambda }m_P`$, where $`\lambda =\frac{3}{4\pi }(M_5^6/m_P^2)`$, with $`M_5`$ being the fundamental, 5-dim Planck mass. Then the evolution equations are:

$`H={\displaystyle \frac{t^1}{3(1+w)}}`$ $`at^{\frac{1}{3(1+w)}}`$ $`\rho ={\displaystyle \frac{\sqrt{6\lambda }}{3(1+w)}}\left({\displaystyle \frac{m_P}{t}}\right)`$ (2)
The modified dynamics affect the inflationary era due to excessive friction on the roll of the scalar field, which allows inflation with steep potential. Then, the slow-roll parameters become: $`ϵ2\lambda m_P^2(V^{})^2/V^3`$ and $`\eta 2\lambda m_P^2V^{\prime \prime }/V^2`$, where the prime denotes derivative w.r.t. $`\varphi `$. The modified slow-roll changes the amplitude and spectral index of the generated density perturbations :

$`{\displaystyle \frac{\delta \rho }{\rho }}{\displaystyle \frac{1}{2\sqrt{6}\pi }}{\displaystyle \frac{V^3}{\lambda ^{3/2}|V^{}|m_P^3}}`$ and $`n_s14m_P^2{\displaystyle \frac{\lambda }{V}}\left[\mathrm{\hspace{0.33em}3}\left({\displaystyle \frac{V^{}}{V}}\right)^2{\displaystyle \frac{V^{\prime \prime }}{V}}\right]`$ (3)
Gravitational reheating creates a thermal bath of temperature $`T_{\mathrm{reh}}=\frac{\alpha }{2\pi }H_{\mathrm{end}}`$ , where $`\alpha `$ is the reheating efficiency and the subscript ’end’ denotes the end of inflation. Due to the inefficiency of reheating, after the end of inflation the Universe becomes dominated by $`\rho _{\mathrm{kin}}`$, i.e. $`\rho \rho _{\mathrm{kin}}a^6`$. As long as $`\rho _{\mathrm{kin}}>2\lambda `$ we have $`at^{1/6}`$ and $`\varphi (t)=\varphi _{\mathrm{end}}+\frac{2}{\sqrt{3}}\sqrt{\lambda /V_{\mathrm{end}}}\left(\sqrt{t/t_{\mathrm{end}}}1\right)m_P`$. However, when $`\rho _{\mathrm{kin}}<2\lambda `$ the brane regime ends and one has $`at^{1/3}`$ and $`\varphi (t)=\varphi _{\mathrm{end}}+\frac{1}{\sqrt{6}}\left[2+\mathrm{ln}\left(t/t_\lambda \right)\right]m_P`$, where $`t_\lambda =\frac{1}{2\sqrt{6}}m_P/\sqrt{\lambda }`$ is the crossover time for which $`\varphi _\lambda =\varphi _{\mathrm{end}}+\sqrt{2/3}m_P`$. The thermal bath eventually dominates the density of the Universe and the Hot Big Bang begins at the temperature:

$$T_{}=\frac{\alpha ^3}{96\pi ^2}\sqrt{\frac{g_{}}{15}}\frac{V_{\mathrm{end}}^{5/2}}{\lambda ^{3/2}m_P^3}$$
(4)
where $`g_{}10^2`$ is the number of relativistic degrees of freedom. In order not to affect Nucleosynthesis (BBN) we require $`T_{}>1`$ MeV. After the onset of the Hot Big Bang $`\rho _{\mathrm{kin}}`$ rapidly decreases and the field freezes at the value:

$$\varphi _F=\varphi _{\mathrm{end}}+\frac{m_P}{\sqrt{6}}\left[14+\frac{3}{2}\mathrm{ln}\left(\frac{30\pi ^2}{g_{}\alpha ^4}\right)4\mathrm{ln}\left(\frac{m_P^4}{\lambda }\right)+5\mathrm{ln}\left(\frac{m_P^4}{V_{\mathrm{end}}}\right)\right]$$
(5)
## 3 Why Branes

Standard–cosmology quintessential inflation needs a potential with two flat regions: the inflationary plateau and the quintessential tail . BBN and coincidence demand that the density scale of the flat regions differs by $`10^{100}`$! To prepare for this abysmal dive, $`V`$ is strongly curved near the end of inflation. As a result, the slow-roll parameter $`\eta `$ is large, leading to spectral index $`n_s`$ far from unity. This makes it difficult to construct single-branch models111Multi-branch models such as  are disfavoured in our minimalistic approach. for quintessential inflation, although solutions do exist . In contrast, this $`\eta `$-problem does not appear when considering brane cosmology because the modified Friedman equation allows for steep inflation . Hence, one may obtain successful quintessential inflation with a simple potential.

## 4 The exponential tail

String theory disfavours eternal acceleration . This immediately rules out quintessential tails milder than exponential and also frozen quintessence (in which $`\varphi =\varphi _F`$ at present). In any case, coincidence is hard to achieve with mild tails (e.g. inverse power-law type ) and frozen quintessence is virtually indistinguishable from the cosmological constant alternative. On the other hand steeper-than-exponential tails have disastrous attractors . Thus, the most reasonable approach is the exponential tail, $`VV_0e^{b\varphi /m_P}`$, where $`b`$ is a positive constant. Then, the $`\varphi `$-field equation gives an exact attractor solution:

$`\varphi _{\mathrm{attr}}={\displaystyle \frac{2m_P}{b}}\mathrm{ln}\left[\sqrt{{\displaystyle \frac{V_0}{2}}\left({\displaystyle \frac{1+w}{1w}}\right)}{\displaystyle \frac{bt}{m_P}}\right]`$ and $`V_{\mathrm{attr}}={\displaystyle \frac{2}{b^2}}\left({\displaystyle \frac{1w}{1+w}}\right)\left({\displaystyle \frac{m_P}{t}}\right)^2`$ (6)
If $`\rho _\varphi <\rho _B`$ then $`w=w_B`$ and $`\rho _\varphi /\rho _B=3(1+w_B)/b^2=`$ const., which means $`w_\varphi =w_B0`$ and $`b^2>3(1+w_B)`$. If, however, $`\rho _\varphi >\rho _B`$ then $`w=w_\varphi `$ and $`\rho =\rho _\varphi a^{3(1+w_\varphi )}`$. This time $`b^2=3(1+w_\varphi )`$. Therefore, dark energy dominates without eternal acceleration when $`2<b^2<3(1+w_B)`$. Brief acceleration is possible at the time when $`\varphi `$ unfreezes to follow the attractor, due to superfreezing . In fact, superfreezing enlarges the range to $`2<b^2<24`$ .

## 5 A toy-model example

Consider the potential (also studied in ):

$$V(\varphi )=\frac{M^4}{\mathrm{cosh}(b\varphi /m_P)1}\{\begin{array}{cc}V2M^4e^{b\varphi /m_P}\hfill & \varphi m_P/b\hfill \\ V2M^4(m_P/b\varphi )^2\hfill & \varphi m_P/b\hfill \end{array}$$
(7)
The slow-roll parameters are: $`\eta =2A[\mathrm{cosh}(b\varphi /m_P)+2]=ϵ+2A`$. Thus, inflation ends when $`\eta (\varphi _{\mathrm{end}})1`$, where $`Ab^2(\lambda /M^4)`$. The spectral index is

$$n_s1=4A\frac{3+\left(\frac{16A}{12A}\right)\mathrm{exp}(4N_{\mathrm{dec}}A)}{1\left(\frac{16A}{12A}\right)\mathrm{exp}(4N_{\mathrm{dec}}A)}$$
(8)
where $`N_{\mathrm{dec}}69`$ is the number of e-folds of remaining inflation corresponding to the scale of the horizon at decoupling. From the observations $`|n_s1|0.1`$, which demands the constraint $`A<\frac{1}{148}`$. Using $`A1`$ the above becomes $`n_s1\frac{4}{N_{\mathrm{dec}}+1}`$, which gives $`n_s0.94`$. The density contrast with $`A1`$ is:

$$\frac{\delta \rho }{\rho }\frac{2b^2}{\sqrt{6}\pi }\left(\frac{M}{m_P}\right)^2\sqrt{A}(N_{\mathrm{dec}}+1)^2$$
(9)
Similarly, for the inflationary scale we find $`V_{\mathrm{end}}2AM^4`$ or, equivalently,

$$V_{\mathrm{end}}\frac{3\pi ^2}{b^4}\left(\frac{\delta \rho }{\rho }\right)^2\frac{m_P^4}{(N_{\mathrm{dec}}+1)^4}$$
(10)
The field freezes at $`\varphi _F=\varphi _{\mathrm{end}}+\frac{1}{\sqrt{6}}(61.74\mathrm{ln}b)m_P`$. Then coincidence demands: $`b14.5`$, which is too large for brief acceleration. Further, Eq.(10) gives $`V_{\mathrm{end}}^{1/4}2\times 10^{13}`$GeV. Using $`V_{\mathrm{end}}2b^2\lambda `$ one finds $`\lambda ^{1/4}6\times 10^{12}`$GeV and so $`M_55\times 10^{14}`$GeV. The bound on $`A`$ suggests: $`M>7\times 10^{13}`$GeV, which enables the identification $`M=M_5`$. The Hot Big Bang begins at temperature:

$$T_{}=\frac{\alpha ^3}{16b(N_{\mathrm{dec}}+1)^4}\sqrt{\frac{2g_{}}{15}}\left(\frac{\delta \rho }{\rho }\right)^2m_P$$
(11)
which satisfies the BBN constraint if $`\alpha \text{ }\stackrel{>}{}\text{ }0.1`$. With $`\alpha 0.1`$, the modified Friedman equation gives $`T_{\mathrm{reh}}10^7`$GeV, which satisfies the gravitino bound.

## 6 Conclusions

Quintessential-inflation model-building is easier in brane-cosmology because the $`\eta `$-problem is overcome by considering steep inflation . The dynamics of the Universe from inflation through to the present have been analysed and employed in the toy-model presented. This model incorporates two natural mass-scales ($`m_P`$ and $`M_5`$) and leads to a spectral index within observational bounds. However, it fails to provide late-time accelerated expansion.

Acknowledgements.  Work supported by the EU network: hprn-ct-00-00152.


text:
# Spectro-Polarimetric search for hidden AGNs in four southern Ultraluminous Infrared Galaxies

## 1 Introduction

With bolometric luminosities exceeding $`10^{12}L_{}`$, Ultra-Luminous Infrared Galaxies (ULIRGs) are the brightest sources in the local universe. Whether they are powered by a hidden quasar, a circum-nuclear starburst, or a combination of both, has been subject of a lively debate since their discovery. Although observational hints on the nature of these galaxies have been increasing in the last decade, their physical understanding is still subject to uncertainties. If mid–infrared ISO spectroscopy highlights clear cases of starburst–dominated activity (Genzel et al. 1998), hard X–ray observations otherwise suggest that the presence of an obscured AGN may be more frequent and significant than predicted on the basis of the IR data alone. Analyses of multi–wavelength (X–ray, optical and IR) observations currently suggest that both AGN and starburst activity happen concomitantly in ULIRGs (Genzel et al. 1998, Veilleux et al. 1999, Braito et al. 2002, Berta et al. 2002).

In the presence of a central engine enshrouded in a dusty torus, an appreciable non–thermal polarized radiation may be detected, produced by either electron or dust scattering. The use of polarimetry for study hidden broad line regions in Seyfert 2 galaxies has been pionereed by Miller & Antonucci (1983). Polarisation produced by scattering or absorptive dichroism may be also detected in the continuum of starburst galaxies spectra, and so the polarisation itself is not a straight indicator of an abscured AGN. The technique of searching for polarised light in IRAS galaxies has been well exploited recently by some authors (see for example Tran et al., 1999; Young et al., 1996). Up till now the only galaxies which show polarised broad H$`\alpha `$ all have seyfert or LINER classifications based on their optical spectra.

Genzel et al. (1998) studied the mid–IR spectral properties of a sample of 15 ULIRGs from a complete IRAS sample with $`S_{60\mu m}=5.4`$ Jy. The far–IR selection allows to avoid selection biases due to dust absorption, making easier the analysis of the relationship between AGN and starburst activity.

To add further constraints on the physical processes happening within these sources, we have performed a spectro-polarimetric survey of four of the Genzel et al. IRAS sample in the southern sky. We report here the results of the polarimetric observations, while the full-light spectra are presented and discussed by Berta et al (2002) and Fritz et al. (2002). A review of the optical and IR properties of these sources can be found in Duc et al. (1997). Section 2 reports on the observations and data analysis. Section 3 presents details about our results, and Section 4 contains a short discussion.

## 2 Observations and data reduction

The spectro-polarimetric data have been collected with EFOSC2 at the 3.6m ESO telescope in La Silla. Observations have been carried out during the nights of October 29th to November 1st, 2000. Observing logs are summarized in Table 1. Sky was not photometric, with a seeing between 0.9 and 1.2 arcsec during the run. We adopted the 236 line/mm grism, with an effective dispersion of 2.77 Å/pixel at the central wavelength, operating between 3690 and 9320 Å, and a 1.2$`\times `$20 arcsec slit. This slit width corresponds to a spectral resolution of 21.2 Å FWHM, or 970 Km/s in the source frame. Position angles for the slit were chosen in order to encompass the most interesting morphological features. In the case of galaxies with double-nuclei (i.e. IRAS 19254-7245 and IRAS 20100-4156), the slit has been positioned along the line joining the two nuclei.

Polarimetric observations were obtained with the combination of a $`\lambda /2`$ retarder plate and a Wollastone prism. Each object was observed with four different retarder angles (0, 22.5, 45 and 67.5), corresponding to different polarization angles $`\theta `$ of the transmitted radiation. In order to check for and eventually remove instrumental polarization, spectra of unpolarized standard stars and light transmitted through a pinhole with a polaroid filter were obtained.

Data reduction of the 2D spectra was performed using the standard IRAF tasks (for details see Berta et al. 2002). Spectra were flux-calibrated by means of the spectrophotometric standard star in NGC7293 (Oke 1990) and then corrected for atmospheric extinction as well as for Galactic reddening. The polarimetric reduction of the 1D extracted spectra was performed using software written by J.R. Walsh under the MIDAS environment and IDL procedures developed by us. Particular care has been taken in order to avoid the problem concerning the biased threshold always present in polarimetric measures. Debiasing has be done following the relation $`PP_{obs}[1(\sigma _P/P_{obs})^2]^{1/2}`$ where $`P_{obs}`$ is the observed degree of polarization, $`P`$ is the bias corrected one and $`\sigma _P`$ is the r.m.s. error on the polarization (Di Serego Alighieri, 1997).

## 3 Results

The spectra of the four observed galaxies are shown in Figures 1 and 2. Each spectrum is obtained by identifying the trace of the main nucleus. On the left panels, from top to bottom, the unpolarized and polarized spectrum, the polarization degree and direction angle are shown in the order, throughout the observed wavelength range. The panels on the right show expansions around H$`\alpha `$. The plots relative to the polarization (flux, percentage and position angle) are spectrally binned in order to have an error bar on the polarization degree of $`\pm `$0.5 % (the error bars shown in the polarization degree spectra). The original spectra are not rebinned. Error bars are shown also in the plots of polarization angle.

When detected (as in the case of IRAS 19254-7245 and IRAS 20551-4250), the polarized lines appear to contain some structure, which however is only due to the poor signal to noise ratio (we have checked in particular that these apparent blends do not correspond to the redshifted H$`\alpha `$ and \[Nii\]$`\lambda `$6584 lines). Table 2 summarizes our data on the polarization degree, H$`\alpha `$ line fluxes and observed line widths.

### 3.1 Note on individual sources

IRAS 19254–7245. the Superantennae (Colina et al., 1991), has been classified as a Seyfert 2 galaxy by Mirabel et al. (1991) on the basis of optical spectral properties. In their polarimetric analysis, Heisler et al. (1997) and Lumsden et al. (2001) fail in detecting scattered light and attribute the lack of polarized components to geometric effects: if the scattering particles lie close to the plane of the obscuring torus, an edge–on observer wouldn’t be able to observe any polarized flux. In particular Lumsden et al. (2001) found only polarised narrow emission line, and this may be explained if one consider that we were looking at a slightly different location (Lumsden, 2002). Moreover our total flux spectrum is essentially identical to theirs (Lumsden, 2002).

The analysis of our EFOSC2 integration on the Southern Nucleus leads to a detection of a significant polarized residual, in particular in the spectral region around H$`\alpha `$ (see Figure 1, top left panel). The polarization degree turns out to be $`1.8\pm 0.4\%`$.

The width of this spectral structure in polarized light, which might be interpreted as the H$`\alpha `$+\[Nii\] blend although we do not have the spectral resolution to resolve it, is $``$70 Å, or 3500$`\pm `$240 km/s. This line width is only marginally wider than that of the H$`\alpha `$+\[Nii\] blend measured in the total-light spectrum by Berta et al. (2002, 2500 km/s).

Figure 3 shows the comparison of the polarized to the total–light H$`\alpha `$ as well as the difference between the two. There is no clear positive residual, but the very low S/N doesn’t allow to be conclusive about the presence of a scattered component. In this figure (as in the following figure 4) we do not have rebinned the polarized flux spectra in order to keep the full spectral resolution around the Hα line.

We have attempted also to fit this polarized feature by taking into account some information on the complex kinematic structure of the source as reported by previously published higher-resolution work. Following Colina et al. (1991), Vanzi et al. (2002) infer that each emission lines are in fact composed of three different kinematical components. They reproduce the observed optical emission lines with a synthetic profile with three different gaussians. Similarly, we have fit our polarized emission structure with three different kinematic components for the triplet \[Nii\]$`\lambda `$6548, H$`\alpha `$ and \[Nii\]$`\lambda `$6584 (hereafter “model A”). The gaussian parameters (recession velocity $`z`$, $`\sigma =10.4`$ Åand percentage of flux) of Vanzi et al. (2002) were adopted. When taking into account our lower spectral resolution $`\sigma =9.0`$ Å, the observed width of the polarized structure is not reproduced in this way.

Instead, we have found a best-fit by adding some intrinsic line-spread, i.e. $`\sigma `$=23.3, 5.8, 4.3 Å for the three components at $`z`$=0.06149, 0.05858, and 0.06414 (see Figure 4, “model B”, Vanzi et al., 2002). Note that the sigma value of the two latter component are lower than our instrumental spectral resolution. In conclusion, our analysis indicate a kinematic component with an intrinsic width of FWHM of 55 Å, or $``$2360 km/s in the velocity space.

Given the complex kinematic structure of the source and our modest spectral resolution and low signal-to-noise, we cannot be conclusive about the presence an a broad-line AGN based on these data, particularly considering that the full-light spectrum displays a similarly broad spectral complex. Dychroic transmission of the H$`\alpha `$+\[NII\] complex through a dust screen may still be consistent with the present data.

IRAS 20100–4156. Lying at $`z=0.129583`$, this object is extremely bright ($`L_{IR}[8÷1000\mu m]510^{12}L_{}`$, adopting $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, H0=65 km/s/Mpc). The top left panel of Figure 2 shows the spectrum of the southern nucleus: only an upper limit of $`210^{17}`$ erg cm-2 s-1 to the polarized H$`\alpha `$ flux can be set. The northern nucleus, of comparable luminosity, does not show any polarization either.

IRAS 20551–4250. A detailed description of the optical properties of ESO 286-IG19 is given by Johansson (1991). We find very faint polarized flux in correspondance of the H$`\alpha `$ line at $`1\%`$ level, with low significance and with a complex structure. The polarized feature is not spectrally resolved (see Table 2). Very similar results on this object has been obtained by Lumsden and his collaborators (2002).

IRAS 22491–1808. The South America galaxy (Carico, 1990) was observed in polarized light by Young et al. (1996) in the R band, but no polarization was detected. We confirm the lack of polarization ($`<1\%`$) and set an upper limit of $`510^{17}`$ erg cm-2 s-1 to the H$`\alpha `$ scattered flux.

## 4 Discussion and conclusions

We have found very faint but significant polarization signals in two out of four of the Ultra-Luminous IR Galaxies observed in our spectro-polarimetric survey: these are the IRAS 19254-7245 (Superantennae) and IRAS 20551-4250. In the former, the polarized $`H_\alpha `$ line may be broadened by $``$2000–3000 km/s. The present observations tend to confirm that ULIRGs display, except in few cases (see, for example, Young et al. 1996, Tran et al. 1999, Hines et al. 1999), faint signals in optical polarized light. These results are consistent with either very weak AGN components in the source nuclei (otherwise dominated by starburst emission), or with a dust distribution covering the nuclear source and producing dichroic transmission. The signal-to-noise ratio of the data does not allow to distinguish between the possible polarization mechanism (e.g. electron scattering, dust scattering or dichroism). In the case of IRAS 19254-7245 a slightly trend of polarization degree with wavelength seems to be present (see Figure 1 left panel, third row) and this should exclude the presence of electron scattering (that is wavelength independent).

It is interesting to compare the present results with those based on independent data. All four of our sources have been observed with XMM-Newton by Braito et al. (2002). Their analysis is revealing clear evidence for AGN activity in IRAS 19254 and 20551. In both sources, strong iron K$`\alpha `$ lines and large hydrogen column densities, indicative of hidden AGNs, are detected. With regard to the other two sources, the presence of non–thermal activity cannot be ruled out but is strongly constrained by the XMM data.

The optical-IR-millimetric spectral energy distributions of these same sources have also been analysed by Berta et al. (2002) and Fritz et al. (2002) adopting combined stellar, starburst and AGN templates. AGN contributions to the bolometric emissions of $``$40% and 18% are indicated for IRAS 19254 and IRAS 20551, and of the order of few % for the other two. We notice interesting agreement between the results of the present polarimetric analysis and those based on entirely independent approaches.

## Acknowledgments

We wish to thank the technical staff of ESO 3.6m telescope at La Silla, and particularly M. Sterzik and R. Athreya, for their efficient support. We thanks also S.L. Lumsden for useful comments on our paper and for kindly furnish us his data on IRAS 19254-7245 and IRAS 20551-4250.