text:
# Quarkonium production and NRQCD matrix elements

## 1 INTRODUCTION

Quarkonium is an interesting system since the mass of the heavy quark is much larger than the QCD scale,

$$m_Q\mathrm{\Lambda }_{\mathrm{QCD}}.$$
(1)
One would hope, therefore, that it would be possible to calculate reliably in perturbation theory. If that were the end of the story, things would be simple. The only difficulty, which makes quarkonium very interesting, is that there is another natural small scale in the problem, the relative velocity of the heavy quarks, $`v`$. For charmonium $`v_c^21/4`$, while for bottomonium $`v_b^21/10`$.

Since we have a multi-scale problem, with the hierarchy

$$m_Qm_Qvm_Qv^2\mathrm{\Lambda }_{\mathrm{QCD}},$$
(2)
it is useful to use an effective field theory. The effective field theory appropriate for a system with two heavy quarks is Nonrelativistic QCD (NRQCD) . As the name implies, there is a relativistic expansion in $`v`$, as well as the usual perturbative expansion in $`\alpha _s`$.

Before the advent of NRQCD, quarkonium production and decay was calculated using the so-called “Color-Singlet Model” (CSM). In this model, the $`\overline{Q}Q`$ pair is treated in a color-singlet state, in the $`v0`$ limit. It was discovered, however, that the CSM did not describe quarkonium production very well in certain circumstances. For example, the CSM prediction for $`\psi ^{}`$ production at the Tevatron is off by about a factor of 30! As we will see, NRQCD can easily explain the discrepancy in this situation .

Whether NRQCD is the correct effective field theory to be used for quarkonium is an interesting and important question. There are some potential problems comparing NRQCD predictions to data, in particular the polarization. NRQCD make a definite prediction that at large transverse momenta, $`\psi ^{}`$ should be nearly 100% polarized . However, it appears that the data is not following that trend . So it may be that the $`\alpha _s`$ and $`v`$ expansions do not converge well for charmonium. This is not to say that NRQCD is incorrect. NRQCD is a valid effective field theory, in that as the mass $`m_Q`$ goes to infinity, the theory correctly reproduces full QCD. It may be that $`v_c`$ is too large, or $`m_Qv_c^2`$ too small, for NRQCD to be useful in describing charmonium production. One could still hope that it will work for bottomonium.

How is it possible to test NRQCD? There is the polarization discussed above . Another way is to compare NRQCD matrix elements extracted from different experiments.

## 2 FACTORIZATION FORMALISM

The physical picture of production in NRQCD begins with a hard scattering, in which a $`\overline{Q}Q`$ pair are produced with any spin, angular and color quantum numbers. This process can be calculated in perturbation theory. Then the $`\overline{Q}Q`$ evolve into the final state quarkonium $`H`$ in some non-perturbative fashion. This is encoded in the NRQCD matrix elements, which must be extracted from experiment.

The NRQCD factorization formula is the mathematical realization of the above picture. A general production process for $`i+jH+X`$ can be written as

$$d\sigma =\underset{n}{}d\widehat{\sigma }(ijQ\overline{Q}[n]+X)O^H(n),$$
(3)
where $`d\widehat{\sigma }(ijQ\overline{Q}[n]+X)`$ is the cross section for producing $`Q\overline{Q}`$ in state $`n`$ by scattering $`i`$ and $`j`$, calculable in a perturbative expansion in $`\alpha _s`$ and perhaps convoluted with parton distribution functions (PDFs). The long distance matrix elements encode the hadronization of the heavy quarks in state $`n`$ into the final quarkonium state $`H`$. The matrix elements can be written as

$`O^H(n)`$ $`=`$ $`0|\psi ^{}\mathrm{\Gamma }^n\chi |{\displaystyle \underset{X}{}}H+X`$ (4)
$`\times {\displaystyle \underset{X}{}}H+X|\chi ^{}\mathrm{\Gamma }_{}^{n}{}_{}{}^{}\psi |0`$
The $`\mathrm{\Gamma }^n`$ contains Dirac matrices, color matrices and derivatives. The NRQCD matrix elements scale with a definite power of $`v`$, determined by $`\mathrm{\Gamma }^n`$, allowing the truncation of the $`v`$ expansion. Again, we expect the velocity expansion will work better in the bottom sector than the charm sector. The factorization described above still needs to be put on firmer ground. We should keep in mind that when we use the NRQCD factorization formalism to predict production rates, we are not only testing the $`\alpha _s`$ and $`v`$ expansions, but also the factorization formalism.

The scaling is determined by looking at what perturbations are necessary to give a non-vanishing result for the time ordered product. To give non-zero overlap, multipole moment interactions may need to be inserted. For example, the matrix element $`O_8^{J/\psi }(^3S_1)`$ scales as $`v^4`$ compared to the $`O_1^{J/\psi }(^3S_1)`$, since we need two $`E1`$ insertions in the amplitude. The first insertion changes $`L`$ and neutralizes the color. The second changes $`L`$ back to an $`S`$ wave. Each insertion cost a factor of $`v`$ in the amplitude, or a total of $`v^4`$ in the rate. The scalings of the most important matrix elements for $`\psi `$ production are

$`O_1^\psi (^3S_1)`$ $``$ $`v^0,`$ (5)
$`O_8^\psi (^3S_1)`$ $``$ $`v^4,`$ (6)
$`O_8^\psi (^3P_J)`$ $``$ $`v^4,`$ (7)
$`O_8^\psi (^1S_0)`$ $``$ $`v^4.`$ (8)
We can relate the color-singlet matrix elements to the wave function at the origin, or derivatives of the wave function at the origin, using the vacuum saturation approximation

$$O^H(n)0|\psi ^{}\mathrm{\Gamma }^n\chi |HH|\chi ^{}\mathrm{\Gamma }^n\psi |0,$$
(9)
dropping the sum over $`X`$. This allows us to use the lattice, data, or models to obtain the color-singlet matrix elements. It also means that the color-singlet model has been incorporated into the NRQCD factorization formalism.

A quick aside. Many people say the NRQCD model, or the “color-octet” model. However, NRQCD is not a model. The CSM, on the other hand, is a model, since there is no limit in which it will reproduce full QCD.

## 3 QUARKONIUM AT THE TEVATRON

As previously mentioned, the CSM fails to describe the data at the Tevatron. We will now discuss how NRQCD improves the situation. We will concentrate on $`J/\psi `$ production. A similar analysis can be done for $`\mathrm{{\rm Y}}`$ production. See Refs.  for recent analyses.

At large transverse momentum, the leading order NRQCD prediction in $`\alpha _s`$ and $`v`$ is $`O(\alpha _s^3v^0)`$. This is the same as the leading order CSM prediction, which grossly underestimates the rate. Even the shape of spectrum is wrong, with the theory dropping as $`1/p_{}^8`$ compared to the data which falls as $`1/p_{}^4`$.

We can go to higher order in the perturbative expansion: $`O(\alpha _s^5v^0)`$. In this case, fragmentation processes occur, which, due to some propagators being close to on-shell, can have large contributions. These fragmentation diagrams dominate over the lower order graphs at large $`p_{}`$, with the theory now falling as $`1/p_{}^4`$, similar to the data. However, the normalization of the theory is still far below the experiment.

This is as far as the CSM can go. To improve the prediction, NRQCD is needed . Instead of just including higher order in $`\alpha _s`$, we can also include higher order in $`v`$: $`O(\alpha _s^3v^4)`$. Now there are color-octet fragmentation processes proportional to the matrix element $`O_8^\psi (^3S_1)`$, which also scales as $`1/p_{}^4`$. Of course, the value of the matrix element is unknown, so the fact that we can get a good fit really means that the extracted matrix element is not abnormally large or small.

The fragmentation contribution can be approximated as

$$d\sigma (p\overline{p}H)=\underset{i}{}𝑑z𝑑\sigma _iD_{iH},$$
(10)
where $`d\sigma _i`$ is the cross section to produce an on-shell parton $`i`$, $`D_{iH}`$ is the fragmentation function for parton $`i`$ to produce quarkonium state $`H`$ with momentum fraction $`z`$, and the sum is over partons. The advantage of writing it in this way is that it is easy to sum large logs of $`\mathrm{log}p_{}/m_Q`$ using Altarelli-Parisi (AP) evolution.

At lower $`p_{}`$, fragmentation no longer dominates, and it is necessary to include all non-fragmentation diagrams as well. Contributions due to other octet matrix elements, $`O_8^\psi (^1S_0)`$ and $`O_8^\psi (^3P_J)`$ are now important. These channels both fall as $`1/p_{}^6`$, making it impossible to determine them independently. Instead, a linear combination of the two

$$M_k^H=O_8^\psi (^1S_0)+\frac{k}{m_Q^2}O_8^\psi (^3P_0)$$
(11)
is usually extracted.

At small $`p_{}`$, higher twist contributions become important. These contributions scale as powers of $`\mathrm{\Lambda }/\sqrt{p_{}^2+m_Q^2v^n}`$, where the value of $`n`$ is not known. Also at small $`p_{}`$, the rates tend to diverge, due to soft gluon effects. To handle these divergences, the low $`p_{}`$ region can just be ignored or modeled by including initial parton showering/intrinsic $`k_{}`$.

In Table 1, different extractions of $`M_k^{J/\psi }`$ are collected. The second set of errors, when present, are due to scale variation. The other errors are statistical. As can be seen, there is large uncertainty due to the choice of PDF. At this point, the best we can hope to say is that this combination of matrix elements is on the order of

$$M_3^{J/\psi }\mathrm{few}10^2\mathrm{GeV}^3.$$
(12)
Table 2 contains different extractions of $`O_8^{J/\psi }(^3S_1)`$. Again the second set of errors, when present, are due to scale variation, while the other errors are statistical. Here the largest error is not due to the PDFs, but due to the scale variation. This is an indication that higher order perturbative corrections may be very large. Again we only know the matrix element is on the order of

$$O_8^{J/\psi }(^3S_1)\mathrm{few}10^3\mathrm{GeV}^3.$$
(13)
## 4 LEP EXTRACTION

We would like to do better than the order of magnitude extractions discussed above at the Tevatron. Since the dominant errors are due to scale or PDF variation, it seems difficult make improvements at the Tevatron. We would also like to test the formalism, by comparing the matrix elements extracted from different experiments.

The natural place to look for cleaner extractions would be at an $`e^+e^{}`$ machine, since there is no PDF nor initial state gluon radiation to worry about. The prediction at CLEO does not have a large dependence on the color-octet matrix elements away from the endpoint. Therefore LEP is the natural choice.

Formally there are two leading order contributions in the $`\alpha _s`$ and $`v`$ expansion, both in the singlet channel, of order $`O(\alpha _s^2v^0)`$. The contribution from gluon radiation in the singlet channel $`Z\psi gg`$ is suppressed by powers of $`M_\psi ^2/E_\psi ^2`$ . The color-singlet charm quark fragmentation process $`Z\psi c\overline{c}`$ , which has no power suppression, dominates over non-fragmentation processes for large $`E_\psi `$. Light quark octet fragmentation (in which the mother parton does not combine to form part of the bound state) is naively of order $`\alpha _s^2v^4`$, down by $`v^41/10`$ compared to charm fragmentation. However, this channel is enhanced due to the presence of large logs and a numerical factor of five due to the number of possible quarks that initiate the process . The same logs that enhance the octet channel also put the convergence of the perturbative expansion into question.

The tree-level calculation of the differential cross section in the color octet production channel  scales as $`z1`$ as

$$\frac{d\mathrm{\Gamma }}{dz}\alpha _s^2\frac{\mathrm{log}(M_Z^2/M_\psi ^2)}{z}O_8^{J/\psi }(^3S_1),$$
(14)
leading to large double logs in the total rate. Since $`\alpha _s\mathrm{log}(M_Z^2/M_\psi ^2)1.5`$, we should treat $`\alpha _s\mathrm{log}(M_Z^2/M_\psi ^2)`$ as order one and resum all powers of the large logarithm. With this counting, the octet channel is $`O(\alpha _s^0v^4)`$, on par with the singlet fragmentation contribution. More practically, the tree-level calculation has a factor of two uncertainty associated with the scale at which $`\alpha _s`$ is evaluated, since $`\alpha _s(M_\psi )/\alpha _s(M_Z)2`$ (this is just a restatement that there is a large logarithm). The resummation of the leading logarithms reduces this uncertainty, so the resummation procedure is essential from both a practical and a formal standpoint.

To resum the logs, we write the rate in the fragmentation limit as

$`{\displaystyle \frac{d\mathrm{\Gamma }}{dz}}`$ $`=`$ $`{\displaystyle _z^1}{\displaystyle \frac{dy}{y}}[2C_q(\mu ^2,y)D_q(\mu ^2,z/y)`$ (15)
$`+C_g(\mu ^2,y)D_g(\mu ^2,z/y)].`$
We can now use AP to sum the logs of $`M_Z^2/M_\psi ^2`$.

However, summing the above mentioned logs will only yield the correct leading order differential rate if $`z`$ is sufficiently large.111At very large $`z1`$ we do not get the correct result since we have neglected the shape function which enters near the endpoint. This shape function has been argued to be the solution of the so called “Hera anomaly”. It is not important here since the rate is very small in this region. When $`z`$ is parametrically small, terms of the form $`\alpha _s\mathrm{log}(z)/z`$ become just as important. Furthermore, these logs will also contribute double logs to the total rate given that the lower limit on $`z`$ is $`2M_\psi /M_Z`$. This second type of log, due to soft gluon emission, is resummed using a formalism familiar from discussions of jet multiplicities , where the color-coherence of the soft gluon emission is very important.

It is possible to do both of these resummations , see Fig. 1. The rate depends on the linear combination

$`{\displaystyle \underset{m}{}}O_8^{H(m)}(^3S_1)\times \mathrm{Br}[H(m)J/\psi X]`$
$`=(1.9\pm \mathrm{\hspace{0.17em}0.5}_{stat}\pm \mathrm{\hspace{0.17em}1.0}_{theory})\times 10^2\mathrm{GeV}^3`$ (16)
since the data includes feeddown from higher charmonium states. It turns out that this leads to an extraction with much smaller theoretical uncertainty than the analogous value from the Tevatron $`1.4\pm 0.2_{stat}\pm 1.4_{theory}`$. For the Tevatron value, the extraction from Ref.  was used for the central value.

It is difficult to tell from Fig. 1 if there is much improvement in the prediction from this resummation. To make the comparison more quantitative, it is useful to take the ratio of the first moment to the zeroth moment

$$\frac{1}{\mathrm{\Gamma }(ZJ/\psi X)}𝑑zz\frac{d\mathrm{\Gamma }}{dz}(ZJ/\psi X)$$
(17)
which for the resummed rate is $`0.30`$. The tree level rate gives $`0.5`$. A very rough estimate of this quantity obtained from the data  suggests a value of $`0.26\pm 0.10`$. This is in sharp contrast to the color singlet prediction. The tree-level color singlet decay rate predicts the ratio of the first moment over the zeroth moment to be $`0.62`$. Resummation softens the color singlet decay rate, but the ratio is still too large, $`0.47`$. The ratio is independent of the color singlet matrix element. A rigorous extraction of the first moment by the experimental groups could provide an extremely clean, quantitative test of the NRQCD approach.

## 5 OTHER EXTRACTIONS

### 5.1 HERA

The NRQCD prediction does not seem to fit the data at HERA very well. As $`z=E_{J/\psi }/E_\gamma `$ goes to the endpoint, $`z1`$, the NRQCD prediction begins to rise. This has been dubbed the “HERA anomaly”, and has been used to argue that NRQCD is incorrect.

However, it is now known that the NRQCD prediction breaks down near the boundary of phase space, precisely where there is a problem at HERA . In this region the velocity expansion is breaking down, and a shape function must be introduced. In fact, one should not compare the prediction with the data above $`z1v^20.7`$. If one looks only below this point, there is no problem with the theory. Unfortunately, it is not really possible to extract the matrix elements from the data in this region.

### 5.2 $`BJ/\psi `$ DECAYS

For $`BJ/\psi `$ decays, the first worry is that the corrections to factorization will be large since there is not a lot of energy in the decay products. There are also $`\alpha _s`$ corrections  and corrections in the Heavy Quark Effective Theory $`1/m_b`$ expansion , the most important being the Fermi motion of the b-quark, which requires the use of another shape function.

Nevertheless, it is possible to obtain the linear combination $`M_k^{J/\psi }`$ from $`B`$ decays, with qualitative agreement with the Tevatron extractions

$`M_{3.1}^{J/\psi }=1.5\times 10^2\mathrm{GeV}^3,`$ (18)
$`M_{3.4}^{J/\psi }=2.4\times 10^2\mathrm{GeV}^3.`$ (19)
## 6 CONCLUSIONS

NRQCD has been used for a number of years to make predictions about quarkonium production. At this time, there are consistent extractions of the NRQCD matrix elements from a number of different experiments. The most analysis has been done on hadronic collisions, in particular at the Tevatron due to the high statistics.

LEP is a very clean place to extract the matrix elements. It is possible obtain smaller theoretical uncertainty compared to the Tevatron extractions. Unfortunately the statistics are not as good. Nevertheless, due to the large theoretical errors in the Tevatron extractions, the LEP extraction is the best extraction to date.

It is also possible to test the NRQCD formalism at LEP by looking at the moments of the spectrum. Without the color octet channels included in the NRQCD formalism, the theory and experiment would not agree for this test. This is one place where NRQCD passes a test. It will be interesting to see whether it will be possible to resolve the polarization problem within the NRQCD factorization formalism.

I would like to thank my collaborators, Glenn Boyd, Peter Cho, and Ira Rothstein, and the organizers of the IV International Conference on Hyperon, Charm and Beauty Hadrons. This work was supported in part by the Department of Energy under grant number DOE-ER-40682-143.


format: text

text:
# Particle number density fluctuations and pressure effects on structure formation

## I Introduction

According to current cosmological models, the large-scale structure, as seen today, formed from the evolution of gravitational instabilities which originated from inflationary processes in the primordial fluid. In studying the evolution of such instabilities, analytically in the fluid approximation***For high accuracy the full Boltzmann equations are used and these require numerical codes, which provide different insights into the process., one usually assumes a dust or radiation energy momentum tensor. The object of this paper is to consider some implications of modifying these assumptions.

First we note that while the two extreme fluid states are reasonable for most of the history of the universe, there are short, but at the same time important, specific periods in which matter is in the non-relativistic regime with a tiny although non-negligible pressure. These are periods during which collisions, or small random motions (velocity dispersion) of collisionless matter, give rise to a small but non-zero matter pressure. Even though in these cases the matter pressure in the background may be negligible, small disturbances (perturbations) of the spacetime induce a tiny pressure term which can modify the evolution of the inhomogeneities. Such a scenario may also have applications when studying weakly self-interacting matter (for details see Spergel and Steinhardt ) or warm dark matter (for more details of this approximation see Hogan and Dalcanton ) in the linear regime. These dark matter models have been the subject of much recent work using numerical integrations in the nonlinear regime, where they have been employed to investigate ways of resolving problems with the standard CDM model on small scales (see ) in galaxy formation.

Second we take the particulate nature of matter seriously. It must be born in mind that in structure formation one really needs to look at the clumping of a set of particles rather than at the evolution of energy density instabilities. Although in taking the dust assumption to be valid such a distinction is not relevant, this is not true for those periods in which matter is undergoing a transition to the non-relativistic regime and, consequently, for which the dust approximation is not accurate enough, even in the linear regime. To model the effects of the particulate nature of the matter we make use of thermodynamic arguments. An alternative Kinetic Theory approach is given in .

The method we use is based on the Ellis-Bruni  covariant and gauge-invariant perturbation formalism, which has proved useful in different cosmological problems , , . It is one of the gauge-independent descriptions of relativistic perturbations, described in detail in the literature, see for instance ,,,.

In our treatment of the fluid models we generalise the usual equations for the growth of cosmological perturbations for matter in a Friedman-Robertson-Walker (FRW) background to include first order pressure effects due to a non-relativistic component, in the radiation and in the matter dominated eras. Thus the simplifying hypothesis that the matter is dust, i.e. has vanishing pressure, is relaxed, in both cases, to include a small pressure due to collisions or thermal motions of the particles which make up the fluid.

## II Matter versus energy density perturbations

The starting point for the the Ellis-Bruni  formalism is the choice of a time-like congruence $`u^a`$ along which typical (comoving) observers move. Once this four-velocity has been selected all the physical quantities are split into their spatial and time-like parts using the spatial projector tensor $`h^{ab}`$. The relevant physical quantities are defined as follows

$$\delta _a=a\frac{\text{D}_a\rho }{\rho },\theta _a=a\text{D}_a\theta ,p_a=a\frac{\text{D}_ap}{\rho },\nu _a=a\frac{\text{D}_an}{n},$$
(1)
where $`\text{D}_ah_a^b_b`$ is the spatial projection of the covariant derivative, $`\theta u_{;a}^a`$ the expansion, $`p`$ the pressure, $`\rho `$ the energy density and $`a`$ is a representative or average length scale defined by $`3\dot{a}/a=\theta `$. In a FRW metric $`a`$ is the background scale factor and we will use the notation $`3H=\theta `$ in the background. The quantities $`\delta _a`$, $`p_a`$, $`\theta _a`$ and $`\nu _a`$ describe the spatial inhomogeneity of quantities which are homogeneous and isotropic in the background. Each of them is spatial (contraction with $`u^a`$ vanishes), covariant and gauge-invariant to first order, as all of them vanish in the FRW background . The scalar parts of perturbations described in (1) are obtained by taking the comoving spatial divergence, e.g., $`\delta =a\text{D}^a\delta _a`$.

The general procedure to obtain evolution and constraint equations for the quantities defined in (1) starts from a suitable spatial projection of the Ricci and Bianchi identities as given in terms of the fluid variables once the field equations have been taken into account (see for example  for a pedagogical account). Here we shall follow the notation of .

Restricting ourselves to a spatially flat FRW background, the equations governing the evolution of the comoving fractional density gradient $`\delta _a`$ and the comoving expansion gradient $`\theta _a`$ are (see  equations 25,26)

$`\dot{\delta }_a=3wH\delta _a(1+w)\theta _a,`$ (2)
$`\dot{\theta }_a{\displaystyle \frac{9}{2}}H^2p_a+2H\theta _aa\text{D}_a\text{D}^b\dot{u}_b+{\displaystyle \frac{3H^2}{2}}(\rho _a+3p_a)=0,`$ (3)
where

$$w=\frac{p}{\rho }.$$
(4)
Note that equation (2) omits second order terms proportional to $`\sigma _{ab}\delta ^b`$ and $`\omega _{ab}\delta ^b`$ where $`\sigma _{ab}`$ and $`\omega _{ab}`$ are the shear and vorticity tensors.

We will now discuss perturbations in the particle number density. This discussion is motivated by the idea that a full understanding of large scale structure formation in the framework of gravitational instabilities requires an analysis of how particles, which comprise the cosmic fluid, clump together due to their mutual attraction. Although dealing with the energy density is the usual procedure, we believe that it is desirable to consider the particle number density $`n`$ because, ultimately, particles are the building blocks of all the cosmic structures we see today. For dust, it does not make any difference whether one chooses $`\rho `$ or $`n`$ because no random motions are present and the internal energy vanishes identically, yielding $`\rho =mn`$, i.e.

$$\delta _a=\nu _a.$$
(5)
The situation is quite different when internal energy associated with thermal motions is introduced. The behaviour of $`\nu _a`$ requires a separate investigation. In  it is shown that from the Gibbs equation

$$T\text{d}s=\text{d}\left(\frac{\rho }{n}\right)+p\text{d}\left(\frac{1}{n}\right),$$
(6)
it follows that

$$e_a=\delta _a(1+w)\nu _a,$$
(7)
where

$$e_a=\frac{anT\text{D}_as}{\rho },$$
(8)
is a dimensionless entropy gradient. Perturbations which are non-dissipative ($`\dot{s}=0`$) and without spatial variation of the entropy ($`e_a=0`$) are called isentropic in  because they have the same entropy at all the points of the spacetime. In the weaker case where no spatial variation of the entropy is allowed $`(e_a=0)`$ the number density perturbations are related algebraically to the energy density perturbations by

$$\nu _a=\frac{\delta _a}{1+w}.$$
(9)
In fact, $`e_a=0`$ is a strong assumption which does not hold in general. Its departure from zero affects the growth of perturbations. From the particle conservation equation

$$\dot{n}+n\theta =0,$$
(10)
we get an evolution equation for $`\nu _a`$

$$\dot{\nu }_a+\theta _a\frac{\theta }{1+w}p_a=0.$$
(11)
If we use the standard, but generally unphysical, equation of state $`w=\mathrm{constant}`$ then together with (2), the integral of (11) yields

$$\nu _a=\frac{\delta _a}{1+w}+k_a,\dot{k}_a=0.$$
(12)
Thus, in general, when this simple equation of state is used, a new stationary (entropy) mode appears in the evolution of $`\nu _a`$ which could have an impact on the structure formation at certain stages of the evolution of inhomogeneities. The problem with (12) is that it makes sense only for the very special, but widely used, cases of dust and radiation. In more general circumstances such as in self interacting dark matter or where there are collisions , the equations of state should involve two independent thermodynamical quantities, i.e. $`n`$ and $`\rho `$. To deal with the inhomogeneities that arise in non-relativistic matter with a non-vanishing pressure we may follow two different approaches.

(i) If we assume that the fluid has reached a collision-dominated equilibriumFor an example of such a model see the paper by Spergel and Steinhardt  where it is pointed out that under certain circumstances the cold dark matter behaves as a collision dominated gas , we learn from kinetic theory that the equation of state in the non-relativistic regime is

$$\rho =mn+\frac{3}{2}p,$$
(13)
which together with the energy balance equation

$$\dot{\rho }+\theta (\rho +p)=0,$$
(14)
and the particle number balance equation (10) leads to

$$\dot{p}=\frac{5}{3}\theta p,$$
(15)
with solution

$$p=p_0\left(\frac{a_0}{a}\right)^5.$$
(16)
Also from (13) we get the equation

$$\delta _a=\frac{mn}{\rho }\nu _a+\frac{3}{2}p_a$$
(17)
which, together with (11) and the conservation equations, determines $`\nu _a`$.

(ii) Without imposing the restriction of a collision-dominated equilibrium fluid, we could tackle the problem of a collisionless gas in which a small pressure arises because of velocity dispersion (random motions) of the particles making up the systemSuch a dark matter model has recently been described by Hogan and Dalcanton . In  it was shown that by neglecting higher powers of the velocity dispersion in the Boltzmann equation, the evolution equation for pressure is given by (15). We remark that in obtaining this equation no assumption about the equation of state was made. Only the absence of collisions and the smallness of the velocity dispersion were used. Despite the evident mathematical similarities to case (i), a crucial difference emerges here: for a collisionless gas the Gibbs equation (6) does not apply and, as we shall show below, this leads to some new results.

In the next section we discuss cosmological settings where the small pressure terms introduced above are included.

## III Non-relativistic matter perturbations with pressure

In this section we investigate the role played by a first order pressure in the non-relativistic regime on an FRW background filled, in the first case, by a dust gas and, in the second, by dust and radiation. The point is that while the idealized background can be described by dust or dust and radiation, as soon as the dust is disturbed (by velocity dispersion or collisions) a small nonzero pressure arises. This means that the pressure is a first order quantity, as it vanishes in the background. This differs from the usual approach which takes an identically zero pressure in the background and in the real (perturbed) spacetime. We relax this assumption by allowing first order corrections to pressure due either to the fluid being in a collision-dominated regime (case (i) above) or through velocity dispersion without collisions (case (ii) above). This enables us to initiate a theoretical discussion of the effects being studied mainly qualitatively or numerically in the papers . One-component and two-component models will be considered in turn. We note here that the perturbations will involve perturbations of the entropy. The precise relation of the covariant description of entropy perturbations to the metric description  is not straightforward and will not be discussed here.

### A Dust background

For simplicity we use an Einstein-de Sitter background. As stated above we focus on perturbations for which pressure is a first order quantity. From the previous section we have that, both in case (i) and case (ii), $`p`$ is given by the evolution equation (15) leading to (16)

$$p=p_0\left(\frac{a_0}{a}\right)^5.$$
(18)
As $`p`$ is now a first order quantity we get from (2) and (3),

$$\ddot{\delta }_a2H\dot{\delta }_a+\frac{3}{2}H^2\delta _a+\frac{a}{\rho }\text{D}_a\text{D}^2p=0,$$
(19)
where the momentum balance equation

$$(\rho +p)\dot{u}_a=\text{D}_ap,$$
(20)
has been used. Note that for large-scale perturbations the Laplacian term is negligible and we recover the standard evolution equation for dust. Hence, as expected, pressure plays a role only at small scales.

### B Collision dominated fluid

At this point we distinguish between the two cases discussed in section two. If we assume that the fluid is in collision-dominated equilibrium, then the Gibbs equation holds which, together with the equation of state (13), leads to

$$\text{D}_ap=0,$$
(21)
at first order, reducing (19) to the standard equation valid for dust at any scale. Expression (21) implies that the spatial gradient of $`p`$ is a gauge-invariant second order quantity since it vanishes at first order. This allows us to go beyond first order and derive a second order equation for the pressure as we will now show. A fully covariant and gauge-invariant treatment for second order perturbations is not possible at the moment. The point is that in order to construct $`n`$-order gauge-invariant variables following Ellis-Bruni, we have to ensure that they are constant (usually zero) at all orders below $`n`$ . In particular it is not clear how to find (if at all possible) a second-order meaningful variable for the perturbed density in such a way that it vanishes at zero and first orders. Bruni and co-workers provide a systematic method to tackle relativistic perturbations beyond the linear order, but one has to pay the price of losing the gauge-invariant character .

It is possible, however to obtain a second order gauge invariant constraint on the perturbations from the above results. Operate on equation (15) with $`a\text{D}_a/\rho `$ and use the identity (26)§§§Since $`p_a`$ is second order, the first order identity gives rise to second order terms and any corrections are of higher order. below to obtain

$$\dot{p}_a+p_a\frac{\dot{\rho }}{\rho }+\frac{5}{3}\theta _a\frac{p}{\rho }+\frac{5}{3}\theta p_a=0$$
(22)
and then use,

$`\dot{\delta }_a`$ $`=`$ $`\theta _a+O(2)`$
$`\dot{\rho }`$ $`=`$ $`3\left({\displaystyle \frac{\dot{a}}{a}}\right)p+O(1)`$
derived from equations (2) and (14) and where $`O(n)`$ indicates quantities of order $`n`$. Note that $`p_a`$ is of order $`O(2)`$. After a change in the independent variable the equations (22) become

$$p_a^{}+2\left(\frac{a_0}{a}\right)p_a=\frac{5p}{3\rho }\delta _a^{},$$
(23)
where a prime denotes a derivative with respect to $`a/a_0`$. According to the argument above and the fact that $`p`$ and $`\delta _a`$ are first order quantities, equation (23) contains only second order quantities. It is one of the set of second order equations. We treat it as a constraint on the gauge invariant pressure perturbations in the sense that a full solution to all the second order perturbation equations will have to satisfy (23) although the actual solution may be more restricted. With this in mind we proceed to solve it on a dust background. In equation (23) we use equations (16) and (29), below, for the dust background and the standard first order solution for dust

$$\delta _a=K_a^+\left(\frac{a}{a_0}\right)+K_a^{}\left(\frac{a}{a_0}\right)^{3/2},$$
(24)
where $`K_a^+`$ and $`K_a^{}`$ are constants of order $`O(1)`$, to obtain

$$p_a=A_1K_a^+\left(\frac{a_0}{a}\right)+A_2K_a^{}\left(\frac{a_0}{a}\right)^{7/2}+K_a^{}\left(\frac{a_0}{a}\right)^2,A_1\frac{5p_0}{3\rho _0a_0}A_2a_0^{5/2}.$$
(25)
Here $`A_1`$, $`A_2`$, $`K_a^+`$ and $`K_a^{}`$ are constants of order $`O(1)`$, so the products are of order $`O(2)`$ and $`K_a^{}`$ are constants of order $`O(2)`$. The solution (25) is a simple and exact solution giving a limitation on the evolution of a second order quantity in the Ellis-Bruni formalism. The point is that we can apply the formalism because $`\text{D}_ap`$ vanishes at the zero and at the first order approximation.

### C Collision free fluid

The more interesting and physically appealing situation is that of a collisionless gas for which the foregoing arguments are not applicable. Applying the operator $`a\text{D}_a\text{D}^2`$ to (15) and using the first order identity

$$(a\text{D}_af)^{}=a\text{D}_a\dot{f},$$
(26)
where $`f`$ denotes any first order scalar or tensor quantity, we get

$$(a\text{D}_a\text{D}^2p)^{^{}}=7\left(\frac{a_0}{a}\right)a\text{D}_a\text{D}^2p,$$
(27)
and hence

$$a\text{D}_a\text{D}^2p=A_a^0\left(\frac{a_0}{a}\right)^7,A_a^0a\text{D}_a\text{D}^2p|_{a=a_0},\dot{A}_a^0=0.$$
(28)
For a dust background

$$H=\sqrt{\frac{\rho _0}{3}}\left(\frac{a_0}{a}\right)^{3/2},\rho =\rho _0\left(\frac{a_0}{a}\right)^3,$$
(29)
so equation (19) reduces to

$$\delta _a^{^{\prime \prime }}+\left(\frac{H^{^{}}}{H}+3\frac{a_0}{a}\right)\delta _a^{^{}}\frac{3}{2}\left(\frac{a_0}{a}\right)^2\delta _a=\frac{A_a^0}{\rho _0H^2}\left(\frac{a_0}{a}\right)^6,$$
(30)
which can be solved to give

$$\delta _a=K_a^+\left(\frac{a}{a_0}\right)+K_a^{}\left(\frac{a_0}{a}\right)^{3/2}\frac{3A_a^0}{\rho _0^2}\left(\frac{a_0}{a}\right)$$
(31)
where $`\dot{K}_a^\pm =0`$. Thus a new non-adiabatic mode is obtained which decays, albeit at a slower rate than the standard decaying mode for dust. Although the effect of the mode will be at small scales, as desired, the fact that it is decaying means that it is unlikely to be significant.

### D Dust-radiation background

The novelty in this subsection is that the dynamics of the background is given by a decoupled mixture of dust and radiation. We will consider perturbations of the matter component which again will be taken as dust in the background but with a non-vanishing first order pressure in the real spacetime. This situation mimics that of the cosmological fluid around the time of the transition from the radiation to the matter dominated era with a matter component which is acquiring a major role in the dynamics. In the following $`p_r`$, $`\rho _r`$ and $`\rho _m`$ denote the zero order pressure and energy density of radiation and matter respectively, whereas $`p_m`$ is the first order pressure of matter as given by (16).

Assuming that the fluids are decoupled and share the same 4-velocityUsing the same velocity is valid when the two fluids are coupled. Here it is a simplifying assumption which enables us to provide an analytical treatment of a problem, that has only been treated numerically up to now, and we are able to obtain new results. The role of the assumption is made explicit, in particular in the comment following equation (40), where it is relevant. we get

$`(\rho _r+p_r)\dot{u}_a+\text{D}_ap_r=0,`$ (32)
$`\rho _m\dot{u}_a+\text{D}_ap_m=0.`$ (33)
From these equations and the equation of state for radiation we find that

$$\frac{\text{D}_ap_r}{\rho _r+p_r}=\frac{\text{D}_ap_m}{\rho _m}\delta _a^r=4p_a^m,$$
(34)
where

$$\delta _a^ra\frac{\text{D}_a\rho _r}{\rho _r},p_a^ma\frac{\text{D}_ap_m}{\rho _m}.$$
(35)
Note that if we were to keep the matter component as dust both in the background and in the real spacetime, we would have ended up with $`\dot{u}_a=0=\text{D}_ap_r`$, i.e. radiation perturbations would not have been allowed.

In order to get an evolution equation for $`\theta _a`$ we start with the Raychaudhuri equation

$$\dot{\theta }+\frac{1}{3}\theta ^2\text{D}^a\dot{u}_a+\frac{1}{2}(\rho +3p)=0,$$
(36)
where $`\rho =\rho _r+\rho _m`$ is the total energy density and $`p=p_r+p_m`$ is the total pressure. We have to be very careful in dealing with $`p`$ since in the background $`p=p_r`$ but when we apply the spatial derivative operator we get $`\text{D}^ap=\text{D}^ap_r+\text{D}^ap_m`$. Bearing all this in mind the evolution equation for $`\theta _a`$ becomes

$$\dot{\theta }_a+2H\theta _a=a\text{D}_a\text{D}^b\dot{u}_b\frac{9}{2}H^2a\dot{u}_a(1+w)\frac{1}{2}(8\rho _rp_a^m+3\rho _mp_a^m+\rho _m\delta _a^m),$$
(37)
where equation (34) has been used.

Before dealing with the acceleration term we recall the identity

$$\text{D}^2(\text{D}_af)=\text{D}_a(\text{D}^2f)+\frac{2}{3}(\rho 3H^2)\text{D}_af+2\dot{f}\text{curl}\omega _a.$$
(38)
From this identity and bearing in mind that $`p_m`$ is a first order quantity, we have

$$a\text{D}_a\text{D}^b\dot{u}_b=a\frac{\text{D}_a\text{D}^2p_m}{\rho _m}=a\frac{\text{D}^2\text{D}_ap_m}{\rho _m}=\text{D}^2p_a^m.$$
(39)
It is important to note that the result obtained from commuting the operators $`\text{D}_a`$ and $`\text{D}^2`$ would have been different if we had used $`p_r`$ instead of $`p_m`$. In the latter case a first order term $`\dot{\rho }_r\text{curl}\omega _a`$ would have arisen. This means that consistency requires

$$\text{curl}\omega _a=0,$$
(40)
i.e., in a radiation-matter decoupled mixture for which (a) both fluids share the same 4-velocity and (b) the matter pressure arises from perturbations so that $`p_m`$ is a first order quantity, $`\mathrm{curl}\omega _a`$ vanishes at first order.

Changing the independent variables from $`t`$ to $`(a/a_0)`$, we get the evolution equation for $`\theta _a`$,

$$\theta _a^{}+2\theta _a\left(\frac{a_0}{a}\right)=\left(\frac{a_0}{a}\right)\left\{\frac{9H}{2}\left(4+\frac{3\rho _m}{\rho _r}\right)p_a^m\frac{\text{D}^2p_a^m}{H}\frac{1}{2H}\left[p_a^m(8\rho _r+3\rho _m)+\rho _m\delta _a^m\right]\right\}.$$
(41)
The background equations for dust plus radiation are

$$\rho =\rho _m+\rho _r=\frac{3}{\beta a^4}(1+a\alpha ),p=p_r=\frac{1}{\beta a^4},H^2=\frac{1}{\beta a^4}(1+a\alpha ),$$
(42)
where

$$\beta \frac{3}{\rho _0^ra_0^4},\alpha \frac{\rho _0^m}{\rho _0^ra_0}.$$
(43)
From the evolution equation for $`p`$ for a collisionless gas

$$\dot{p}_m=\frac{5}{3}\theta p_m,$$
(44)
and using the equation (26) we have

$$\dot{p}_a^m=2Hp_a^m(p_a^m)^{}=\frac{2}{a}p_a^m,$$
(45)
which gives

$$p_a^m=K_a^m\left(\frac{a_0}{a}\right)^2,\dot{K}_a^m=0.$$
(46)
Using similar reasoning for $`\text{D}^2p_a^m`$ and with the help of the identity

$$\left(\text{D}^2f\right)^.=\text{D}^2\dot{f}2H\text{D}^2f+\dot{f}\text{D}^a\dot{u}_a,$$
(47)
from , we get

$$\text{D}^2p_a^m=M_a^m\left(\frac{a_0}{a}\right)^4,\dot{M}_a^m=0.$$
(48)
We now define the scalar parts of the perturbations by

$$\delta a\text{D}^a\delta _a,Za\text{D}^a\theta _a.$$
(49)
From (46) it follows that

$$a\text{D}^ap_a^m=K\left(\frac{a_0}{a}\right)^2,\dot{K}=0,$$
(50)
where $`K=a\text{D}^aK_a^m`$, and that

$$a\text{D}^a[a\text{D}_a(\text{D}^b\dot{u}_b)]=M\left(\frac{a_0}{a}\right)^4,\dot{M}=0,$$
(51)
where $`Ma\text{D}^aM_a^m`$. From equations (33) and (34) it follows that

$$a\text{D}^a(a\dot{u}_a)=K\left(\frac{a_0}{a}\right)^2.$$
(52)
Operating on equation (37) with $`a\text{D}^a`$ gives

$$\dot{Z}+2HZ+\frac{1}{2}\rho _m\delta =M\left(\frac{a_0}{a}\right)^4K\left(\frac{a_0}{a}\right)^2[4\rho _r+\frac{3}{2}\rho _m\frac{9}{2}(1+w)H^2].$$
(53)
For the matter scalar perturbations, $`w=0`$ and equation (2) implies $`Z=\dot{\delta }`$. The evolution equation for $`\delta `$ written in terms of the scaled independent variable

$$x=\alpha a=\frac{\rho _0^m}{\rho _0^r}\frac{a}{a_0}$$

is

$$\delta ^{^{\prime \prime }}+\frac{3x+2}{2x(x+1)}\delta ^{^{}}\frac{3}{2}\frac{1}{x(x+1)}\delta =F(x),$$
(54)
where the prime now denotes a derivative with respect to the new independent variable $`x`$ and

$$F(x)=\frac{3}{1+x}\left[\frac{M}{\rho _0^rx^2}+\frac{5}{2}K\left(\frac{\rho _0^m}{\rho _0^r}\right)^2\frac{1}{x^4}\right]$$
(55)
is the inhomogeneous part. It is immediately apparent that

$$\delta _1=3x+2$$
(56)
is a solution to the homogeneous part of equation (54). The method of reduction of order  then leads to a complete solution. The second solution of the homogeneous part is

$$\delta _2=\delta _1\frac{dx}{(3x+2)^2x\sqrt{x+1}},$$
(57)
which can be integrated using partial fractions to obtain

$$\delta _2=\frac{3}{2}\sqrt{x+1}+\frac{(3x+2)}{4}\mathrm{ln}\left[\frac{\sqrt{x+1}1}{\sqrt{x+1}+1}\right].$$
(58)
The homogeneous form of equation (54) was obtained by Mészáros  in 1974 who obtained a solution by transforming the equation to a hypergeometric form. Groth and Peebles  obtained the solutions (56 , 58) and gave the asymptotic behaviour

$`\begin{array}{cccccc}\delta _2\hfill & \hfill & \frac{1}{2}\mathrm{ln}(\tau /8),\hfill & x\hfill & \hfill & 1,\hfill \\ \delta _2\hfill & \hfill & 8/(45\tau ),\hfill & x\hfill & \hfill & 1,\hfill \end{array}`$

where

$$\tau =\frac{1}{\sqrt{3}}\frac{(\rho ^m)^2}{(\rho ^r)^{3/2}}t.$$

The results are further described in Peebles’s book .

The full solution to the inhomogeneous equation (54) is given by

$$\delta =C_1\delta _1+C_2\delta _2+\delta _1\frac{1}{E\delta _1^2}\left(E\delta _1F𝑑x\right)𝑑x,$$
(59)
where $`E=\mathrm{exp}\left((3x+2)[2x(x+1)]^1𝑑x\right)=x\sqrt{x+1}`$.

The expression $`E\delta _1F`$ can be written in the form

$$E\delta _1F=\frac{1}{\sqrt{x+1}}\left\{A_0+\frac{A_1}{x}+\frac{A_0}{3x^2}+\frac{A_1}{3x^3}\right\}$$
(60)
where $`A_0=9M/\rho _0^r`$, $`A_1=\frac{45}{2}K\left(\rho _0^m/\rho _0^r\right)^2`$. From this the particular integral can be calculated term by term. After integration, the first term $`\delta _{p1}`$ of the particular integral is given by

$`\delta _{p1}`$ $`=`$ $`\delta _1{\displaystyle \frac{1}{E\delta _1^2}\left(\frac{A_0}{\sqrt{x+1}}𝑑x\right)𝑑x}`$
$`=`$ $`A_0+{\displaystyle \frac{A_0}{2}}(3x+2)\mathrm{ln}\left({\displaystyle \frac{x}{3x+2}}\right).`$
The remaining terms are more difficult to integrate. From the first term it appears that the full inhomogeneous solution will contribute new modes. Among these it will contribute a constant mode and may contribute a mode similar to the second term in the solution $`\delta _2`$ and hence modify the coefficient of the factor $`(3x+2)`$ which could affect the asymptotic behaviour. The significance of this has to be viewed with caution because the remaining, still to be determined terms, may contribute further modes or cause cancellations. From the physics one is led to expect that the contribution of a small pressure perturbation will only be decaying modes.

Two other features of equation (54) are worthy of note. First, for large $`x`$ the asymptotic form of the homogeneous equation is the usual equation for first order perturbations in dust. Second, even for initial conditions in which $`\delta `$ and its first derivative with respect to $`x`$ are zero, perturbations will arise from the influence of the function $`F(x)`$. In other words small pressure effects or number density fluctuations can source density perturbations. This is possibly more significant, for the formation of structure, than the modification of the modes.

## IV Conclusion

For analytical as opposed to numerical modelling of the evolution of inhomogeneities in the universe, it is conventional to use a fluid approximation. At a detailed level this is at variance with reality because the matter is generally more particulate than hydrodynamics allows. Also most of the analytical literature assumes that the cosmological non-relativistic matter is dust. If we take the particulate nature into account, then at early times weak self-interactions or small random motions of collisionless matter give rise to a small but non-zero pressure in a fluid model. With this motivation we have discussed some of the implications of including a first order pressure in the matter distribution for fluid models.

We have derived the relation between the energy density perturbations and the number density perturbations for the isentropic case and we show that in the simple non-isentropic case with $`p/\rho =\mathrm{constant}`$ a new stationary mode appears. This new mode will affect the formation of structure at certain stages in the evolution of the inhomogeneity.

In the main body of the paper we derive an exact solution to the energy density perturbation equation when the matter pressure is non-zero at first order. The implications of assuming a non-zero first order matter pressure are quite deep as is illustrated in equation (25), which is an equation for second order gauge invariant quantities, and in equation (33), which together with (32) shows that if $`\text{D}_ap_m`$ were equal to $`0`$, then $`\dot{u}_a=0=\text{D}_ap_r`$. This would reduce the problem to the standard case.

The complete solution (59) is different from the usual solutions used in CDM approximations and from the Mészáros solution. Having such a solution may be useful as an analytical tool for understanding some features of recent numerical work on weakly self-interacting or warm dark matter models in the non-linear regime . This is a subject of further investigation.

Acknowledgements

The authors thank Roy Maartens and Marco Bruni for useful discussions on this work. J. Triginer was supported in part by a Royal Society Joint Project Grant.


text:
# Abstract

## Abstract

Within the so called $`𝒫𝒯`$ symmetric version of quantum mechanics a brief review of the exactly solvable models is given. Distinction is made between the curved and straight coordinate lines, between their unbounded (aperiodic) and bounded (periodic) choices, and between the completely and partially solvable cases.

## \[KEYWORDS\]

\[ $`𝒫𝒯`$ symmetry, Schrödinger equation in complex domain, exact solutions\]

## \[AMS 1991 Mathematics Subject Classification\]

81P10 81Q05 81R40 47B50 34A20 15A57

## 1 Introduction

$`𝒫𝒯`$ symmetric quantum mechanics  has independently been proposed and used as a methodical laboratory in quantum physics by several groups of authors: By Caliceti et al  in perturbation theory, by Bessis et al  in field theory and by Andrianov et al  in supersymmetric context.

This short review will pay attention to the (partially as well as completely) exactly solvable models within this framework, with emphasis on the results obtained by the present author.

## 2 Complete solvability on curved paths

One of the first exactly and completely solvable examples of a $`𝒫𝒯`$ symmetric system has been found by Cannata et al  and re-discovered by Bender et al  more than one year later. Its modified Schrödinger bound state problem is defined on certain curved, left-right symmetric “generalized coordinate” lines in the complex plane. Mathematically, it is defined via exponential potential and proves exactly solvable in terms of Bessel functions. Due to its relationship to a power-law forces in the large exponent limit, it can be most simply interpreted as a certain smooth and non-Hermitian $`𝒫𝒯`$ symmetric analogue of the current square well.

Recently, a double-well counterpart of the latter set of models has been shown exactly solvable, in terms of Laguerre polynomials, in ref. . In a way similar to the above “single well” example its paths of integration are the same, curved complex lines again. Their spectra exhibit a puzzling and highly unexpected feature of certain coupling-dependent re-arrangements mediated by “unavoided” crossings at critical points. The phenomenon reflects the non-Hermiticity of the Hamiltonian .

Via a suitable Liouvillean change of the variables in the above double-well-like differential Schrödinger equation one can immediately obtain another Laguerre-related solvable system with the potential of the Coulombic single-pole form. In this case the deformation of the integration path plays the beneficial role of a natural regularization prescription. At the same time it also leads to the necessity of working with the complex charges in a way described in ref. . Also the related spectrum of energies exhibits certain unexpected features: Positivity, a coexistence of the growth and decrease with the increasing coupling strength, etc.

A clear mathematical explanation of behaviour as well as a consistent and/or possible physical interpretation of the counterintuitive models of this type remain unclarified up to now. Much better situation emerges in the case of the integration curves defined as the left-right-symmetric straight lines.

## 3 Straight paths

Using the language of the so called Kustaanheimo - Steifel transformation (cf.  for detailed references) both the above curvilinear examples can be shown equivalent to the $`𝒫𝒯`$ symmetric harmonic oscillator of ref.  with a centrifugal term regularized by the mere downward complex shift of the (full) real axis,

$$\left(\frac{d^2}{dx^2}+x^22icx+\frac{\alpha ^21/4}{(xic)^2}\right)\phi (x)=(E+c^2)\phi (x),\phi (x)L_2(\mathrm{},\mathrm{}).$$
(1)
This oscillator with the Laguerre-polynomial normalizable solutions

$$\phi _{(\pm n)}(x)=𝒩(xic)^{\pm \alpha +1/2}e^{(xic)^2/2}L_n^{(\pm \alpha )}\left[(xic)^2\right],n=0,1,\mathrm{}$$

possesses the non-equidistant spectrum of energies $`E=E_{(\pm n)}=4n+2\pm 2\alpha `$ and represents a certain unperturbed limit of the quartic oscillator models of Buslaev and Grecchi . These authors, unfortunately, did not notice the existence of the “quasi-even”, $`{}_{(n)}{}^{}`$signed half of the spectrum. This omission can be easily corrected. One just introduces a “two-to-one” isospectrality correspondence between the respective Hermitian and non-Hermitian anharmonic oscillator models of ref. .

Interpretation of the models “living” on the straight lines becomes significantly facilitated by the easier identification and interpretation of their complex components . Immediate purely analytic constructions recover, e.g., the existence of models which are in a one-to-one correspondence to the so called shape-invariant real forces in one dimension (cf. their presentation in ref. ) and on the half line (their $`𝒫𝒯`$ symmetric counterparts were described and listed in ref. ).

The situation is reviewed in ref. . A fully general form of this type of analytic constructions dates back to the introduction of the so called Natanzon potentials and, in the present context, is thoroughly analyzed and described in ref. .

## 4 Models with periodic boundary conditions

A new and promising development of the $`𝒫𝒯`$ symmetric considerations has been recently inspired by the study of the two- and three-particle models . The $`𝒫𝒯`$ symmetrization of the Hamiltonians has been again conjectured to be sufficient for keeping their spectrum real. The related “weakening of the Hermiticity” finds a natural generalization in the new context.

Particular attention has been paid to the possible non-Hermitian generalizations of the well known Calogero model . In this setting the separability of the underlying partial differential Schrödinger equation in the hyperspherical coordinates helps us to reduce the problem to the “hyperangular” ordinary differential equation defined on a finite interval. In this way, in the simplest cases one has to solve the complexified ordinary differential equations of the generalized Pöschl-Teller type,

$$\left(\frac{d^2}{d\varphi ^2}+\frac{\mathrm{}(\mathrm{}+1)}{\mathrm{sin}^2\varphi }+\frac{\lambda (\lambda +1)}{\mathrm{cos}^2\varphi }\right)\chi (\varphi )=E\chi (\varphi )$$
(2)
on an interval $`\varphi (M\pi /2,M\pi /2)`$ with a suitable integer $`M`$. These equations can be solved exactly in terms of the hypergeometric functions .

The strongly repulsive singularities at $`\varphi _j=j\pi /2`$ are currently not penetrable . Here, they become regularized in a $`𝒫𝒯`$ symmetric manner which parallels a few older constructions on the unbounded intervals . The quasi-symmetric and quasi-antisymmetric solutions arise from certain ad hoc boundary conditions .

For the most elementary illustration let us now choose $`\lambda =0`$ and $`M=2`$. Then, the differential equation (2) possesses the two independent hypergeometric solutions,

$$\chi ^{(\pm )}(\varphi )=(\mathrm{sin}\varphi )_2^{1/2\pm \alpha }F_1(u^{(\pm )},v^{(\pm )};1\pm \alpha ;\mathrm{sin}^2\varphi ),\alpha =\mathrm{}+1/2>0$$

where $`2u^{(\pm )}=1/2\beta \pm \alpha `$ and $`2v^{(\pm )}=1/2+\beta \pm \alpha `$. On the boundary of convergence $`\mathrm{sin}^2\varphi =1`$ the matching of the logarithmic derivatives is equivalent to the termination of this series to the Gegenbauer polynomials,

$$\chi (\varphi )=\chi _{(\pm k)}(\varphi )=(\mathrm{sin}\varphi )^{1/2\pm \alpha }C_k^{1/2\pm \alpha }(\mathrm{cos}\varphi ),k=0,1,\mathrm{}.$$

The construction also quantizes the energies and gives them in the closed form,

$$E=E_{(\pm k)}=(k\pm \alpha +1/2)^2,k=0,1,\mathrm{}.$$

This set of eigenvalues is composed of the two subsets in a way which resembles the above-mentioned non-equidistant spectrum of the $`𝒫𝒯`$ symmetrized singular harmonic oscillator (1.

## 5 Conclusions and outlook

Historically, one of the first persuasive manifestations of the merits and power of the $`𝒫𝒯`$ symmetry has been offered by Bender and Boettcher  who discovered the quasi-exact (i.e., incomplete) solvability of the most common and popular quartic polynomial oscillators.

The plausible reasons of the unexpected delay of such an “obvious” observation are closely related to the above-mentioned “forgotten” energies. One has to keep in mind the “spontaneous” regularity of the singular potentials within the new formalism. This observation was made explicit in our paper  where the presence of the two additional singular terms has still been shown compatible with the quasi-exact solvability of quartic potentials.

In the quasi-exact context the changes of variables can play the same role as in the completely solvable models. This has been illustrated by the particular constructions of the decadic model  and of the harmonic + Coulomb superposition . Further work in this direction is in progress .

## Acknowledgements

Research assisted by the Grant Agency of the Academy of Sciences of the Czech Republic, contract Nr. A 1048004.