text:
# Abstract

## Abstract

Schrödinger equation with imaginary $`𝒫𝒯`$ symmetric potential $`V(x)=ix^3`$ is studied using the numerical discretization methods in both the coordinate and momentum representations. In the former case our results confirm that the model generates an infinite number of bound states with real energies. In the latter case the differential equation is of the third order and a square-well, solvable approximation of kinetic energy is recommended and discussed. One finds that in the strong-coupling limit, the exact $`𝒫𝒯`$ symmetric solutions converge to their Hermitian predecessors.

PACS 03.65.Ge

### Acknowledgements

Work supported by the GAČR grant Nr. 202/07/1307, by the MŠMT “Doppler Institute” project Nr. LC06002 and by the Institutional Research Plan AV0Z10480505.

## 1 Introduction

An interest in imaginary cubic anharmonic oscillators dates back to their perturbation analysis by Caliceti et al . The simplified homework example with the mere two-term non-Hermitian Hamiltonian

$$H_{BZ}=p^2+ix^3$$

has been proposed by D. Bessis and J. Zinn-Justin who had in mind its possible applicability in the context of statistical physics . The example has been re-vitalized by C. Bender et al due to its possible methodical relevance in relativistic quantum field theory . They emphasized that the apparent reality of the spectrum of energies $`E_{BZ}`$ is quite puzzling.

The conjecture of a full compatibility of similar Hamiltonians with the postulates of quantum mechanics  opened many new and interesting questions. The current Hermiticity of Hamiltonians was replaced by a weaker condition of their commutativity with a product $`𝒫𝒯`$ of the spatial parity $`𝒫`$ and of the complex conjugation $`𝒯`$. The latter factor is to be understood as a one-dimensional version of the operator of time reversal.

The discussion between C. Bender and A. Mezincescu  pointed out that one of the key problems of the new studies lies in the ambiguity of the spectrum which depends quite crucially on our choice of the boundary conditions which can be, in general, complexified . The fragile character of the reality of the energies has been confirmed by the WKB and perturbative studies  and by the quasi-exact and exact models  where the admissible unavoided level crossings  prove sometimes followed by the spontaneous breakdown of $`𝒫𝒯`$ symmetry .

In such a context we propose here an extremely elementary approach to similar models replacing interactions which admit just a numerical treatment (typically, $`V(x)=ix^3`$) by their exactly solvable square-well analogues.

## 2 Models in coordinate representation

In a search for analogies between the solvable and unsolvable models in one dimension, all of the possible forms of a confining well are often being approximated by the ordinary real and symmetric square well

$$V^{(SQW)}(x)=\{\begin{array}{cc}S^2,\hfill & x(\mathrm{},\pi )(\pi ,\mathrm{}),\hfill \\ 0,\hfill & x(\pi ,\pi ).\hfill \end{array}$$
(1)
In this spirit one can also replace the antisymmetric and imaginary homework potential $`V_{BZ}(x)=ix^3`$ by its elementary square-well analogue

$$V^{(ISQW)}(x)=\{\begin{array}{cc}iT^2,\hfill & x(\mathrm{},\pi ),\hfill \\ 0,\hfill & x(\pi ,\pi ),\hfill \\ +iT^2,\hfill & x(\pi ,\mathrm{}).\hfill \end{array}$$

Schrödinger equation which appears in such a setting,

$$\left[\frac{\mathrm{}}{2m}\frac{d}{dx^2}+V^{(ISQW)}(x)\right]\psi (x)=E\psi (x)$$
(2)
will be complemented by the standard $`\mathrm{L}^2(lR)`$ boundary conditions

$$\psi (\pm \mathrm{})=0.$$
(3)
The well known $`𝒫𝒯`$ symmetric normalization convention will be employed, with a free real parameter $`G`$ in the unbroken $`𝒫𝒯`$-symmetry requirements

$$\psi (0)=1,_x\psi (0)=iG.$$
(4)
Putting $`\mathrm{}=2m=1`$ and using the ansatz

$$\psi (x)=\{\begin{array}{ccc}\mathrm{cos}kx+B\mathrm{sin}kx,\hfill & x(0,\pi ),\hfill & k^2=E,\\ (L+iN)\mathrm{exp}(\sigma x),\hfill & x(\pi ,\mathrm{}),\hfill & \sigma ^2=iT^2k^2,\end{array}$$
(5)
we shall guarantee the full compatibility of such a convention with the symmetry requirements (4) by the choice of the purely imaginary constant $`B=iG/k`$ in wave functions (5).

## 3 Matching conditions at $`x=\pi `$

Let us split $`\sigma =p+iq`$ in its real and imaginary part with a fixed sign, $`p,q0`$. This gives $`p^2+k^2=q^2`$ and $`2pq=T^2`$. These rules are easily re-parameterized in terms of a single variable $`\alpha `$,

$$p=q\mathrm{cos}\alpha ,k=q\mathrm{sin}\alpha ,q=\frac{T}{\sqrt{2\mathrm{cos}\alpha }},\alpha (0,\pi /2).$$
(6)
The standard matching at the point of discontinuity is immediate,

$$\mathrm{cos}k\pi +B\mathrm{sin}k\pi =(L+iN)\mathrm{exp}(\sigma \pi ),$$

$$\mathrm{sin}k\pi +B\mathrm{cos}k\pi =\frac{\sigma }{k}(L+iN)\mathrm{exp}(\sigma \pi ).$$

After we abbreviate $`\sigma /k=\mathrm{tan}\mathrm{\Omega }\pi `$, we get an elementary complex condition of matching of logarithmic derivatives at $`x=\pi `$,

$$G=ik\mathrm{tan}(k+\mathrm{\Omega })\pi .$$
(7)
The real part defines our first unknown parameter, $`G=G(\alpha )`$. Due to our normalization conventions, the imaginary part of the right-hand-side expression must vanish, $`\mathrm{Re}[\mathrm{tan}(k+\mathrm{\Omega })\pi ]=0`$. An elementary re-arrangement of such an equation acquires the form of an elementary quadratic algebraic equation for $`X=\mathrm{tan}k\pi `$. Its two explicit solutions read

$$X_1=\frac{p+q}{k},X_2=\frac{pq}{k}$$
(8)
or, after all the insertions,

$$\mathrm{tan}\left[\frac{\pi T\mathrm{sin}\alpha ^{(+)}}{\sqrt{2\mathrm{cos}\alpha ^{(+)}}}\right]=\mathrm{tan}\left[\frac{\pi \alpha ^{(+)}}{2}\right],$$
(9)
$$\mathrm{tan}\left[\frac{\pi T\mathrm{sin}\alpha ^{()}}{\sqrt{2\mathrm{cos}\alpha ^{()}}}\right]=\mathrm{tan}\left[\frac{\alpha ^{()}}{2}\right].$$
(10)
In implicit manner these equations specify the two respective infinite series of appropriately bounded real roots $`\alpha =\alpha _n^{(\pm )}(0,\pi /2)`$.

## 4 Energies

For $`\alpha (0,\pi /2)`$ the left-hand-side arguments $`[\mathrm{}]`$ in eqs. (9) and (10) run from zero to infinity. Their tangens functions oscillate infinitely many times from minus infinity to plus infinity. Within the same interval, the limited variation of the argument $`\alpha `$ makes both the eligible right-hand side functions monotonic, very smooth and bounded, $`\mathrm{tan}[(\pi \alpha ^{(+)})/2](1,\mathrm{})`$ and $`\mathrm{tan}[\alpha ^{()}/2](0,1)`$. A priori this indicates that our roots $`k=k(\alpha _n^{(\pm )})`$ will all lie within well determined intervals,

$$k_n^{(+)}(n+\frac{1}{4},n+\frac{1}{2}),n=0,1,\mathrm{},$$

$$k_m^{()}(m+\frac{3}{4},m+1)m=0,1,\mathrm{}.$$

After such an approximate localization of the roots, an unexpected additional merit of our parametrization (6) manifests itself in an unambiguous removal of the tangens operators from both eqs. (9) and (10). This gives the following two relations,

$$k_n^{(+)}=n+\frac{1}{2}\frac{\omega _n^{(+)}}{4},k_m^{()}=m+1\frac{\omega _m^{()}}{4},\omega _n^{(\pm )}=\frac{2\alpha _n^{(\pm )}}{\pi }(0,1).$$

After an elementary change of notation with $`\omega _n^{(+)}=\omega _{2n}`$ and $`\omega _n^{()}=\omega _{2n+1}`$, we may finally combine the latter two rules in the single secular equation

$$\mathrm{sin}\left(\frac{\pi }{2}\omega _N\right)=\frac{2N+2\omega _N}{4T}\sqrt{2\mathrm{cos}\left(\frac{\pi }{2}\omega _N\right)}N=0,1,\mathrm{},$$
(11)
In a graphical interpretation this equation represents an intersection of a tangens-like curve with the infinite family of parallel lines. This is illustrated in Figure 1. The equation generates, therefore, an infinite number of the real roots $`\omega _N(0,1)`$ at all the non-negative integers $`N=0,1,\mathrm{}`$.

## 5 Wave functions in the weak coupling regime

Equation (7) in combination with eqs. (9) and (10) determines the real parameter

$$G=G^{(\pm )}=\frac{k^2}{q\pm p}$$
(12)
responsible for the behavior of the wave functions near the origin \[remember that $`B=iG/k`$ in eq. (5)\]. For its analysis let us introduce an auxiliary linear function of $`\omega `$ and $`N`$,

$$\sqrt{R(\omega _N,N)}=\frac{2N+2\omega _N}{4T}(\frac{N+1/2}{2T},\frac{N+1}{2T}).$$

Our secular eq. (11) can be then read as an algebraic quadratic equation with the unique positive solution,

$$\mathrm{cos}\left(\frac{\pi }{2}\omega _N\right)=\frac{1}{R(\omega _N,N)+\sqrt{R^2(\omega _N,N)+1}}.$$
(13)
This is an amended implicit definition of sequence $`\omega _N`$. As long as the right hand side expression is very smooth and never exceeds one, the latter formula re-verifies that the root $`\omega _N`$ is always real and bounded as required.

In the domain of the large and almost constant $`R1`$ (i.e., for small square-well heights $`T`$ or at the higher excitations), our new secular equation (13) gives a better picture of our bound-state parameters $`\omega _N=1\eta _N`$ which all lie very close to one. The estimate

$$\frac{\pi }{2}\eta _N=\mathrm{arcsin}\frac{1}{R+\sqrt{R^2+1}}\frac{1}{2R}\frac{5}{48R^3}+\mathrm{}$$

represents a quickly convergent iterative algorithm for the efficient numerical evaluation of the roots $`\omega _N`$. One can conclude that in a way compatible with our a priori expectations, the value of $`p=p_N=\mathrm{Re}\sigma q/2R`$ lies very close to zero. As a consequence, the asymptotic decrease of our wave functions remains slow. We have $`q=q_N=\mathrm{Im}\sigma k`$ so that, asymptotically, our wave functions very much resemble free waves $`\mathrm{exp}(ikx)`$. In the light of eq. (12) we have $`\psi (x)\mathrm{exp}(ikx)`$ near the origin.

## 6 Wave functions in the strong coupling regime

For the models with a very small $`R`$ (i.e., for the low-lying excitations in a deep well with $`T1`$) we get an alternative estimate

$$\frac{\pi }{4}\omega _N=\mathrm{arcsin}\sqrt{\frac{1}{2}\left[R\left(\sqrt{1+R^2}1\right)\right]}\frac{1}{2}R\frac{1}{4}R^2+\mathrm{}\frac{\pi }{4}.$$

In the limit $`R0`$ the present spectrum of energies moves towards (and precisely coincides with) the well known levels of the infinitely deep Hermitian square well of the same width $`I=(\pi ,\pi )`$ (cf. eq. (1) with $`S\mathrm{}`$). The complex-rotation transition from the Hermitian well $`V^{(SQW)}(x)`$ of eq. (1) (with $`S1`$) to its present non-Hermitian $`𝒫𝒯`$ symmetric alternative $`V^{(ISQW)}(x)`$ of eq. (2) (with $`T1`$) proves amazingly smooth.

Wave functions exhibit a similar tendency. In outer region, they are proportional to $`\mathrm{exp}(px)`$ and decay very quickly since $`p=𝒪(R^{1/2})`$. Parameter $`G^{(\pm )}`$ becomes strongly superscript-dependent,

$$G^{(+)}=\frac{k^2}{q+p}=𝒪(R^{3/2}),G^{()}=(q+p)=𝒪(R^{1/2}).$$

In the interior domain of $`x(\pi ,\pi )`$ the wave functions with superscripts (+) and (-) become dominated by their spatially even and odd components $`\mathrm{cos}kx`$ and $`\mathrm{sin}kx`$, respectively. The superscript mimics (or at least keeps the trace of) the quantum number of the slightly broken spatial parity $`𝒫`$.

We can summarize that our present $`𝒫𝒯`$ symmetric model is quite robust. Independently of the coupling $`T`$ the spectrum is unbounded from above and remains constrained by inequalities

$$\frac{(N+1/2)^2}{4}E_N\frac{(N+1)^2}{4}.$$
(14)
The analogy between our exactly solvable square-well model and the standard or “paradigmatic” $`𝒫𝒯`$ symmetric Hamiltonian $`H_{BZ}`$ appears closer than expected.

## 7 Transition to the momentum representation

Let us turn our attention to one-dimensional harmonic oscillator $`H^{(HO)}=p^2+x^2`$ which is exactly solvable and which appears in virtually any textbook on quantum mechanics. In $`𝒫𝒯`$ symmetric quantum mechanics a similar guiding role can be and has been attributed to the non-Hermitian cubic Hamiltonian $`H^{(CO)}=p^2+ix^3`$ of Bender et al . We have seen that a straightforward numerical and semi-classical analysis of the related Schrödinger equation

$$H^{(CO)}|\psi _n=E_n^{(CO)}|\psi _n,n=0,1,\mathrm{}$$
(15)
supports a highly plausible conjecture that the spectrum of energies is real, discrete and bounded below. The conjectured absence of its imaginary components is indicated by the Hilbert-Schmidt analysis  and by the perturbation calculations in both the weak-coupling regime  or in its strong-coupling, purely numerically generated re-arrangement .

Certainly, the problem deserves a change of the traditional perspective. Let us, therefore, move now to its momentum representation. This would give the momentum as a mere number, $`p(\mathrm{},\mathrm{})`$ while the coordinate $`x`$ becomes represented by the differential operator $`\widehat{x}=i_p`$. Equation (15) then acquires a new form containing the purely real differential Schrödinger operator $`H^{(CO)}E`$ of third order,

$$\left[\frac{d^3}{dp^3}+p^2\right]\psi (p)=E\psi (p),\psi (p)L_2(IR).$$

This gives an unusual formulation of our bound-state problem where the quadratic $`p`$dependence of the kinetic term $`T(p)=p^2`$ does not seem to make the equation any easier to solve. For this reason we shall drastically simplify the kinetic term and deduce some consequences.

## 8 Piecewise constant approximate kinetic energy

In a way proposed by Prüfer  many wave functions can be visualized as certain deformations of solutions which correspond to a locally constant potential, $`\psi (x)c_1\mathrm{sin}[\varrho (x)]+c_2\mathrm{cos}[\varrho (x)]`$. In the standard quantum mechanics such a trick found immediate applications in numerical computations  while it still admits an easy interpretation via some traditional Sturm Liouvillean oscillation theorems . Using this idea as a methodical guide let us now replace the kinetic energy operator $`T(p)=p^2`$ by the most elementary square well of a finite depth $`Z>0`$,

$$T(p)=\{\begin{array}{cc}Z,\hfill & p(\mathrm{},1),\hfill \\ 0,\hfill & p(1,1),\hfill \\ Z,\hfill & p(1,\mathrm{}).\hfill \end{array}$$

In the bounded range of energies $`E<Z`$ this splits our toy model in the two separate differential equations,

$$\left[\frac{d^3}{dp^3}8\alpha ^3\right]\psi (p)=0,p(1,1),$$

$$\left[\frac{d^3}{dp^3}+8\beta ^3\right]\psi (p)=0,p(\mathrm{},1)(1,\mathrm{}).$$

The two auxiliary parameters $`\alpha =\alpha (E)>0`$ and $`\beta =\beta (E)>0`$ are defined in such a way that $`Z=E+8\beta ^3>E=8\alpha ^3`$. They appear in the three independent (exponential) solutions of our equation. Their general superpositions are complex but they may be given the real, trigonometric form. Near the origin we have

$$\psi _0(p)=de^{2\alpha p}+fe^{\alpha p}\mathrm{cos}(\stackrel{~}{\alpha }p+\theta ),p(1,1),\stackrel{~}{\alpha }=\sqrt{3}\alpha $$

where the symbols $`d`$, $`f`$ and $`\theta `$ stand for the three undetermined constant parameters. In the right and left asymptotic regions we obtain the similar formulae. After we omit their exponentially growing and normalization-violating unphysical components we get the one-parametric family

$$\psi _+(p)=ge^{2\beta p},p(1,\mathrm{}).$$

The two-parametric left-barrier counterpart of this formula reads

$$\psi _{}(p)=ce^{\beta p}\mathrm{cos}(\stackrel{~}{\beta }p+\eta ),p(\mathrm{},1),\stackrel{~}{\beta }=\sqrt{3}\beta .$$

At the right discontinuity $`p=1`$ we have to guarantee the continuity of $`\psi (p)`$, $`_p\psi (p)`$ and $`_p^2\psi (p)`$. This is equivalent to the three matching conditions

$$\begin{array}{c}de^{2\alpha }+fe^\alpha \mathrm{cos}(\stackrel{~}{\alpha }+\theta )=ge^{2\beta },\\ 2\alpha de^{2\alpha }\alpha fe^\alpha [\mathrm{cos}(\stackrel{~}{\alpha }+\theta )+\sqrt{3}\mathrm{sin}(\stackrel{~}{\alpha }+\theta )]=2\beta ge^{2\beta },\\ 4\alpha ^2de^{2\alpha }2\alpha ^2fe^\alpha [\mathrm{cos}(\stackrel{~}{\alpha }+\theta )\sqrt{3}\mathrm{sin}(\stackrel{~}{\alpha }+\theta )]=4\beta ^2ge^{2\beta }.\end{array}$$
(16)
The weighted sum of these equations re-scales and interrelates the unknown coefficients $`d`$ and $`g`$ in terms of a new energy parameter $`t=t(E)=\beta (E)/\alpha (E)>0`$ or rather $`R=R(E)=(1t+t^2)^{1/2}>0`$,

$$de^{2\alpha }D(E)=\frac{G(E)}{3R^2(E)},G(E)ge^{2\beta }.$$

After we eliminate $`G`$ from the last two equations (16) which are linear in $`d`$ and $`f`$ we obtain the elementary formula which defines the shift $`\theta =\theta (E)`$,

$$\mathrm{tan}(\stackrel{~}{\alpha }+\theta )=\frac{\sqrt{3}}{2/t1}.$$

The trigonometric factors become fixed up to their common sign $`\epsilon =\pm 1`$,

$$\mathrm{cos}(\stackrel{~}{\alpha }+\theta )=\epsilon [1t(E)/2]R(E),\mathrm{sin}(\stackrel{~}{\alpha }+\theta )=\frac{\sqrt{3}}{2}\epsilon t(E)R(E).$$

The same sign enters our last decoupled definition

$$fe^\alpha F(E)=\frac{2[t(E)+1]G(E)}{3\epsilon R(E)}.$$

Up to an overall normalization (say, $`g=1`$) in our wave functions matched at $`p=1`$, all the free parameters become specified as functions of the energy.

At the second, $`p=1`$ discontinuity we have to satisfy three matching conditions as well. In terms of abbreviations

$$\begin{array}{c}L_1=de^{2\alpha }+fe^\alpha \mathrm{cos}(\stackrel{~}{\alpha }+\theta ),\\ L_2=2\alpha de^{2\alpha }\alpha fe^\alpha [\mathrm{cos}(\stackrel{~}{\alpha }+\theta )+\sqrt{3}\mathrm{sin}(\stackrel{~}{\alpha }+\theta )],\\ L_3=4\alpha ^2de^{2\alpha }2\alpha ^2fe^\alpha [\mathrm{cos}(\stackrel{~}{\alpha }+\theta )\sqrt{3}\mathrm{sin}(\stackrel{~}{\alpha }+\theta )]\end{array}$$

they read

$$\begin{array}{c}L_1(\alpha ,d,f,\theta )=ce^\beta \mathrm{cos}(\stackrel{~}{\beta }+\eta ),\\ L_2(\alpha ,d,f,\theta )=\beta ce^\beta [\mathrm{cos}(\stackrel{~}{\beta }+\eta )\sqrt{3}\mathrm{sin}(\stackrel{~}{\beta }+\eta ),\\ L_3(\alpha ,d,f,\theta )=2\beta ^2ce^\beta [\mathrm{cos}(\stackrel{~}{\beta }+\eta )+\sqrt{3}\mathrm{sin}(\stackrel{~}{\beta }+\eta ).\end{array}$$

They determine the values of $`c=c(E)`$ and $`\eta =\eta (E)`$. Their properly weighted sum gives

$$\sqrt{3}F(E)\mathrm{sin}(\stackrel{~}{\alpha }+\theta )+(2t1)F(E)\mathrm{cos}(\stackrel{~}{\alpha }+\theta )+2\frac{1t+t^2}{t+1}D(E)e^{6\alpha }=0.$$

After its slight re-arrangement we arrive at the amazingly transparent relation

$$(14t+t^2)\mathrm{cos}(2\sqrt{3}\alpha )+\sqrt{3}(1t^2)\mathrm{sin}(2\sqrt{3}\alpha )=\left(\frac{1t+t^2}{1+t}e^{3\alpha }\right)^2.$$
(17)
It should be read as an implicit definition of physical energies $`E`$.

## 9 The $`Z`$dependence of the spectrum

It is quite instructive to search for the physical energies numerically. Starting from the very-shallow-well extreme in eq. (17) we find the two clearly distinguished energy roots. The qualitative features of the graph of secular determinant remain unchanged in a broad interval of the inverse strengths $`1/Z`$. Its shape is sampled in Figure 2 at $`Z=10^3`$. We have tested that even the approximate height $`5`$ of its left plateau stays virtually unchanged between $`Z=10^5`$ and $`Z=10^3`$. Within the same interval of the shallowest wells the left zero grows from the value $`0.00280`$ till $`0.0241`$. Beyond the broad, downwards-oriented peak one finds the second, right zero moving from the value $`0.01047`$ (found very close to the instantaneous threshold $`0.01077`$) up the value $`0.1056`$ (not very far from its threshold $`0.1077`$, either) within the same interval of $`1/Z`$.

A new feature emerges around $`Z=10^1`$ (with the left zero at $`0.0445`$ and with the right zero $`0.222`$ still quite close to the threshold $`0.232`$) and $`Z=1`$ (with the left zero at $`0.072`$ and with the right zero $`0.446`$, not that close already to the threshold $`1/2`$) in the left half of the picture. The plateau develops a local, safely negative maximum.

In the subsequent domain of $`Z>1`$ we have to switch our attention back to the right half of our graph. Immediately before the coupling reaches the integer value of $`Z=5`$, the end of the curve returns to the negative half-plane near the maximal (i.e., threshold) energy. This means that there emerges the third energy level there. The total number of bound states grows to $`N=3`$ (cf. the leftmost items in our Table 1).

### Table 1.
Number of levels $`N`$ and its changes $`\mathrm{\Delta }`$ with growing $`Z`$.

$$\begin{array}{cccccccccccccccccccccccc}& & & & & & & & & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & & & & & & & & & \\ N& 2& & 3& & 5& & 3& & 4& & 5& & 7& & 8& & 6& & 7& & 9& & 10\\ & & & & & & & & & & & & & & & & & & & & & & & \\ \mathrm{\Delta }& & 1& & 2& & 2& & 1& & 1& & 2& & 1& & 2& & 1& & 2& & 1& \\ & & & & & & & & & & & & & & & & & & & & & & & \end{array}$$

###

Beyond $`Z=5`$, our attention has to return quickly to the left half of the picture where the very slow growth of the local maximum creates a new quality at last. The top of the local bump touches and crosses the horizontal axis at $`Z5.3003`$ and $`E0.6244`$. At $`Z=5.3005`$ a new doublet of energies is formed in a way illustrated in Figure 3. The number of levels jumps to $`N=5`$.

A smooth deformation of the graph takes place when the value of $`Z`$ grows on. During this evolution we discover that our (originally broad), downward-oriented peak shrinks quite quickly and moves comparatively slowly to the right. It gets close to the rightmost and, to its bad luck, slightly more slowly moving zero number five. The magnified picture of the resulting “collision” is displayed here in Figure 4. At $`Z=35`$ it shows that in the threshold region of the energies,

* the wavy motion of the threshold end of our graph still did not manage to reach the zero axis;
* the downwards-oriented peak has already left the positive part (and moved to the negative part) of the curve in question.

As a consequence, the number of levels drops, quite unexpectedly, down to $`3`$ again (cf. Table 1).

In the vicinity of $`Z=40`$ the new, rightmost energy root emerges at last. Up to $`Z=100`$ and beyond, the number of levels stays equal to $`4`$. Then it increases to $`5`$, due to the emergence of the next threshold zero. Only after that, the slowly moving downward peak reaches the domain of the fourth zero. At almost exactly $`Z=190`$ its left (and temporarily negative) local maximum reaches the zero value again (cf. Figure 5). At this moment the number of states jumps up by two to seven. The magnified graphical proof is offered by Figure 5 at $`Z=200`$.

The latter Figure illustrates nicely the rapid shrinking of our peak with $`Z`$. The numerical detection of its position becomes more and more difficult. Although this position plays a crucial role in the practical determination of the number of levels at a given $`Z`$, we must be very careful in distinguishing the subgraph of Figure 6 (with three zeros) from a simple straight line with the single zero.

The pattern is deceitful and the standard software which searches for roots has to be used with due care. Vice versa, the above analysis enables us to take into account all the specific features of the $`Z`$dependence of the graph in eq. (17). We get a regular pattern summarized in Table 1 and exhibiting a certain regularity of the $`Z`$dependence of the number of levels $`N=N(Z)`$.

## 10 Energies in the square-well approximation

The main consequence of the presence of the above-mentioned narrow peak is an unusual irregularity observed in the emergence of the new levels in the deeper wells. We may conclude that this irregularity is not an artifact of the computation method. The energy formula (17) for our square-well toy model is exact and the seemingly unpredictable emergence of its roots just reflects the fact that our Schrödinger equation is of the third order. In particular, there exists no symmetry/antisymmetry with respect to the parity $`pp`$ etc.

Methodical consequences of our analysis are a bit discouraging. Firstly, the very symbolic-manipulation derivation of our present formulae proved unexpectedly complicated even in comparison with the multiple standard square wells in textbooks. That’s why we did not move to any further piece-wise constant approximations of $`T(p)`$.

Secondly, even our use of the most elementary solvable example revealed quite clearly a very real danger of the possible loss of certain levels. For an illustration let us imagine that our numerical study would have been started in the deep-well domain, i.e., at the large $`Z`$. It is quite easy to generate the graphs of eq. (17) there. In the standard and routine finite-precision computer arithmetics one discovers that the results are very smooth and look virtually the same, say, in the interval of $`Z(1000,1200)`$.

Let us pick up, for definiteness, the larger sample $`Z=1200`$. We get a picture (cf. our last Figure 7) which is regular and, deceptively, indicates that $`N(1200)=7`$. Unfortunately, the correct answer (appearing at the right end of our Table 1) is $`N(1200)=10`$. Its derivation requires the use of a significantly enhanced precision. Otherwise, whenever we use just the standard 14 digits and Figure 7, we would have missed as much as three (i.e., cca 30 % of all) energy levels.

In the light of our preceding considerations, an easy explanation of the latter numerical paradox lies in the presence of the narrow peak. A priori, it is hardly predictable of course. It is necessary to spot it by brute force. One finds that at $`Z=1000`$, this anomalous peak still lives safely below the sixth energy level. The related number of levels is reliably confirmed as equal to seven, indeed.

In between $`Z=1000`$ and $`Z=1200`$, it is necessary to work in an enhanced precision arithmetics. One finds that the upper, threshold end of the curve crosses the horizontal axis only slightly above $`Z=1100`$. Due to the very steep slope of the curve in this region, this crossing is not visible even at $`z=1200`$ in Figure 7.

One has to trace the narrow peak carefully. It overtakes the sixth energy level at $`Z1190`$ (and $`E=4.217`$), in an arrangement resembling our Figure 5 above. Thus, one concludes, finally, that the new, almost degenerate pair of the energy levels emerges immediately beyond this point.

## 11 Wave functions and their zeros

A marginal merit of our use of the square-well-shaped $`T(p)`$ lies in the availability of the explicit wave functions. For the lack of space we have to omit illustrative pictures, mentioning just a few of their most characteristic features.

In the first step we notice that in the rightmost interval of $`p`$ the absence of any nodal zero in the wave function is in fact very similar to the usual Sturm Liouville behavior. Less expectedly, at the exact energy value one encounters an infinity of the nodal zeros in the leftmost subinterval of $`p(\mathrm{},1)`$. In this domain we are fortunate in studying the exactly solvable case. The very presence of this infinite “left” set of nodal zeros is extremely sensitive to the numerical level of precision we use. Indeed, the errors are proportional to the unphysical $`\psi ^{(unphys.)}(p)exp(2\beta p)`$ which is growing rapidly at $`p1`$.

After the smallest deviation of the energy $`E`$ from its absolutely precise bound-state value even the non-numerical and absolutely precise wave functions will be dominated by the growing asymptotics $`\psi ^{(unphys.)}(p)exp(2\beta p)`$ near the left infinity. The change of sign of the asymptotics is a reliable source of information about the fact that the energy crossed it physical value. This observation survives in the shooting numerical algorithms  as well as in the rigorously proved versions of the method of Hill determinants ). In this context our present numerical experiments could be perceived as opening a number of new questions. Some of them emerge in the purely numerical context of an appropriate generalization of the Prüfer-type algorithms. Especially in the vicinity of the correct physical energies they could lead to reliable and robust right-to-left shooting numerical recipes.

## 12 Outlook

In the $`x`$representation of our problem our main emphasis has been put on the exact solvability of its replacement by the purely imaginary square well model. New light has been thrown on some properties of wave functions. One can expect that the further detailed study of the $`𝒫𝒯`$ symmetric square wells will give new answers to the puzzles concerning the irregular behavior of the nodal zeros in the complex plane as formulated in ref. . Our present study indicates that some complexified versions of the Sturm Liouville oscillation theorems should be developed for the study of zeros of the separate real and imaginary parts of $`𝒫𝒯`$ symmetric wave functions.

After the standard Fourier-transformation transition to the $`p`$representation of our imaginary cubic oscillator the underlying eigenvalue problem can be seen from a different perspective. Its Hamiltonian is being replaced by a real differential expression. On a suitable Hilbert space this specifies the Hamiltonian operator with the numerical range (and, hence, spectrum) which is, obviously, real. This complements the extensive discussion of this topic in . Among several immediate constructive consequences of the latter observation we underlined the consistency of the approximations imposed directly upon the kinetic term $`T(p)`$.

As long as the behavior of wave-function asymptotics at large $`|p|1`$ differs in the left and right infinity, several new qualitative aspects of the problem emerge and became clarified by our schematic piece-wise constant approximation of $`T(p)`$. At small $`p0`$ the emergence and motion of the nodal zeros can be interpreted in a graphical manner explaining some features of the $`N(Z)`$ dependence. In particular, the puzzling loss of its monotonicity seems confirmed by our solvable model.

The use of the momentum representation proved able to throw a new light on the counterintuitive bound states in $`𝒫𝒯`$ symmetric quantum mechanics. The emergence/disappearance of our quasi-degenerate doublets should be emphasized as, perhaps, analogous to the unavoided level crossings in harmonic oscillators  and/or to the anomalous doubling of levels in the models of Natanzon type . Similar irregularities in the spectra could be, perhaps, attributed to a peculiar combination of the analyticity and non-Hermiticity in $`𝒫𝒯`$ symmetric systems.

In a brief summary of our numerical observations let us point out the regularity of the $`Z`$dependence of the number $`N(Z)`$ of the bound states. This indicates that one should search for an improved application of Sturm-Liouville theory in complex domain . The possibility of deduction of new oscillation-type theorems exists, first of all, in the middle interval of $`p(1,1)`$ where, for a continuously growing energy parameter $`E`$, a steady right-ward movement of the nodal zeros competes with the exponential terms which are varying slowly.


format: text

text:
# Luminosity and mass function of galactic open clusters: I. NGC 4815 Based on observations made at the European Southern Observatory, La Silla, Chile.

## 1 Introduction

Galactic open clusters are loose and weakly bound systems, which only in exceptional cases survive more than $`10^8`$ yrs berg80. Their mean lifetime is determined essentially by two factors: the birth-place (or the orbit) and the initial mass. The dissolution of open star clusters is mainly due to close encounters with giant molecular clouds which are very frequent in the galactic thin disk theu92a; theu92b. For this reason, the age distribution is strongly biased towards young objects wiel71. Out of about 1,200 open clusters which are known in the thin disk, only 80 are older than the Hyades frie95.
The dynamical evolution of open star clusters in the Galactic disk environment has been studied by several authors terl87; theu92a; theu92b To summarize, it has been shown that encounters with dense molecular clouds are catastrophic, whereas diffuse clouds do not completely destroy the cluster. Moreover, tidal effects induced by the general Galactic gravitational field disturb significantly cluster evolution, producing distorted debris, at the center of which a few closely bound stars remain, with a large fraction of binaries.
On the other hand, the internal dynamics of open clusters has been investigated by means of N-body models aars96; fuen97, which include the effects of stellar evolution and a fine treatment of binaries and/or multiple stars.
Among the others, a very important result is that the cluster evolution crucially depends on the Initial Mass Function (IMF), and on the cluster richness. To go in some detail, de la Fuente Marcos fuen97 finds that power-law salp55 like IMFs accelerate the cluster disruption if the initial population is small, and that at the end of the evolution a huge number of binary systems with the same mass is formed.

An important constraint for all these models would be the possibility to derive the IMF from observations of open star clusters. Unfortunately, this is a very difficult task for several reasons.
First, most open clusters are poorly populated ensembles of star, containing from tens to hundreds of stars, and only in rare cases up to some thousands. Secondly, they are located in the inner regions of the thin disk. This fact renders it very difficult to detect cluster members (at odds with globular clusters), due to the high contamination of the field stars in the same line of sight as the clusters. Third, they can be strongly obscured due to interstellar absorption lying between us and the cluster jane82.

Probably, the most challenging task is to obtain cluster membership. This is feasible for objects close to the Sun, for which we can obtain high quality radial velocity and proper motion measurements. Going farther away one has to rely on statistical corrections which usually are made by comparing the cluster with the field in its outskirts. Ideally, when a correction for field star contamination has been performed, a reasonable actual (Present Day) mass function (PDMF) can be derived from the observed LF. By modeling the dynamical evolution of the cluster the PDMF can be converted into an IMF which can be compared with the MFs of other clusters to look for universality or deviations.

In this paper, the first of a series dedicated to obtaining luminosity and mass functions for Galactic open clusters, we present results for NGC 4815, a rather rich, distant open cluster located at about 2.5 kpc from the Sun .
The ultimate goal is to provide a sample of PDMFs as a function of cluster age, population, and location to serve as template for N-body models.
We plan to study other clusters with the same method in forthcoming papers. In order to obtain a good representation of the field stellar population, in addition to the cluster we observed four nearby Galactic fields, whose population has been averaged and used for the statistical subtraction of the foreground/background objects from the cluster. An adaptive kernel estimator was then applied to obtain a continuous LF, free of binning bias, together with the standard histogram technique. Finally, the MF has been derived from $`m=4m_{}`$ down to $`m=0.35m_{}`$ by using the luminosity-mass relation from the Girardi et al. gira00 models.

The layout of the paper is as follows. Section 2 summarizes our present knowledge of the open cluster NGC 4815. In Section 3 we present the observations and data reduction technique; Section 4 is dedicated to the discussion of the CMDs obtained for the cluster and the field. Section 5 is dedicated to a detailed analysis of the star counts, while Section 6 deals with the LFs derived with the two different techniques. Section 7 presents the MF of NGC 4815 and the comparison with other MFs. Finally Section 8 summarizes our results.

## 2 Previous investigations

NGC 4815 has been studied in the past by several authors.
The first photometric data have been collected by Moffat & Vogt moff73, who studied 9 stars in the cluster region, but did not recognize any sequence.
Kjeldsen & Frandsen kjel91 presented a catalogue of UBV photometry for 559 stars.
Carraro and Ortolani carr94 firstly made an accurate analysis of the CMD by obtaining BV photometry of 2500 stars down to V=19.5. Their results led to the determination of the cluster parameters: NGC 4815 turns out to be of intermediate age (about half a Gyr old, like the Hyades), with a reddening $`E(BV)0.70`$, and a distance of about 2.5 kpc from the Sun. The metallicity is half Solar (see Table 1). From a preliminary analysis of the MS mass function they obtained an index $`\alpha =2.55`$, slightly higher than the classical Salpeter value.
These data have been used to study the cluster structure and luminosity function with a different method by Chen et al. chen98, who obtained for the first time the cluster mass function by subtracting the field contamination with a Galaxy structure model. In the mass range $`1.1<\frac{m}{m_{}}<2.2`$ they obtain a slope $`\alpha =1.97\pm 0.17`$, and a total mass of 912 $`m_{}`$ to the photometric limit (V = 19). They finally found evidence for mass segregation.

## 3 Observations and reduction

### 3.1 Observations

NGC 4815 has been observed with the ESO-NTT telescope on April 14-16, 1999. We used the EMMI red arm, equipped with the 2048$`\times `$2048 ESO$`\mathrm{\#}36`$ Tektronics CCD. The scale on the sky is 0.27 arcsec/pixel, for a total field of $`9.2\times 9.2`$ arcmin2. The observing log is in Table 1. Almost all the images have been collected in sub-arcsec seeing conditions and in photometric conditions. During the same nights, four Galactic fields have been observed. These fields are located at 30 arcmin from the center of NGC 4815, on the North, South, East and West sides, respectively. The exposure times and the seeing of the Galactic field images are the same as for the cluster data, with the only exception of the southern one, for which only short exposures could be obtained. We opted to exclude this field in the statistical analysis of luminosity and mass function.

### 3.2 Data reduction and calibration

All the images have been pre-processed in the standard way with IRAF, using the sets of bias and sky flat-field images collected during the observing nights. The instrumental magnitudes, and the positions of the stars for each frame have been derived by profile-fitting photometry with the package DAOPHOT II and ALLSTAR stet87. Then we used ALLFRAME stet94 to obtain the simultaneous PSF-fitting photometry of all the individual frames. In order to obtain the transformation equations relating the instrumental ($`vi`$) magnitudes to the standard V (Johnson), $`I`$(Kron-Cousins) system we followed the procedure already described in Bed20. Seven Landolt (1982) fields of standards have been observed. Specifically: the PG0918+029, PG0942-029, PG1047+003, PG1525-071, PG1657+078, SA104334, and SA110-496 regions. In each of these, there are other secondary standard stars by Stetson stet00, which extend the previous Landolt sequence. For all of these stars, aperture photometry was obtained on all the images after the digital subtraction of neighboring objects from the frames. We used DAOGROW stet90 to obtain the aperture growth curve for each frame and other related codes (by PBS) to determine the aperture corrections, the transformation coefficients, and the final calibrated photometry (e.g., Stetson 1993). We used transformation equations of the form:

$$v=V+A_0+A_1X+A_2(VI),$$

$$i=I+B_0+B_1X+B_2(VI),$$

where the values of the coefficients are listed in Table 3 for the two nights. In these equations $`v`$ and $`i`$ are the aperture magnitude already normalized to 1 sec, and $`X`$ is the airmass. Second order color terms were tried and turned out to be negligible in comparison to their uncertainties. It is a reasonable hypothesis that the color-dependent terms imposed by the telescope, filter, and detector should be reasonably constant over the course of a few nights, so after $`A_2`$ and $`B_2`$ had been independently determined from each night’s data, we obtained weighted mean values of these parameters and imposed them as fixed quantities for the definitive calibration. The NGC 4815 images enabled us to define local standards that met the following conditions: each was well separated from its neighbors (according to a criterion described in Stetson stet93 §4.1), each was observed in all frames, each had a mean value of the goodness-of-fit statistic, CHI, less than 1.5. Once calibrated, this local standard sequence enables us to place all of observations of NGC 4815, and of the Galactic fields on the same magnitude system with an uncertainty less than 0.005 magnitudes (internal error).

### 3.3 Artificial star tests

In order to obtain the LF, we estimated the completeness of our samples. The completeness corrections have been determined by standard artificial-star experiments on our data. For each field, we performed five independent experiments adding on the original frames as large a number as possible of artificial stars in a spatial grid as in Piotto & Zoccali piot99. We performed the photometry of the frames with the added artificial stars following all the steps and using the same parameters as for the original images. We converted the new instrumental magnitudes to the standard photometric system using the coefficients of the calibration equations of Section 2.2.

In the LFs we include only the values for which the completeness corrections, defined as the ratio between the number of found artificial stars to that of the original added ones, were 50% or higher. This sets the limiting magnitude to V = 24.

## 4 The Color-Magnitude Diagrams

We obtained the V and I photometry with the aim of determining the LF (and MF) of the main sequence stars in NGC 4815. We used two colors as we wanted also to derive the CMD, which allows to: i) discriminate stars from false detections; ii) better identify the cluster population, in particular the MS stars.
The CMD for NGC 4815 is shown in Fig. 1, whereas the CMDs for the four nearby fields are presented in Fig. 2. Clearly, the CMD of the cluster field does not differ much from those of the surrounding Galactic disk. This means that NGC 4815 is actually badly contaminated by background/foreground stars, as expected from its particular position in the Galactic thin disk. The cluster field (Fig. 1) is much more populated than the neighboring fields (Fig. 2), confirming that NGC 4815 is a real open cluster (see also the discussion in Section 5). Our photometry reaches $`V25`$ at the base of the MS.
The MS extends almost vertically for more than 10 magnitudes, from V = 15 to V=25, although the completeness analysis prevents us from using stars dimmer than V= 24. The width of the main sequence gets greater at larger magnitudes. This is partly due to the increasing photometric errors at increasing magnitudes (see Fig. 1). However the MS is much wider than expected simply from photometric errors. Between the various causes which concur to enlarging the natural MS width we can envisage the presence of unresolved binary stars, common in open clusters, and a possible spread in metals. The red vertical sequence, which detaches from the MS at about V = 22.5 represents the giant branch of the Galactic disk population.
At odd with the previous photometry, the present data (see Fig. 1) do not allow us to recognize clearly the presence of a red giant clump, better visible in the BV photometry of Carraro and Ortolani carr94 at $`V14.3`$. This is due to the saturation of many bright stars above V = 14.5, which were not of interest for the purposes of this work.

## 5 Star counts and cluster size

The aim of this section is to obtain the surface density distribution of NGC 4815, and derive the cluster size and extension in the magnitude space by means of star counts. The cluster radius is indeed one of the most important cluster parameters, useful (together with cluster mass) for the determination of cluster dynamical parameters. Star counts allow to determine statistical properties of clusters (as visible star condensations) with respect to the surrounding stellar background.

Star counts have been performed following the methods described in dan86, sel94 and sele98 by evaluating the functions $`N(r)`$ and $`N(V)`$ ($`N(I)`$). $`N(r)`$ represents star number in the circle of radius $`r`$ and $`N(V)`$ is star number up to apparent magnitude $`V`$, i.e. integral luminosity function. In order to obtain $`N(r)`$ we divided the cluster field (and comparing fields) in concentric rings 40 pixels (11 arcsec) wide and counted the number of stars falling inside each ring. Then $`N(r)`$ is obtained by summing the corresponding ring star numbers. Similarly, in order to obtain $`N(V)`$, we divided stars in the various fields in magnitude bins $`0.2^m`$ wide and summed up corresponding bin star numbers. Star counts with $`N(r)`$ give us cluster radius and cluster star number, star counts with $`N(V)`$ give us the magnitude of faintest cluster stars (up to the limiting magnitude).

Due to the saturation of the brightest stars ($`V<14.5`$), some regions of the cluster cannot be used to perform star counts. In these cases, for $`N(r)`$ determination we used the stars located in the same regions from the work of Carraro & Ortolani (1994).

A smoothed surface density $`F(r)`$ is then obtained by differentiating numerically the best fit least squared polynomial approximation for $`N(r)`$ according to the equation:

$$N(r)=2\pi _0^rr^{}F(r^{})𝑑r^{}$$
(1)
(see Seleznev (1994), for details). The results are shown in Fig.3, where the logarithm of the surface density profile is plotted as a function of the limiting magnitude. From this figure it is evident that the cluster significantly emerges from the background only for $`V`$ less than 21. Below this value NGC 4815 is completely confused with the Galactic field. Moreover we note a significant density decrease at the cluster position.

Integral distribution functions are then obtained both for cluster field and for comparing fields. Let us consider below a theoretical integral distribution function $`N_i(x)`$, where $`x`$ can be $`r`$, $`V`$ and $`I`$, and $`i=0`$ indicates cluster field, while $`i=1,2,3,\mathrm{}`$ indicate comparing fields.

Following dan86 and sele98 we consider the normalized deviation $`z`$ of $`N_i(x)`$ from the corresponding average number of field stars:

$$z_i(x)=\frac{\left|N_i(x)\overline{N(x)}\right|}{\sigma (x)}(i=0,1,\mathrm{},k),$$
(2)
where

$$\overline{N(x)}=\frac{1}{k}\underset{i=1}{\overset{k}{}}N_i(x),$$
(3)
and

$$\sigma ^2(x)=\frac{1}{k1}\underset{i=1}{\overset{k}{}}\left[N_i(x)\overline{N(x)}\right]^2.$$
(4)
The normalized deviations $`z_i`$ are connected with the probability that the i-th field significantly deviates from other fields. Then $`z`$ values are usually called ”significance levels”.

There are two criteria to determine cluster boundaries in $`x`$-space. The first criterion is as follows. One has to look for a value $`x_c`$ such that $`(N_0(x)\overline{N(x)})=N_c=const`$ for all $`xx_c`$ and $`z_0(x)z_i(x)`$ ($`i=1,2,3,\mathrm{}`$) for all $`xx_c`$. In this case we can conclude that the cluster has the boundary at the point $`x=x_c`$ and $`N_c`$ is cluster star number, and we can estimate an error for $`x_c`$ and $`N_c`$ (Danilov et al. 1986).

The second criterion is as follows. The cluster has no boundary in the sense of the first criterion, but for $`x>x_{c2}`$ $`z_0(x)<z_i(x)`$ for some $`i`$, i.e. i-th field deviates from mean $`\overline{N(x)}`$ more significantly than cluster field. This situation means that we can determine only a lower estimate of the cluster boundary, i.e. $`x_c>x_{c2}`$ and $`N_c>N(x_{c2})`$. The reason is related to large-scale fluctuations of integral distribution function in the comparing fields. Unfortunately this situation frequently occurs for open star clusters located in a strongly fluctuating background.

The star counts in $`V`$ and $`I`$ spaces have been performed both for central part of NGC 4815 (circle with $`r2arcmin`$ in the central field, see Table 2) with twelve comparing fields (four circles in each East, West and North fields) and for all cluster field (i.e. the entire central field) with three comparing fields (East, West and North).

The star counts for the central part of the cluster confirm that it is less populated by faint stars ($`V>18^m20^m`$) than the other fields (see Fig.3). We obtained cluster boundary $`V_c=20.8^m\pm 0.2^m`$ following the first criterion (see above). The trend of $`z_0(V)`$ is shown in Fig.4, where it is compared with mean $`z`$ for North, West and East fields. The luminosity function of cluster stars (see following section), that is differential distribution function, has negative value for $`V>V_c`$. Presumably, the reason is the higher light absorption along line of sight beyond the central part of the cluster.

By performing star counts in $`V`$ and $`I`$ for all the central field we did not obtain the cluster boundary according to first criterion. In the $`I`$ passband $`N_0(I)`$ significantly deviates from mean $`\overline{N(I)}`$ up to $`I22^m`$, but in the $`V`$ $`N_0(V)`$ significantly deviates from mean $`\overline{N(V)}`$ only up to $`V20^m`$. However $`(N_0(V)\overline{N(V)})`$ and $`(N_0(I)\overline{N(I)})`$ increase for much fainter magnitudes, and the corresponding luminosity functions are positive (see Fig.5 and 6 in the following section). This contradictory result is due to the very complicate character of the stellar background in the vicinity of NGC 4815 and beyond it. It means that mean integral distribution function $`\overline{N(V)}`$ in comparing fields increases more slowly than $`N_0(V)`$, but the dispersion $`\sigma (V)`$ is too large.

Because the cluster center is slightly shifted with respect to the center of the central field, star counts in $`r`$-space can be carried out in the entire circle only to radius $`r=880pixel`$ (4 $`\mathrm{𝑎𝑟𝑐𝑚𝑖𝑛}`$). Due to this reason star counts were made also in sector of $`199^{}`$ to radius $`r=1080pixel`$ (4.9 $`\mathrm{𝑎𝑟𝑐𝑚𝑖𝑛}`$). The results of star counts for different limiting magnitudes are listed in Table 4.

First column contains limiting magnitude, second and third columns contain cluster radius and cluster star number from counts in the entire circle, fourth and fifth columns contain the same values from counts in $`199^{}`$ sector. Star number in fifth column is obtained as result of star counts multiplied by (360/199).

The comparison of results of star counts in the entire circle and in the sector shows that the asymmetry of the cluster increases at increasing limiting magnitude. Beginning from $`V_{lim}=21^m`$ the cluster disappears as a statistically significant stellar density fluctuation, and below $`V_{lim}=24^m`$ the cluster region appears as stellar density minimum - statistically significant comparing with three fields. The most probable reason is the presence of absorbing cloud just beyond the cluster.

Star counts with stars up to $`V_{lim}=19^m`$ gave us the estimate of cluster radius, which turns out to be $`3.6\pm 0.3`$ arcmin, significantly lower than $`4.6`$ $`\mathrm{𝑎𝑟𝑐𝑚𝑖𝑛}`$ estimated by Carraro & Ortolani (1994). However, result of Carraro & Ortolani (1994) can be regarded only as preliminary, because cluster radius was determined only by surface density profile as distance to a point where the slope of this curve changes, being the slope negative along all the profile. It is difficult to compare our results with previous results directly, because now we are using a quantitative statistical method. Unfortunately, the cluster radius of $`3.6`$ $`\mathrm{𝑎𝑟𝑐𝑚𝑖𝑛}`$ is very close to the limit of investigated field, and we can not regard this result as definitive. Furthermore, $`3.6`$ $`\mathrm{𝑎𝑟𝑐𝑚𝑖𝑛}`$ corresponds to a radius of $`2.6pc`$ at the distance of $`2.5kpc`$. This value lies in the lower limit of cluster radii interval according to dan94. Then it is quite possible that we have determined the radius of the cluster core.

It is necessary to note that the central field contains much more stars than any comparing field (see Fig.1 and Fig.2). This excess in star number is mainly due to the periphery of the central field which is very rich in faint stars. Unfortunately our star counts can not unambiguously answer the question - is the cluster really disappearing at $`V_{lim}=21^m`$ and has radius of about 3.6 arcmin or faint stars in the outskirts of central field belong to cluster and form its halo ?

Firstly, star counts in wider fields are necessary (both for cluster and comparing fields). Desirable radius of investigated field is of about 10-15 pc (i.e. 14-21 arcmin). Secondly, it is necessary to study color extinction for stars of different magnitudes, i.e. to observe in three colors.
On the basis of the previous discussion, we are going to estimate the luminosity function for all the central field in the magnitude range $`V[14^m;24^m]`$.

## 6 Luminosity Function

The luminosity function of NGC 4815 was obtained as follows. In the CMDs, two lines, to the left and right of the main sequence, defined by the two relations $`V=6.8(VI)+12`$ and $`V=6.8(VI)+3.84`$, were used in order to select only the main sequence stars (Fig. 1). For each field, an adaptive kernel estimator was then applied to obtain a continuous LF. The LF was also derived by the standard histogram technique. A detailed description of the application of the kernel function method for estimating distribution functions can be found in silv86, and in sele98 and sele00.

The advantage of this method is twofold: 1) first of all it is not necessary to define a bin, whose size usually is arbitrary and does not allow to detect unambiguously all the statistically significant features present in the LF; the estimate of the LF by a kernel method with its confidence interval is an estimate of an unknown density of a stellar population, while the histogram alone gives only the distribution of a specific sample; 2) the form of a histogram - crucial for the determination of the LF slope - depends both on the bin width and on the initial bin location. The distribution function provided by the kernel estimate only depends on the kernel width. Moreover criteria have been developed to fix the kernel width with high degree of confidence silv86.

The adaptive kernel estimate of the density distribution in one-dimensional case is calculated by the formula:

$$\widehat{f}(x)=\frac{1}{n}\underset{i=1}{\overset{n}{}}h\lambda _iK\left(\frac{xX_i}{h\lambda _i}\right),$$
(5)
where n is the size of the observing sample $`X_1\mathrm{}X_n`$ of the random variable $`x`$, $`K`$ is the kernel function which must satisfy the condition:

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}K(x)𝑑x=1,$$
(6)
and $`\lambda _i`$ are the local window width factors. The adopted kernel has the form

$$K(x)=\frac{3}{4\sqrt{5}}\left(1\frac{x^2}{5}\right)$$
(7)
when $`|x|\sqrt{5}`$, and zero elsewhere.

Basically eq.5 spreads each observed value $`X_i`$ of the random variable $`x`$ over an interval (window) centered on $`X_i`$ with the probability density $`K\left(\frac{xX_i}{h\lambda _i}\right)`$. The kernel width $`h\lambda _i`$ can differ from point to point. The advantage of having an adaptive kernel is that it is easy to properly sample the tails of the distribution, just playing with the kernel width, which is set larger in low-density regions automatically (see silv86). Parameter $`h`$ was set to be $`0.3^m`$.

To determine the confidence interval for the resulting LF, we used the method described in Seleznev et al. (2000) and references therein.

For each field, the LFs have been corrected for incompleteness. The LFs of the external fields were used to subtract the field-star contamination from the LF of the central field.

To this end we defined $`\widehat{f}_i(V)`$ as the density of the distribution of the stellar $`V`$ magnitudes, where $`i`$ is the number of the field (here we assume that $`i=0`$ corresponds to the cluster field or central field and $`i=1,2,3`$ correspond to the external fields). The luminosity function of the cluster stars is given by

$$\widehat{f}(V)=\widehat{f}_0(V)\frac{1}{k}\underset{i=1}{\overset{k}{}}\widehat{f}_i(V).$$
(8)
The same procedure was repeated to get the $`I`$ LF. The $`V`$ and $`I`$ LFs obtained both with the histogram and the adaptive kernel techniques for all the central field are shown in Figs.5 and 6. Note the large uncertainties below $`V=21`$ ($`I19`$), where fluctuations of stellar background density play a significant role. LFs obtained by different methods are close each other due to large volume of the sample. Note that the stellar content below $`V=20`$ ($`I=18`$) is due to outskirts of the central field (see previous section).

## 7 Mass Function

The LFs can be transformed into MFs using a mass-luminosity relation (MLR). Since we cannot obtain an empirical transformation, we must rely on the theoretical models. Therefore, we used the ZAMS isochrones by Girardi et al. (2000) with a metallicity Z=0.008. The distance modulus and the reddening that result from the fit with the model are in agreement, within the errors, with the values found by Carraro and Ortolani (1994). We adopted $`E(BV)=0.70`$ and $`(mM)_V=14.10`$. From these values we get $`(mM)_I=13.20`$ (see Table 1). Fig. 7 shows the MF from the V LFs obtained by the adaptive kernel estimator and the standard histogram technique. The two MFs are rather similar. There is also a substantial agreement between the MFs from the V and I LFs (Fig. 8), showing an internal consistency of both the models and the method.

By looking at Fig. 7, it appears that there is a dip in MF at $`m0.9M_{}`$. This dip in MF is a the consequence of the corresponding dips in LFs at $`V19.5^m20^m`$ and $`I17.5^m18^m`$ (see Fig.s 5 and 6). To estimate the statistical signifincance of this dip, we concentrate on ”Adaptive kernel estimator” method, which is is more sensitive to peculiar properties of LF and MF than the histogram technique (see previous sections).
We can use the ratio $`\eta `$ of the depth of this minimum (relatively to the mean smooth MF level) to the dispersion of MF estimate in the vicinity of MF minimum as a measure of the statistical significance of the MF minimum. As for the MF minimum at $`m0.9M_{}`$ we obtain $`\eta 1.98`$ from the ”adaptive kernel” MF estimate, which corresponds to a reality probability of $`95\%`$. On the other hand, $`\eta 0.78`$ comes out of the ”histogram” MF estimate, which corresponds to reality probability of $`56\%`$. Then we can conclude that this dip of MF seems to be statistically significant.

Finally, it is possible to obtain an estimate of the cluster mass by integrating the luminosity function. In detail,

$$M=<m>\times N=N\underset{M_{V_{min}}}{\overset{M_{V_{max}}}{}}m(M_V)\varphi (M_V)𝑑M_V$$
(9)
where $`M_V`$ is absolute stellar magnitude $`M_V=V(mM)_V`$, $`\varphi (M_V)`$ is the normalized luminosity function

$$\varphi (M_V)=f(M_V)/N,$$
(10)
N is the total number of stars, $`m(M_V)`$ is the adopted mass-luminosity relation. The luminosity function and mass-luminosity relation have been approximated by a spline. Limits in the integral are $`M_{V_{min}}=0.1`$ and $`M_{V_{max}}=9.9`$. The integral turns out to be

$$M1620\pm 690M_{}$$

for kernel LF estimate and

$$M1820\pm 650M_{}$$

for histogram LF estimate.

This mass estimate refers to total central field. Unfortunately it is very difficult to evaluate what part of total cluster mass is constituting by this estimate. It depends on both unknown cluster size and unknown upper limit of magnitude. If we assume that NGC 4815 has a limiting magnitude of $`V_c=20.8^m`$, we get the mass estimate

$$M880\pm 230M_{}$$

.

The error is determined from mass estimate for upper and lower boundaries of confidence interval for LF (see Fig.5).

## 8 Discussion and conclusions

The MF can not be represented by a unique power law. In the mass range $`2.5<m/m_{}<0.8`$, the PDMF is reproduced by a power-law with $`\alpha =3.1\pm 0.3`$ ($`\alpha =2.35`$ for the Salpeter MF in this notation). In the same mass range the slope of the MF from the $`I`$ LF is $`\alpha =2.9\pm 0.3`$.
Below $`0.8m_{}`$, the MF is much flatter ($`\alpha =0.68\pm 0.49`$).
It is difficult to correctly interpret the PDMF in this low mass regime. However we believe that the flattening of the PDMF below $`0.8m_{}`$ can be due to basically two effects.
First of all, we have already noted how the cluster vanishes (but not disappears) into the Galactic field below $`V=21`$, i.e. $`m=0.8m_{}`$. This is likely originated by irregularities in the stellar background. As a consequence low mass stars are preferentially missed due to absorption.
Secondly, we cannot however exclude cluster dynamical evolution, which causes a loss of preferentially low mass stars. Also, stripping of low mass stars can occur along the cluster orbit, due to the interaction of the cluster with the dense Galactic field. NGC 4815 is indeed located close to the Sagittarius-Carina spiral arm, rich in HII regions and molecular clouds. Estimating the relative importance of these two effects is very difficult and goes beyond the aim of this paper.
Therefore, the slope of the MF below 0.8 $`m_{}`$ can not be considered representative of the cluster IMF.

On the other hand, whether the PDMF above 0.8 $`m_{}`$ is representative of the IMF is a matter which needs a careful model of the cluster too.
A recent review on the IMF derived from the star counts is in Scalo scal98. From the comparison of several studies, he found that the spread in the slope at all masses above about $`1m_{}`$ is too large for adopting some average value, and in particular, there is no evidence for a clustering of points around the Salpeter value (2.35). Scalo scal98 concluded that either the systematic uncertainties are so large that the IMF cannot yet be estimated, or that there are real and significant variations of the IMF slope. For sure, the dynamical evolution, which differs from cluster to cluster, can modify the MF slope, and also this effect contributes to the MF slope dispersion. It has emerged that the MF can not be reproduced with a single power low: at low-masses the IMF from the field stars, open and globular clusters piot99 is relatively flat, and it becomes steeper above about $`1m_{}`$. Also NGC 4815 seems to show this break in the MF shape at $`m1m_{}`$ (Fig. 7), though it is not possible to disentangle the real differences in the IMF from the effects of the dynamical evolution. A final point is the effect of unresolved binaries and mass segregation. As discussed by Scalo (1998, and references therein), unresolved binaries and mass segregation can only flatten the apparent IMF derived by not accounting for these, so the real IMF must be steeper than the quoted slope by an unknown amount.

In Fig. 9, the MF of NGC 4815 is compared with the MFs for other 15 open clusters. Most of them cover masses above $`1m_{}`$. We have derived the MFs directly from the published LFs using the same set of isochrones from Girardi et al. gira00 used to get the NGC 4815 MF. The MFs in Fig. 9 are shifted by arbitrary constants for display purposes. The MF of NGC 4815 is one of the most extended (in mass) MFs obtained so far for an open cluster. Only for the better known and nearby open clusters (Pleiades and M35), and for objects observed with HST (NGC 2516, NGC 2477, and NGC 2420) do the MFs reach $``$ 0.1 $`m_{}`$.

Fig. 10 shows the slopes $`\alpha `$ of the MFs plotted in Fig. 9 as a function of the mean mass. The slopes have been obtained by fitting a power law in the $`m<1m_{}`$ and $`m>1m_{}`$ mass ranges. Clusters with ages greater and less than or equal to 0.5 Gyr are plotted with different symbols. As noted by Scalo, the spread in slope is very large. However, it is interesting to note that all the clusters with ages $`0.5`$ Gyr have very similar slopes ($`\alpha =2.25\pm 0.88`$) for $`m>1m_{}`$.

NGC 4815, which has an age of 0.5 Gyr (Carraro and Ortolani carr94), has a MF slope fully consistent with this value. The MFs of the clusters older than 0.5 Gyr have a similar or flatter slope. This may indicate that their PDMF may have not or may have been changed by the dynamical evolution which has had more time to operate. A dynamical model of each of these clusters is necessary to confirm this interpretation of Fig. 10. At the moment, the similarity of the MF slopes for all the youngest clusters in Fig. 10 can be interpreted as an evidence of a similarity of the IMF, which seems to not depend much on metallicity, stellar density, or galactocentric radius.

The interpretation of the MF slopes below $`0.8m_{}`$ is more difficult, both because there are less (and more uncertain) data, and because we expect that in this mass regime the dynamical evolution should be more effective in altering the IMF. Still, all the MFs seems to have the same flat general trend. We also note that these slopes represent an upper limit for the slopes of the MF of the Galactic globular clusters (Piotto and Zoccali piot99).

###### Acknowledgements.

We acknowledge the financial support by the Ministero della Ricerca Scientifica e Tecnologica (MURST) under the program Treatment of large field astronomical images. PBS an AFS gratefully acknowledge the generosity of the Universitá di Padova and MURST for supporting their visit to the Dipartimento di Astronomia and the Osservatorio di Padova.


text:
# CMB Constraints on a Baryonic Dark Matter-Dominated Universe

## 1. Introduction

Since the pioneering work of Zwicky (1933), astronomers have found compelling evidence to suggest that the major contribution to the overall mass density of the universe is in the form of dark matter. Even if its precise nature remains unknown (see e.g. Bergstrom (2000) for a recent review), the currently favoured hypothesis of cold dark matter (CDM) is in good agreement with an enormous body of data: anisotropies in the cosmic microwave background (CMB), large-scale structure surveys, cluster abundances, structure of the Lyman-$`\alpha `$ forest and a host of other measurements (see e.g. Bahcall et al. (1999)). Furthermore, in the past 30 years, upper limits on the CMB anisotropies and predictions made by theories of galaxy formation have been extremely effective at ruling out models based on purely baryonic or hot dark matter (see e.g. Doroshkevich, Zel’dovich & Sunyaev (1978); Wilson & Silk (1981); Kaiser (1983); White, Frenk & Davis (1983)).

The quality of astrophysical data has improved and it currently seems likely that the present fractional overall energy density $`\mathrm{\Omega }`$ includes contributions from a cosmological constant $`\mathrm{\Lambda }`$ ($`\mathrm{\Omega }_\mathrm{\Lambda }`$) and from cold ($`\mathrm{\Omega }_{cdm}`$), hot ($`\mathrm{\Omega }_{hdm}`$) and baryonic ($`\mathrm{\Omega }_b`$) dark matter. Also, it appears that the densities of these components are all within one or two orders of magnitude of each other. While the similar densities of dark ‘energy’ and ‘matter’ can possibly be explained by advocating a quintessence scalar field which ‘tracks’ the matter density (Steinhardt, Wang & Zlatev (1999)), the reason why the densities of the $`3`$ matter components are so similar is not so apparent.

Predictions made by structure evolution theories based on the $`\mathrm{\Lambda }`$CDM model disagree with some galaxy observations, which suggests that, for the model to hold, assumptions made about the properties of CDM particles must be modified (see e.g. Spergel & Steinhardt (1999); Burkert (2000); Madsen (2000); Moore et al. (2000)). If this reasoning is adopted, currently discussed versions of the present cosmological scenario, incorporating a fine-tuned cosmological constant and dark matter particles with ad hoc properties, lose the compelling aesthetic simplicity of the original CDM model (see Sellwood & Kosowsky (2000) for a recent review). Of course, other, more astrophysical, resolutions may be possible Binney, Gerhard and Silk (2001).

The recent CMB anisotropy measurements from the BOOMERanG-98 (de Bernardis et al. (2000)) and MAXIMA-1 (Hanany et al. (2000)) experiments have provided new insights into the cosmological parameters. For the first time, the anisotropy angular power spectrum has been measured over a wide range of angular scales, from multipole $`l50`$ up to $`l800`$ with errors of the order of $`10\%`$. The data sets confirm that there is a peak in the angular power spectrum at $`l200`$ with a steep decline in power from $`l200`$ to $`l300`$. While the presence of such a peak is consistent with the predicted acoustic oscillations in the adiabatic inflationary scenario, the absence of prominent secondary peaks after $`\mathrm{}300`$ is an unexpected result, suggesting that the value for the physical baryonic density is $`50\%`$ higher than that predicted by Big Bang Nucleosynthesis (BBN) (see e.g. Tegmark & Zaldarriaga (2000); Jaffe et al. (2000)). We note however that this discrepancy with the standard model is no more than 2 $`\sigma `$ when the most recent CMB anisotropy results are included Mason et al. (2001).

In recent months, several solutions have been advocated to reconcile BBN with the larger values of the baryon density suggested by the CMB data. In a recent letter by McGaugh (2000), the simple and controversial solution of considering only baryonic dark matter (BDM) plus a cosmological constant has been proposed. The lower CMB second peak can thereby be explained by the collisional damping expected in a BDM universe, while maintaining the baryon density at a level compatible with the BBN constraint. As is well known, purely baryonic adiabatic models are in severe contradiction with the measured amplitude and shape of galaxy and cluster correlations. However, there remain serious issues about the possibly scale-dependent bias parameter. Hence McGaugh’s result is interesting for various reasons: it resolves the disagreement with BBN by using a simpler model ($`\mathrm{\Omega }_{cdm}=\mathrm{\Omega }_{hdm}=0`$), it demonstrates the model-dependence of theoretical fits to the CMB data by exploring a previously uncovered corner of parameter space and, finally, it questions our theoretical understanding of the structure formation process.

If it turns out that the CDM hypothesis does not work well at some cosmological epochs and scales, then even the exotic scenario of a purely baryonic model with modifications to standard gravity could become competitive, especially in view of the insistence of MOND proponents that this alternative theory of gravity merits serious examination (McGaugh & de Block (1998); Sanders (2000)). Thus, it is timely to study how compatible the CMB data is with a BDM adiabatic universe, what consequences this assumption can have on the remaining cosmological parameters, and the need for CDM from independent observations of galaxy clustering.

In this Letter, we perform an analysis in cosmological parameter space. Using the recent CMB data alone we show that a purely baryonic adiabatic model appears unlikely if the universe is flat ($`\mathrm{\Omega }=1`$). By further combining these results with supernova type Ia (SNIa) data, we rule out baryonic adiabatic models with high confidence.

## 2. Method

The structure of the $`C_{\mathrm{}}`$ spectrum depends essentially on 3 cosmological parameters (see e.g. Hu, Sugiyama & Silk (1998); Efstathiou & Bond (1999); Melchiorri & Griffiths (2000) and references therein): the physical baryonic density $`\omega _b=\mathrm{\Omega }_bh^2`$ and the overall matter density $`\omega _m=\mathrm{\Omega }_mh^2=(\mathrm{\Omega }_{cdm}+\mathrm{\Omega }_b)h^2`$ which define the size of the acoustic horizon at decoupling, and the ‘shift’ parameter $`R`$ related to the geometry of the universe through $`R=2\sqrt{|\mathrm{\Omega }_k|/\mathrm{\Omega }_m}/\chi (y)`$, with $`y=\sqrt{|\mathrm{\Omega }_k|}_0^{z_{dec}}[\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_k(1+z)^2+\mathrm{\Omega }_\mathrm{\Lambda }]^{1/2}𝑑z`$, where the function $`\chi (y)`$ depends on the curvature of the universe and is $`y`$, $`\mathrm{sin}(y)`$ or $`\mathrm{sinh}(y)`$ for flat, closed and open models respectively.

A decrease in $`\omega _{cdm}`$, with $`\omega _b`$ and $`R`$ remaining constant, will reduce $`\omega _m`$. This decreases the redshift of equality causing the peak in the spectrum to be shifted towards smaller angular scales. Constant $`\omega _m`$ can be maintained by increasing $`\omega _b`$. However, an increase in $`\omega _b`$ will decrease the sound speed at decoupling $`c_s1/\sqrt{3(1+3\omega _b/(4\omega _{rad}(1+z_{dec})))}`$ again reducing the acoustic horizon size and shifting the peak to higher $`\mathrm{}`$’s.

To keep the position of the first peak fixed while $`\omega _b`$ and $`\omega _m`$ are varied, $`R`$ has to be appropriately tuned. This can be achieved by increasing $`\mathrm{\Omega }_\mathrm{\Lambda }`$. The shift parameter $`R`$ is sensitive to variations in $`\mathrm{\Omega }_\mathrm{\Lambda }`$ as $`\mathrm{\Omega }_m0`$, and for $`\mathrm{\Omega }_m=\mathrm{\Omega }_b0.05`$ a $`10\%`$ increase in $`\mathrm{\Omega }_\mathrm{\Lambda }`$ can produce a $`50\%`$ decrease in $`\mathrm{}_{peak}`$. Thus, a purely baryonic model can produce the observed peak structure at $`\mathrm{}200`$ only if $`R`$ is increased with respect to the corresponding flat CDM model. We can therefore anticipate that viable purely baryonic models will be $`\mathrm{\Lambda }`$-dominated with closed geometries.

Motivated by these considerations, we compare recent CMB and SNIa observations with a set of models with parameters sampled as follows: $`\mathrm{\Omega }_m=\mathrm{\Omega }_b=0.015,\mathrm{},0.5`$; $`\mathrm{\Omega }_\mathrm{\Lambda }=0.80,\mathrm{},1.04`$ and $`h=0.40,\mathrm{},0.95`$. We vary the spectral index of the primordial density perturbations within the range $`n_s=0.50,\mathrm{},1.50`$ and we rescale the fluctuation amplitude by a pre-factor $`C_{10}`$, in units of $`C_{10}^{COBE}`$. The theoretical models are computed using the publicly available cmbfast program (Seljak & Zaldarriaga (1996)) and are compared with the recent BOOMERanG-98 and MAXIMA-1 results. The power spectra from these experiments were estimated in $`12`$ and $`10`$ bins respectively, spanning the range $`25\mathrm{}785`$. In each bin, the spectrum is assigned a flat shape, $`\mathrm{}(\mathrm{}+1)C_{\mathrm{}}/2\pi =C_B`$.

Following Bond, Jaffe & Knox (1999) we approximate the signal $`C_B`$ inside the bin to be an offset lognormal distribution, such that the quantity $`D_B=\mathrm{ln}(C_B+x_B)`$ (where $`x_B`$ is the offset correction) is a Gaussian variable. The likelihood for a given cosmological model is then defined by $`2\mathrm{l}\mathrm{n}L=(D_B^{th}D_B^{ex})M_{BB^{}}(D_B^{}^{th}D_B^{}^{ex})`$ where $`C_B^{th}`$ ($`C_B^{ex}`$) is the theoretical (experimental) band power, $`x_B`$ is the offset correction and $`M_{BB^{}}`$ is the Gaussian curvature of the likelihood matrix at the peak. We consider $`10\%`$ and $`4\%`$ Gaussian distributed calibration errors for the BOOMERanG-98 and MAXIMA-1 experiments respectively. We also include the COBE data using Lloyd Knox’s RADPack packages. Proceeding as in Dodelson & Knox (2000), we attribute a likelihood to a point in the ($`\mathrm{\Omega }_m=\mathrm{\Omega }_b`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$) and ($`\mathrm{\Omega }_b`$, $`h`$) planes by finding the remaining $`3`$ parameters that maximise it. We then define our $`68\%`$, $`95\%`$ and $`99\%`$ contours to be where the likelihood falls to $`0.32`$, $`0.05`$ and $`0.01`$ of its peak value, as would be the case for a two dimensional multivariate Gaussian.

## 3. Results

In the top panel of Figure 1, we plot likelihood contours in the ($`\mathrm{\Omega }_m=\mathrm{\Omega }_b`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$) plane by applying the maximization/marginalization algorithm described above and using only the BOOMERanG-98 data and limiting our analysis to models with age $`t_0>10G\mathrm{yr}`$. As expected, the low sound speed at decoupling due to the high baryon content makes only closed models compatible with the observations. The deviation from flatness becomes less and less important as $`\mathrm{\Omega }_b0`$, as one would expect from the expression for $`R`$. However, the decrease in the redshift at equality causes positive curvature models to be preferred.

In the parameter range we are sampling, we find that flat models are excluded at $`2\sigma `$. In order to further test the discrepancy with a flat universe we also include CDM models in the analysis with $`\mathrm{\Omega }_{cdm}=0.05,\mathrm{},0.45`$. Comparing our models to the BOOMERanG-98 data alone and restricting the analysis to $`\mathrm{\Omega }=1`$, we find that $`h^2\mathrm{\Omega }_{cdm}>0.04`$ at the $`95\%`$ C.L. Thus flat BDM models are excluded with significance.

For the BOOMERanG-98 analysis with $`\mathrm{\Omega }_{cdm}=0`$, the best-fitting model is a closed model with $`\mathrm{\Omega }_B=0.135`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=1.0`$, $`n_S=0.94`$, $`h=0.45`$, $`C_{10}=0.5`$ and the BOOMERanG-98 calibration left untouched. This provides a good fit to the data (see Figure 2, top panel). Including the MAXIMA-1 data points produces a different best-fitting model: $`\mathrm{\Omega }_B=0.035`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=1.00`$, $`n_S=0.92`$, $`h=0.85`$, $`C_{10}=0.5`$ with a $`12\%`$ upward calibration

FIG. 1.— The likelihood contours in the ($`\mathrm{\Omega }_B`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$) plane, with the remaining parameters taking their best-fitting values. Only the BOOMERanG-98 data is considered here. The contours correspond to 0.32, 0.05 and 0.01 of the peak value of the likelihood, which are the 68%, 95% and 99% confidence levels respectively. The dashed line going from left to right diagonally across the plot indicates flat models. The top panel is with a prior on the age of the universe $`t>10`$ Gyr. The bottom panel includes the BBN prior $`\mathrm{\Omega }_bh^2=0.019\pm 0.002`$ and the SNIa data. The results of these two priors are incompatible at $`3\sigma `$.

for the BOOMERanG-98 points and a $`5\%`$ downward calibration for MAXIMA-1.

As the low value of $`C_{10}`$ already suggests, these best-fitting models are in contradiction with the COBE data. Including COBE in the BOOMERanG-98 analysis does not significantly affect the best-fitting cosmological parameters which remain $`\mathrm{\Omega }_B=0.135`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=1.0`$, $`n_S=0.94`$, $`h=0.45`$, but the model now requires $`C_{10}=0.8`$ and an upward shift in the calibration for BOOMERanG-98 of $`30\%`$. This model is ruled out by the COBE and BOOMERanG-98 data at the 95% C.L., although this discrepancy can be resolved by adding a gravitational wave (GW) contribution on large scales. The top panel of Figure 2 indicates that the inclusion of a GW component such that $`C_{10}^T=C_{10}^S`$ leaves the small scale behaviour practically unchanged but provides enough CMB power on large scales to match the COBE normalization.

A GW contribution as large as $`C_2^T/C_2^S1`$ is incompatible at the $`95\%`$ C.L. with flat $`\mathrm{\Lambda }`$CDM models when the BBN prior is assumed in CMB analyses (Kinney, Melchiorri & Riotto (2000)). GWs leave a characteristic imprint on polarization power spectra, via fluctuations in the $`B`$ magnetic-type-parities channel (Kamionkowski, Kosowsky & Stebbins (1996), Seljak & Zaldarriaga (1996)) that vanish in the case of scalar fluctuations. However, the expected polarization amplitude for BDM models is very low, below $`0.5\mu K`$, and future experiments will not have the sensitivity to detect the $`B`$ channel signal unless a substantial amount of reionization occurred in the past (Seljak & Zaldarriaga (1996)).

FIG. 2.— In the top panel, the BOOMERanG-98 and the COBE data are plotted together with the best-fitting baryonic model to the BOOMERanG-98 data. In order to fit the COBE data, a tensor component must be added. For comparison, the $`\mathrm{\Lambda }`$-CDM concordance model is also plotted. The middle panel shows the expected polarization signal of the best-fitting model. Including the tensor components produces $`B`$ modes in the polarization spectra. In the bottom panel it can be seen that the best-fitting model to the BOOMERanG-98 data is in contradiction with the matter power spectrum measured by the PSCz survey.

One can further combine the CMB constraints with those obtained from the luminosity–distance measurements of high-$`z`$ supernovae (Perlmutter et al. (1998); Schmidt et al. (1998)) as we do in the bottom panel of Figure 1. The inclusion of the SNIa data does not greatly affect the constraints on $`\mathrm{\Omega }_\mathrm{\Lambda }`$ but, as expected, it rules out any models with low values of $`\mathrm{\Omega }_b`$. We find that $`\mathrm{\Omega }_b0.13`$ at the $`95\%`$ C.L. and an increasing inconsistency of BDM models with a flat universe.

Also plotted in the bottom panel of Figure 1 are the constraints obtained including a prior $`\mathrm{\Omega }_bh^2=0.019\pm 0.002`$ (Burles, Nollett & Turner (2000)) inferred by measurements of primordial elements and assuming standard BBN. The results of the $`2`$ priors are incompatible at more than $`3\sigma `$.

It is interesting to further explore the compatibility of our results with independent observations by considering likelihood contours in the ($`\mathrm{\Omega }_b`$, $`h`$) plane. Figure 3 plots the results of our analysis with no priors, the results including the SNIa data and the results

FIG. 3.— The likelihood contours in the ($`\mathrm{\Omega }_B`$, $`h`$), with the remaining parameters taking their best-fitting values to the BOOMERanG-98 data. The contours correspond to 0.32, 0.05 and 0.01 of the peak value of the likelihood, which are the 68%, 95% and 99% C.L. respectively. The filled contours correspond to an analysis with a $`t>10`$ Gyr prior. The dashed, dotted and solid contours to the right of the plot are the result when the SNIa data is included. The bold black line to the left of the plot is the $`95\%`$ C.L. region from BBN. with the BBN constraint. We can see from this figure that purely baryonic models are in quite good agreement with BOOMERanG-98 and with the BBN constraints. However, when the SNIa data is included we are restricted to models with values of $`\mathrm{\Omega }_b`$ that are too high and values of $`h`$ that are too low to be consistent at $`<3\sigma `$ with both the BBN constraint and the recent HST result $`h=0.72\pm 0.08`$ (Freedman et al. (2000)).

Quite independently of the prior assumed, we find that the age of the universe in this scenario is constrained to be $`t_0=24\pm 2`$ Gyr, nearly double that expected from standard CDM. The age of the universe must exceed the ages of the oldest globular clusters $`t_{GC}=14\pm 2`$ Gyr by an amount $`\mathrm{\Delta }t0.52`$ Gyr, so the age of the universe required by these models seems to be too high unless globular cluster formation is delayed for a very long time.

Collisional damping erases fluctuations in the matter power spectrum on scales $`<100h^1Mpc`$. The predicted matter power spectrum from the best-fitting model (see Figure 2, bottom panel) is in contradiction with the decorrelated linear power spectrum extracted from the recent PSCz catalogue (Hamilton & Tegmark (2000)). Allowing for a shift in the overall amplitude with a bias factor $`b`$ where $`P(k)_{PSCz}=b^2P(k)`$, we obtain a best-fitting bias $`b=4.5`$ with a $`\chi ^2=155`$ which is an extremely poor fit to 22 data points with one free parameter. The theoretical variance in baryonic matter fluctuations over a sphere of size $`8h^1\mathrm{Mpc}`$, for example, is $`\sigma _810^2`$ to be compared with the observational value $`\sigma _8=0.56\mathrm{\Omega }_m^{0.47}1.4`$ (Viana & Liddle (1999)). So, rescaling the galaxy data with a massive bias parameter may not be physically meaningful since no galaxies would have formed in order to allow there to be a predicted galaxy power spectrum.

Various phenomenological mechanisms can be proposed in an attempt to resolve the discrepancy with the matter power spectrum and render BDM models viable. One such possibility is to increase the spectral index of the primordial density fluctuations, $`n_s`$. However, the ’blue’ models that are able to produce $`\sigma _80.5`$ require $`n_s2`$ which is outside the range $`0.7\stackrel{<}{}n_s\stackrel{>}{}1.2`$ allowed by inflation (Steinhardt (1995)). Also, such a high tilt will cause the CMB spectrum to become incompatible with observations.

Another option that could increase small scale density perturbations whilst leaving the CMB spectrum unchanged is to have a spectral index that varies with scale, $`\mathrm{d}n_s/\mathrm{d}\mathrm{ln}k0`$, as predicted by some inflationary models (Kosowsky & Turner (1995)). Unfortunately, the variation needed, $`n_s\stackrel{>}{}1`$, is much greater than that predicted by the most viable models.

More recently, there has been some investigation into features in the primordial power spectrum generated during inflation (Griffiths, Silk & Zaroubi (2000); bar; Knebe, Islam & Silk (2001)). It is possible that such features could increase the perturbations on small scales without affecting the CMB spectrum, although it is doubtful that a single feature would be able to resolve the observed discrepancy and such a mechanism is unlikely to be preferable to the assumption of CDM itself.

## 4. Discussion

We have examined CMB anisotropies and large-scale structure observations in a purely baryonic dark matter universe. Our results suggest that a BDM adiabatic universe with a power-law primordial power spectrum can only reproduce the sub-degree CMB measurements if the universe is closed. Although the premise that $`\mathrm{\Omega }=1`$ is considered to be one of the basic predictions of the inflationary scenario, it is possible to construct inflationary scenarios that predict closed universes (see e.g. Linde (1995); Starobinsky (1996)). It is interesting to note, however, that removing CDM immediately forces us to construct a more elaborate inflationary model, or, in other words, the most simple model needs non-baryonic matter to work.

When flatness was assumed, we found that $`\omega _{cdm}>0.04`$ at $`2\sigma `$ with $`L(\omega _{cdm}=0)/L(\omega _{cdm}=0.1)10^3`$. Therefore, the flatness constraint with BBN renders $`\mathrm{\Lambda }BDM`$ models less consistent with the CMB observations than the standard $`\mathrm{\Lambda }CDM`$ model. Also, $`\mathrm{\Lambda }BDM`$ models that fit the BBN constraints are degenerate with a non-physical region of the parameter space (with very high baryon content and small $`h`$). When any cosmological prior other than BBN is assumed in the analysis, the unphysical region turns out to be preferred. Including the SNIa constraints (or $`\mathrm{\Omega }_m=\mathrm{\Omega }_b>0.25`$), for example, makes the parameter space incompatible at $`>95\%`$ confidence with both the HST constraint $`h=0.72\pm 0.08`$ and with BBN.

As expected, all the BDM adiabatic models fail to reproduce both the CMB observations and the observed amount of galaxy clustering on $`8h^1\mathrm{Mpc}`$ scales. In order to solve all the discrepancies, that are not present in the standard $`\mathrm{\Lambda }`$CDM scenario, we need to introduce ad hoc mechanisms which are unlikely to be preferable to the assumption of CDM itself. Nonetheless, BDM models leave a set of characteristic imprints such as large scale $`B`$-mode polarization and no third acoustic peak in the anisotropy spectrum that will allow future experiments to further scrutinize this hypothesis.

We thank Pedro Ferreira, Arthur Kosowsky, Andrew Liddle, Francesco Melchiorri, Gary Steigman for useful conversations and the High-Z Supernova Search Team for providing the SNIa likelihoods.