text:
# A Hagedorn temperature in post-inflationary dynamics?

## Abstract

We consider the possibility of regarding the maximum temperature $`T_{\text{max}}`$ of a plasma of inflaton decay products as a Hagedorn temperature connected with the transition from an inflaton-dominated epoch to a radiation-dominated universe. Employing a “bottom-up” approach to inflaton decay, we show that the limiting temperature $`T_\mathrm{H}m_\varphi `$, as obtained from a numerical analysis of the evolution equations, can also be understood theoretically in a Hagedorn picture of the phase transition.

The concept of Hagedorn temperature hagedorn65  in statistical systems is commonly invoked as a scale that signals the onset of a phase transition Atick88 . While it does not directly relate to the usual thermodynamic definition of temperature, it can be regarded as a limiting temperature where the population of energy states begins to exceed the Boltzmann weight, implying that the original degrees of freedom strongly overlap in momentum space and lose their identity. One example of this is the melting of hadronic and glueball states in relativistic heavy-ion collisions, where the Hagedorn temperature $`T_\text{H}`$ can be derived from the Bag model Chodos74  or a stringy model Huang70  of hadrons, and can be used as an estimate of the deconfinement temperature, though the exact equivalence of the two is hard to establish, at least in four-dimensional gauge theories Sundborg . However, the idea of a number density of states that grows exponentially with temperature is quite natural coming from a theory with confining forces (like Regge theory), and may be useful when the exact nature of these microscopic forces is far from clear. A closely related example is Hagedorn-type behavior in some string theories at weak coupling Atick88 . In the context of early string dominated cosmology, a Hagedorn-type phase could be associated with open strings attached to D-branes which has recently been shown to provide a negative pressure necessary to drive inflation at high temperatures near the string scale Kogan . This is accompanied by an exponential growth in the number of string modes with energy.

In this letter, we explore the possibility to have a Hagedorn limit of a quite different nature in cosmology, namely, in post-inflationary dynamics that culminate in reheating of the universe. Inflation is a well known contender for making the Universe flat, besides generating adiabatic density perturbations with an almost flat spectrum Liddle . However, inflation must end in order to provide a thermal bath of radiation in order to facilitate synthesis of light nuclei. The standard paradigm is that the inflaton field, $`\varphi `$, enters the oscillatory phase and decays perturbatively into radiation quanta, which subsequently thermalize completely through interactions, yielding a reheat temperature $`T_{\text{reh}}`$, when the universe is completely radiation dominated Turner . It has also been proposed that an initial non-perturbative decay of the inflaton, dubbed as preheating, is a plausible scenario Robert . Nevertheless there has to be a subsequent stage of thermalization after preheating quite similar to perturbative reheating Sarkar . Yet another interesting possibility is to have a non-topological defect formation from the fragmentation of the inflaton condensate, such as $`Q`$-balls Mazumdar , where reheating occurs via surface evaporation of $`Q`$-balls.

Irrespective of a particular scenario the final reheating temperature plays a very crucial role in the early universe, for e.g. see Kari . The reheat temperature can be constrained by bounds coming from thermal and non-thermal gravitino production during reheating/preheating in supersymmetric theories Nanopoulos . Typically, $`T_{\text{reh}}<10^9`$ GeV is required, for gravitino mass $`m_{3/2}100`$ GeV, in order to avoid constraints coming from successful predictions of Big Bang Nucleosynthesis (BBN), for e.g. see Sarkarrep .

The reheat temperature is not, however, the maximum temperature of a thermal plasma, which may be two orders of magnitude higher (but less than $`m_\varphi `$Turner ; Riotto . Physically, this maximum temperature, denoted by $`T_{\text{max}}`$, arises due to the competition between early entropy production within light relativistic degrees of freedom at the earliest stages of inflaton decay and continued expansion of the universe which dilutes the energy density. The main objective of this work is to show that $`T_{\text{max}}`$ may be regarded as a limiting temperature, if we trace the inflaton decay process backwards in time, up to a point when coherent inflaton production from the radiation fields begins. We take into account thermal factors that capture the essential feature of Bose enhancement when the inflaton population builds up. We support our theoretical conjecture about this limiting temperature as a Hagedorn temperature with a numerical analysis of the evolution equations for matter (inflaton) and radiation energy densities run backwards in time.

After inflation ends, the energy density stored in coherent oscillations of the inflaton field is transferred to radiation quanta over a period of time. Several features of this transition, such as thermalization Sarkar , associated heavy-particle production Riotto , and non-perturbative processes (preheating) have been widely investigated Robert ; Mazumdar .

However, little is known about the dynamical process that results in inflaton decay. We may invert this problem, and ask the question as to how the inflaton field may be generated coherently from the radiation fields? Although the coupling between the two is small, our view is that Bose enhancement of the inflaton in the final state can lead to coherence, similar to a condensate effect. Incorporating this feature in the Boltzmann equations that describe the perturbative decay of the inflaton is not hard. We may run the system of equations (see below) backwards in time (decreasing scale factor) to study how the energy densities evolve. We show that an exponential enhancement of the inflaton density can overcome the Boltzmann weight, very similar to the Hagedorn picture. The temperature at this instant is the limiting temperature, beyond which the radiation energy density drops rapidly to zero. This may be interpreted as $`T_{\text{max}}`$, the maximum or limiting temperature. Thus, we have an alternate physical description of the maximum temperature, and we will quantify it both theoretically and numerically in following.

We consider the interaction $`L_{int}=g\varphi \overline{f}f`$ where $`f`$ denotes a fermion 111For Standard model quarks or leptons, such an interaction cannot be written down since it would not be a $`SU(2)`$ gauge invariant monomial. For the sake of exemplification, we assume that these fermions are pure singlet.. The evolution equation for the inflaton energy density $`\rho _\varphi `$ is given by

$`\dot{\rho _\varphi }+3H\rho _\varphi `$ $`=`$ $`\mathrm{\Gamma }_\varphi \rho _\varphi (1+n_B(m_\varphi )),`$ (1)
$`\mathrm{\Gamma }_\varphi `$ $`=`$ $`{\displaystyle \frac{g^2m_\varphi }{32\pi }}=\alpha _\varphi m_\varphi ,`$ (2)
where $`H`$ is the Hubble expansion rate and $`\mathrm{\Gamma }_\varphi `$ is the inflaton decay rate. Note that the source term for $`\rho _\varphi `$ is enhanced by the Bose factor $`(1+n_B(m_\varphi ))=1/(1\mathrm{e}^{m_\varphi /T})`$, since the inflaton is a scalar field. Since the expression for $`\mathrm{\Gamma }_\varphi `$ is derived for zero total momentum (center-of-mass frame), the $`p=0`$ mode of the inflaton is being populated. $`\rho _\varphi `$ as defined here is the thermal energy density

$$\rho _\varphi =\frac{m_\varphi ^4}{(\mathrm{e}^{m_\varphi /T}1)}.$$
(3)
The thermal factor comes out naturally when one considers fermions described by a thermal distribution fusing to form the inflaton. The analogous equation for $`\rho _\text{R}`$, the energy density of radiation is

$$\dot{\rho _\text{R}}+4H\rho _\text{R}=\mathrm{\Gamma }_\varphi \rho _\varphi (1+n_B(m_\varphi )).$$
(4)
We evolve eqns. (1) and (4) numerically, for which purpose it is convenient to define rescaled quantities as

$$x=a/m_\varphi ,\mathrm{\Phi }=\rho _\varphi a^3,R=\rho _\text{R}a^4/m_\varphi .$$
(5)
where $`a`$ is the scale factor of the Friedman-Robertson-Walker metric. In terms of rescaled quantities, the evolution equations take the simple form

$`\left({\displaystyle \frac{d\mathrm{\Phi }}{dx}}\right)`$ $`=`$ $`c{\displaystyle \frac{x\mathrm{\Phi }}{\sqrt{x\mathrm{\Phi }+R}}}(1+n_B(x,R)),`$ (6)
$`\left({\displaystyle \frac{dR}{dx}}\right)`$ $`=`$ $`c{\displaystyle \frac{x^2\mathrm{\Phi }}{\sqrt{x\mathrm{\Phi }+R}}}(1+n_B(x,R)).`$ (7)
where $`c=\sqrt{\frac{3}{8\pi }}\alpha _\varphi M_{Pl}m_\varphi ^{5/2}`$. These equations are similar in form to those in Riotto , but the rescaling is different since we wish to study the build-up rather than the decay of the inflaton, for which time runs backwards. Furthermore, we include the Bose enhancement factor that was not required in Riotto . For the numerical analysis, we choose $`m_\varphi =10^{13}`$ GeV, $`T_{\text{reh}}10^9`$ GeV. The traditional expression for the reheat temperature corresponding to a thermal bath dominated solely by relativistic particles is given by Turner

$$T_{\text{reh}}=0.2\left(\frac{200}{g_{}}\right)^{1/4}\sqrt{\alpha _\varphi m_\varphi M_{Pl}},$$
(8)
where $`g_{}`$ accounts for the relativistic degrees of freedom in the plasma. Using this relation, and with $`g_{}=200`$, we obtain $`\alpha _\varphi =2.3\times 10^{13}`$ and $`c=8.7\times 10^{37.5}`$. We begin the evolution from $`T_{\text{reh}}10^9`$ GeV, which fixes the radiation energy density through the Stefan-Boltzmann law

$$\rho _\text{R}=\left(\frac{g_{}\pi ^2}{30}\right)^{1/4}T^42.85\times 10^{36}\mathrm{GeV}^4$$
(9)
However note that in our case there is no inflaton oscillation dominated phase and radiation is the only component we start with. In this respect the traditional notion of reheat temperature, which applies to a largest temperature of a thermal bath only dominated by relativistic particles, does not hold in this scenario. Nevertheless the numerical significance of eqn. (8) continues to be important as an estimation of a temperature of a relativistic bath when the entire inflaton quanta has been converted into radiation.

The qualitative features of the result of the numerical analysis are independent of the choice of initial value of scale factor $`a`$, although $`x,\mathrm{\Phi },R`$ are quantitatively dependent on it. We choose $`a=10^8(x=10^5)`$ at the beginning, and the results of the evolution across nearly 4 e-foldings are displayed in Fig. 1.

While the radiation curve $`R(x)`$ (long dashed line) is remaining constant, the temperature $`Ta^1`$, as expected for the radiation-dominated epoch. This is because both $`\rho _\text{R}`$ and $`a`$ are changing. As the inflaton energy density (solid line) builds up, the radiation energy density drops rapidly to zero. At this point, it is easy to see from eqn.(6) that $`\mathrm{\Phi }(x)x^3`$, which is the usual initial condition obtained in standard evolution of the inflaton decay equation. Physically, our picture is quite different from that in Riotto , in that we can interpret $`T_{\text{max}}`$ as a limiting temperature, where the inflaton fields acquires a coherent nature, and the phase transition to an inflaton-dominated universe is underway. This interpretation, though novel, can be supported by a simple analytical estimate of $`T_{\text{max}}`$ which we presently derive.

It is also interesting to observe that in our case the process of thermalization is actually occurring very close to $`T_{\text{max}}`$ in such a way that entropy density of a radiation dominated bath follows an adiabatic expansion law $`sa^3=\mathrm{constant}`$, where $`s`$ is the entropy density. In this respect we are tempted to associate $`T_{\text{max}}`$ as an actual reheat temperature of the universe.

At the limiting temperature, when the inflaton energy density is dominant, $`\mathrm{\Gamma }_\varphi H`$ holds, and we may neglect the $`3H\dot{\rho _\varphi }`$ term in eqn.(1). Assuming the temperature to be approximately constant over a small range of time, we may integrate the equation from time $`t_0`$ to $`t`$ ($`t_0<t`$ in the absolute sense) to obtain

$$\mathrm{ln}\left(\frac{n_\varphi (t)}{n_\varphi (t_0)}\right)=\frac{\alpha _\varphi m_\varphi (tt_0)}{1\mathrm{e}^{m_\varphi /T}},$$
(10)
where $`n_\varphi =\rho _\varphi /m_\varphi `$. This equation governs the rate of increase in the inflaton number density in a small interval $`tt_0`$ near the transition. This time interval must be of the order of (or slightly less than) $`1/\mathrm{\Gamma }_\varphi 1/(\alpha _\varphi m_\varphi )`$, because subsequently, the inflaton density drops off much more rapidly than the temperature. This means we may approximate as

$$n_\varphi (t)n_\varphi (t_0)\mathrm{e}^{1/(1\mathrm{e}^{m_\varphi /T})}.$$
(11)
Now, since $`n_\varphi (t_0)`$ includes a statistical weight factor, $`\mathrm{e}^{m_\varphi /T}`$, (we assume $`Tm_\varphi `$ so that it is approximately a Boltzmann weight), the growth in number density of inflatons begins to overwhelm the exponential suppression from the statistical weight when

$$Tm_\varphi (1\mathrm{e}^{m_\varphi /T}),$$
(12)
which is uniquely satisfied for $`T0.74m_\varphi `$. We regard this as a maximum or limiting temperature, reminiscent of the Hagedorn temperature from hadronic physics. This is the main result of this work. Beyond this temperature, the radiation fields overlap strongly and populate the zero momentum mode of the inflaton at a rate that is sufficient to overcome the thermal suppression. The universe then has only an inflaton component to its energy density. We see that the analytical estimate of $`T_{\mathrm{max}}`$ or $`T_\text{H}`$ is close to the numerical value obtained from evolving the Boltzmann equations for inflaton decay backwards in time. This numerical estimate is robust and is not affected by inclusion of inverse reaction terms in the collision kernel of the Boltzmann equations, which only become significant at temperatures exceeding the limiting temperature. In Kolb:2003ke , modulo extra assumptions on the gauge interactions of the fermions and their coupling to the gauge fields, it was argued that $`T_{\mathrm{max}}`$ could exceed $`m_\varphi `$. We find that in the Hagedorn picture, including thermal masses for the fermions, as in Kolb:2003ke , does not affect the bound $`T_{\mathrm{max}}0.74m_\varphi `$.

In summary, we have shown that it is possible to interpret the maximum temperature reached in the evolution of inflaton decay as a Hagedorn temperature. Viewed backwards in time, or with decreasing size of the universe, at this temperature, which is close to $`m_\varphi `$, coherent population of the zero momentum state of the inflaton takes place. The rate of population exceeds the thermal suppression factor, and it is not appropriate to use the evolution equations any longer, since a clear distinction between radiation and inflaton fields cannot be made. This is similar to the purported hadron gas to quark-gluon plasma (QGP) transition at finite temperature $`Tm_\pi `$, where the Hagedorn temperature signifies when the quark fields start to percolate through the overlapping hadrons. At this point, hadrons are no longer genuine physical degrees of freedom, and the deconfinement phase transition is approached. We have applied the same idea to the early stage of inflaton decay, and supported the conjecture with a numerical analysis of the Boltzmann equations.

We find that the maximum temperature, or the corresponding Hagedorn temperature of the post-inflationary universe is very close to the inflaton mass, $`T_{\text{max}}0.74m_\varphi `$. We regard this is as an actual reheat temperature of the universe because below the phase transition, the subsequent evolution of radiation density obeys an adiabatic expansion law with $`Ta^1`$. We consider this result as having an impact on particle-cosmology which we shall explore in a separate publication.

We thank Jim Cline for important remarks. P. J. is supported in part by the Natural Sciences and Engineering Research Council of Canada and in part by the Fonds Nature et Technologies of Québec. A. M. is a CITA national fellow.


format: text

text:
# Nature of excited 0+ states in 158Gd described by the projected shell model

## Abstract

Excited 0+ states are studied in the framework of the projected shell model, aiming at understanding the nature of these states in deformed nuclei in general, and the recently observed 13 excited 0+ states in 158Gd in particular. The model, which contains projected two- and four-quasiparticle states as building blocks in the basis, is able to reproduce reasonably well the energies for all the observed 0+ states. The obtained B(E2) values however tend to suggest that these 0+ states might have a mixed nature of quasiparticle excitations coupled to collective vibrations.

Dynamic perturbations of nuclear shapes around the equilibrium can give rise to physical states at low to moderate excitation energies. Classical examples of such motion are $`\beta `$ and $`\gamma `$ vibrations BM75 ; Ia84 , in which nucleons undergo vibrations in a collective manner. Traditionally, the first excited $`K^\pi =0^+`$ states and the first excited $`2^+`$ states are interpreted, respectively, as the $`\beta `$ and $`\gamma `$ vibrational states. While the $`2^+`$ collective excitations are better understood theoretically, the nature of the lowest $`0^+`$ excitation of deformed nuclei still remains under debate CB94 ; Burke ; Kumar ; Bonn ; Garrett . The physics of higher $`0^+`$ states is even more complex because, on one side, they can predominantly be multi-phonon states based on the single-phonons Ani02 , and on the other side, they can be quasiparticle (qp) excitations in nature. The real situation is, perhaps, that the two aspects, collective excitations and qp states, are mixed by residual interactions.

Data on $`K^\pi =0^+`$ states have been relatively sparse. In a very recent work by Lesher et. al. Exp , a remarkable $`(p,t)`$ experiment revealed a total of 13 excited 0+ states in 158Gd, below an excitation energy of approximately 3.1 MeV. This abundance of 0+ states in a single nucleus provides significant new information on the poorly understood phenomenon, which has immediately sparked off theoretical interest. For this energy range, one may think about an explanation through collective modes. In fact, Zamfir, Zhang, and Casten GCM  suggested that many of the observed 0+ states may be of two-phonon octupole character. Nevertheless, these authors warned also that, although the mechanism was excluded in their collective models, many of the 0+ states in this excitation energy range may be predominantly two-qp in character.

The purpose of the present paper is to investigate whether one can explain the nature of the observed $`0^+`$ states in terms of qp excitations. In contrast to Ref. GCM , our calculation here does not emphasize the aspect of collective excitation. Although, similar to the work of Zamfir, Zhang, and Casten GCM , we may provide only a partial image, our results may shed light on the importance of considering the qp aspects in understanding the nature of these $`0^+`$ states.

Our study is based on the projected shell model (PSM) review . The PSM is the spherical shell model built on a deformed basis. The PSM calculation usually begins with the deformed Nilsson single-particle states at a deformation $`\epsilon `$. Pairing correlations are incorporated into the Nilsson states by BCS calculations. The consequences of the Nilsson-BCS calculations provide us with a set of qp states that define the qp vacuum $`|\varphi (\epsilon )`$. One then constructs the shell model bases by building multi-qp states. The broken symmetry in these states is recovered by angular momentum projection review  (and particle number projection, if necessary) to form a shell model basis in the laboratory frame. Finally a two-body shell model Hamiltonian is diagonalized in the projected space.

To determine the deformation at which the shell model basis should be built, we first study the bulk properties of 158Gd including the deformation. In a microscopic calculation, one searches for energy minima by varying the deformation parameter. Here, we calculate the angular-momentum-projected energies having the form

$$E^I(\epsilon )=\frac{\varphi (\epsilon )\left|\widehat{H}\widehat{P}^I\right|\varphi (\epsilon )}{\varphi (\epsilon )\left|\widehat{P}^I\right|\varphi (\epsilon )},$$
(1)
where $`P^I`$ is the angular-momentum projection operator RS80  which projects the mean-field vacuum $`|\varphi (\epsilon )`$ onto states with good angular momentum. As many previous PSM calculations for the rare earth nuclei, particles in three major shells ($`N=4,5,6`$ for neutrons and $`N=3,4,5`$ for protons) are activated in the present calculation for 158Gd. For comparison, unprojected energies

$$E(\epsilon )=\frac{\varphi (\epsilon )\left|\widehat{H}\right|\varphi (\epsilon )}{\varphi (\epsilon )|\varphi (\epsilon )}.$$
(2)
are also calculated.

As one can see in Fig. 1, the lowest energy for a given angular momentum $`I`$ is well localized at deformations varying from $`\epsilon 0.24`$ at $`I=0`$ to $`\epsilon 0.28`$ at $`I=12`$. The ground-state deformation calculated by Möller et. al. MN95  for this nucleus yields $`\epsilon =0.25`$. These angular-momentum-projected minima lie at deformations that are slightly larger than the mean-field minimum ($`\epsilon 0.23`$). Fig. 1 indicates that 158Gd is a stably deformed nucleus against rotation, with a pronounced energy minimum corresponding to a deformed prolate shape. In the following calculations, we thus construct the shell model basis at the deformation $`\epsilon =0.26`$. The vacuum state $`|\varphi (\epsilon =0.26)`$ is hereafter written as $`|0`$.

Following the spirit of the Tamm-Dancoff method RS80 , we build the shell model space by including 0-, 2- and 4-qp states:

$$|\mathrm{\Phi }_\kappa =\{|0,\alpha _{n_i}^{}\alpha _{n_j}^{}|0,\alpha _{p_i}^{}\alpha _{p_j}^{}|0,\alpha _{n_i}^{}\alpha _{n_j}^{}\alpha _{p_i}^{}\alpha _{p_j}^{}|0\},$$
(3)
where $`\alpha ^{}`$ is the creation operator for a qp and the index $`n`$ ($`p`$) denotes neutron (proton) Nilsson quantum numbers which run over the low-lying orbitals. Thus, the projected multi-qp states are the building blocks of our shell model basis:

$$|\mathrm{\Psi }_M^I=\underset{\kappa }{}f_\kappa ^I\widehat{P}_{MK}^I|\mathrm{\Phi }_\kappa .$$
(4)
Here, $`\kappa `$ labels the basis states and $`f_\kappa ^I`$ are determined by configuration mixing.

We then diagonalize the Hamiltonian in the projected multi-qp states given in (4). In the calculation, we employ a quadrupole plus pairing Hamiltonian, with inclusion of quadrupole-pairing term review

$$\widehat{H}=\widehat{H}_0\frac{1}{2}\chi \underset{\mu }{}\widehat{Q}_\mu ^{}\widehat{Q}_\mu G_M\widehat{P}^{}\widehat{P}G_Q\underset{\mu }{}\widehat{P}_\mu ^{}\widehat{P}_\mu ,$$
(5)
where $`\widehat{H}_0`$ is the spherical single-particle Hamiltonian which contains a proper spin-orbit force. The quadrupole-quadrupole interaction strength $`\chi `$ is determined by the self-consistent relation with deformation $`\epsilon `$. The monopole pairing strength $`G_M`$ is taken to be $`G_M=\left[2013(NZ)/A\right]/A`$ with “$``$” sign for neutrons and “$`+`$” sign for protons. Finally, the quadrupole pairing strength $`G_Q`$ is assumed to be proportional to $`G_M`$, the proportionality constant being fixed to 0.20 in the present work. These interaction strengths are the same as the values used in the previous PSM calculations for the rare earth nuclei review .

In Fig. 2, we show the partial theoretical spectrum in 158Gd up to 4 MeV in energy and $`18\mathrm{}`$ in spin. We emphasize that all these states have been obtained by a single diagonalization, without any adjustment for individual states. Out of these many states, let us concentrate on the lowest one at each spin (the yrast band, denoted by diamonds), and on all the $`0^+`$ states (denoted by filled tri-angles). Although it is not a focus of our discussion in this paper, we believe that a comparison of the yrast band with known data can provide a strict test of the model and can provide a useful constraint to the calculations of the $`0^+`$ states.

Fig. 3 presents the PSM results for the yrast band in 158Gd (the same values shown in Fig. 2 as diamonds), which are compared with the known data isotopes , in a plot of $`\gamma `$-ray energy versus spin. As can be seen, the data are described very well. The calculations predict a sudden drop in the curve at spin $`I=16`$, corresponding to a backbending in the moment of inertia. This sudden change occurs just at the upper part of the band where the current measurement stops. Extension from the current measurement should see this phenomenon and will provide a strict test of our model prediction. In this regard, we note that in the isotonic chain of nuclei, 160Dy, 162Er, 164Yb, and 166Hf, similar backbending effects have been observed and successfully described elsewhere by the PSM backb1 ; backb2 .

We now turn our discussion to the excited $`0^+`$ states. Lesher et. al. Exp  observed 13 $`0^+`$ states in the nucleus 158Gd, within an excitation energy range of 1.2 MeV to 3.1 MeV. In Fig. 4, we plot all the theoretical $`0^+`$ states (up to 3.2 MeV in energy; the same values shown in Fig. 2 as filled triangles) in the order of excitation energies. The experimental data Exp  are shown in the same plot for comparison. The predicted $`0^+`$ states are found to be in the right energy range, although deviations can be clearly seen between theory and experiment.

What we found impressive is the number of $`0^+`$ states predicted by the calculation. The PSM produces a sufficient number of $`0^+`$ states to be compared with data. As described in Eq. (4), the total wave function $`|\mathrm{\Psi }_M^I`$ is a linear combination of the (projected) basis states given in (3). The basis states in (3) are not arbitrarily selected but are taken from all the neutron and proton Nilsson orbitals that lie close to the Fermi surface. In 158Gd, the relevant orbitals are: $`\frac{3}{2}[521]\nu `$, $`\frac{5}{2}[523]\nu `$, $`\frac{11}{2}[505]\nu `$, $`\frac{3}{2}[651]\nu `$, and $`\frac{5}{2}[642]\nu `$ for neutrons, and $`\frac{1}{2}[411]\pi `$, $`\frac{3}{2}[411]\pi `$, $`\frac{5}{2}[413]\pi `$, $`\frac{5}{2}[532]\pi `$, and $`\frac{7}{2}[523]\pi `$ for protons. In each of these 10 near-Fermi orbitals, nucleons having opposite signs for the $`K`$ quantum number can couple to a 2-qp state with total $`K=0`$. Combination of a pair of 2-qp states can further give $`K=0`$ 4-qp states. If one neglects the coupling of these qp states to the collective states, the number of low-lying $`0^+`$ states in deformed nuclei obtained in this way depend solely on the single-particle level density and the level distribution near the Fermi surface. Since similar conditions can also be found in many other rare earth nuclei, we expect such an abundance of $`0^+`$ states as found in 158Gd not to be an isolated case. We predict that such an abundance of $`0^+`$ states should also be observed in many other nuclei.

The large number of $`0^+`$ states is difficult to obtain within collective models. In Ref. GCM , one could obtain at most 5 excited $`0^+`$ states up to 3.2 MeV in calculations with the geometric collective model or the interacting boson model (with only $`s`$\- and $`d`$-boson). The reason is that within such collective models, the number of degrees of freedom of collective motion is limited. Only if one considered the odd-parity bosons in an extended boson space, could the authors in Ref. GCM  get more excited $`0^+`$ states.

In Table I, we list the 18 calculated $`0^+`$ states below an excitation energy of 3.25 MeV (slightly higher than the highest experimental $`0^+`$ state). Their leading configurations are also given. Note that each of the configurations listed here can not be pure, but is the dominant part in the corresponding wave function. One sees that most 2-qp states coming from the near-Fermi Nilsson levels have lower energies. Due to the varying responses of the residual interactions through configuration mixing, one finds that some 4-qp states are lower in excitation energy than 2-qp states.

Up to now, we have seen that the PSM calculations, which explicitly include two- and four-qp states built on the basis but no vibrational degrees of freedom, can reasonably produce a sufficient number of $`0^+`$ states within the right excitation energy range. This fact suggests that these states are qp states in nature contrary to the work of Zamfir, Zhang, and Casten GCM . However, before we do that, we should study the transition properties of these states.

The excited $`0^+`$ states can decay to the $`2_g^+`$ state in the ground state band through E2 transition. One such a transition was measured in Ref. Borner . Preliminary results of Lesher et. al. BE2  show 12 such transitions. The PSM calculations of the transition from the $`i`$-th $`0^+`$ state to the $`2_g^+`$ state, B(E2, $`0_i^+2_g^+`$), are listed in Table I. Comparing the theoretical B(E2) values with the experimental data BE2 , we found that the numbers (in W.u.) listed in Table I are one or two orders of magnitude too small for many of the transitions.

We note that the B(E2, $`0_i^+2_g^+`$) values (a few W.u.) of Lesher et. al. BE2  are in average of two orders of magnitude smaller than the in-band B(E2) values of the ground-state band (typically, a few hundred W.u.). These small B(E2) values usually indicate that the observed $`0^+`$ states in Ref. Ani02  should carry significant quasiparticle components. However, our calculated B(E2) results are much smaller than the values of Lesher et. al., which might indicate an insufficient mixing of collectivity to the qp states in our model. This tends to suggest that although most, if not all, of these $`0^+`$ states have significant qp components, they are also mixed with correlations induced by the collective motion QPNM . These correlations can be introduced in the PSM framework by inclusion of interactions of higher orders of the multipole type CG01 , and by addition of the $`D`$-pair operators to the vacuum state Heavy , which takes both quasiparticle and collective degrees of freedom explicitly into account in a shell model basis. Generator Coordinate Method, which consists of a construction of a linear superposition of different product wave functions, can also be adopted by the PSM.

In summary, the projected shell model was employed to understand the nature of excited $`0^+`$ states in deformed nuclei. The shell model space consists of projected 2- and 4-qp states on top of the deformed BCS vacuum state. Therefore, the calculation emphasized the quasiparticle character of these states, in contrast to previous work GCM  using collective models. Our results were compared with the remarkable example of recently observed 13 excited $`0^+`$ states in 158Gd Exp . After performing exact angular momentum projection and configuration mixing calculations (by two-body residual interactions) for all the possible low-lying qp configurations that can give rise to $`K^\pi =0^+`$ states, we found that the obtained energy levels as well as the number of the states can reasonably explain the data. Preliminary measurements of B(E2) values suggest that mixing of these qp states with the collective, vibrational motion may not be neglected.

Research at University of Notre Dame is supported by the NSF through contract No. PHY-0140324. Research at the University of Tennessee is supported by the U. S. Department of Energy through Contract No. DE–FG05–96ER40983.


text:
# The 𝑞𝑞̄ spectra and the structure of the scalar mesons

## 1 Introduction

Since Gell-Man conjecture, most of the meson experimental data were classified as $`q\overline{q}`$ states according to $`SU(N)`$ irreducible representations. Nevertheless a number of interesting issues remains still open as for example the understanding of the structure of the scalar mesons or the new $`D_s`$ states measured on B factories.

The theoretical study of charmonium and bottomonium made clear that heavy-quark systems are properly described by nonrelativistic potential models reflecting the dynamics expected from QCD eich . The light meson sector has been studied by means of constituent quark models, where quarks are dressed with a phenomenological mass and bound in a nonrelativistic potential, usually a harmonic oscillator isgu . Quite surprisingly a large number of properties of hadrons could be reproduced in this way godf . In this talk we present the meson spectra obtained by means of a chiral constituent quark model in a trial to interpret some of the still unclear experimental data in the scalar sector.

## 2 SU(3) chiral constituent quark model

Since the origin of the quark model hadrons have been considered to be built by constituent (massive) quarks. Nowadays it is widely recognized that because of the spontaneous breaking of chiral symmetry in the light quark sector at some momentum scale a constituent quark mass $`M(q^2)`$ appears. Once a constituent quark mass is generated such particles have to interact through $`SU(3)`$ Goldstone modes \[pion, kaon, eta and sigma (which simulates the two-pion exchange)\]. Explicit expressions of these potentials can be found elsewhere vij1 . In the heavy quark sector, chiral symmetry is explicitly broken and therefore these interactions will not appear.

For higher momentum transfer quarks still interact through gluon exchanges. Following de Rújula et al. ruju  the one-gluon-exchange (OGE) interaction is taken as a standard color Fermi-Breit potential. In order to obtain a unified description of light, strange and heavy mesons a scale dependent strong coupling constant has to be used bad . We parametrize this scale dependence by

$$\alpha _s(\mu )=\frac{\alpha _0}{ln\left(\frac{\mu ^2+\mu _0^2}{\mathrm{\Lambda }_0^2}\right)},$$
(1)
where $`\mu `$ is the reduced mass of the $`q\overline{q}`$ system and $`\alpha _0`$, $`\mu _0`$ and $`\mathrm{\Lambda }_0`$ are fitted parameters vij1 . This equation gives rise to $`\alpha _s0.54`$ for the light quark sector, a value consistent with the one used in the study of the nonstrange hadron phenomenology ente , and it also has an appropriate high $`Q^2`$ behavior pre0 , given a value of $`\alpha _s0.127`$ at the $`Z_0`$ mass pre1 . The $`\delta `$ function appearing in the OGE has to be regularized in order to avoid an unbound spectrum from below. To solve numerically the Schrödinger equation with this potential we use a flavor-dependent regularization vij1 .

The other nonperturbative property of QCD is confinement. Lattice QCD studies show that $`q\overline{q}`$ systems are well reproduced at short distances by a linear potential that is screened at large distances due to pair creation bali . One important question which has not been properly answered is the covariance property of confinement. While the spin-orbit splittings in heavy quark systems suggest a scalar confining potential luca , a significant mixture of vector confinement has been used to explain the decay widths of $`P`$-wave $`D`$ mesons suga . Such property, being irrelevant for the central part of the interaction, determines the sign and strength of the spin-orbit Thomas precession term which is important for the scalar mesons. Therefore, we write the confining interaction as an arbitrary combination of scalar and vector terms $`V_{CON}^{SO}(\stackrel{}{r}_{ij})=(1a_s)V_V^{SO}(\stackrel{}{r}_{ij})+a_sV_S^{SO}(\stackrel{}{r}_{ij})`$ where $`V_V^{SO}(\stackrel{}{r})[V_S^{SO}(\stackrel{}{r})]`$ is the vector (scalar) spin-orbit contribution.

## 3 Results

With the quark-quark interaction described above we have solved the Schrödinger equation for the different $`q\overline{q}`$ systems. Most of the parameters of the Goldstone boson fields are taken from the $`NN`$ sector. The eta and kaon cutoff masses are related with the sigma and pion one as explained in Ref. muli : $`\mathrm{\Lambda }[u(d)s]\mathrm{\Lambda }(ud)+m_s`$, $`m_s`$ being the strange quark current mass. The confinement parameters $`a_c`$ and $`\mu _c`$ are fitted to reproduce the energy difference between the $`\rho `$ meson and its first radial excitation and the $`J/\psi `$ and the $`\psi (2S)`$. The parameters involved in the OGE are obtained from a global fit to the hyperfine splittings well established in the Particle Data Group (PDG) pdgb . Finally, the relative strength of the scalar and vector confinement is fitted to the energy of $`a_1(1260)`$ and $`a_2(1320)`$, ordering that cannot be reproduced with a pure scalar confinement. We obtain $`a_s=0.777`$.

The spectra for the light pseudoscalar and vector mesons and for heavy mesons have been reported in Ref. dub . The agreement with experimental data is remarkable. Let us emphasize that with only 11 parameters we are able to describe more than 110 states.

Recently Belle and BaBar bel1 ; bel2 ; bel3  collaborations have reported new experimental measurements for the mass of the $`\eta _c(2S)`$. The average value from the different measurement is significantly larger than most predictions of constituent quark models and the previous experimental value of the PDG: $`M[\eta _c(2S)]=3594\pm 5`$ MeV. Such value cannot be easily explained in the framework of constituent quark models because the resulting $`2S`$ hyperfine splitting (HFS) would be smaller than the predicted for the $`1S`$ states eber ; hfs3 . Based on this fact some authors have claimed for an $`\alpha _s`$ coupling constant depending on the radial excitation.

We predict a value $`M[\eta _c(2S)]=`$ 3627 MeV, within the error bar of the last two Belle measurements, the ones obtained with higher statistics. Moreover the ratio 2S to 1S HFS is found to be 0.537, in good agreement with the experimental data, 0.479. The reason for this agreement can be found in the shape of the confining potential that also influences the HFS, the linear confinement being not enough flexible to accommodate both excitations pedr .

Other puzzling state is a narrow resonance around 2317 MeV, known as $`D_{S_J}^{}(2317)`$ reported by BaBar baba . This state has been confirmed by CLEO cleo  together with another possible resonance around 2460 MeV. Both experiments interpret these resonances as $`J^P=0^+`$ and $`1^+`$ states. This discovery has triggered a series of articles presenting alternative hypothesis arta . The most striking aspect of these two resonances is that their masses are much lower than expected. We obtain a mass of 2470.9 MeV for the $`J^P=0^+`$ and 2565.5 MeV for the $`J^P=1^+`$. They are far from the experimental data although the rest of the states ($`1^+`$,$`2^+`$,$`1^{}`$ and $`0^{}`$) agree reasonably well with the values of the PDG for both for the $`D`$’s and $`D_s`$’s states.

Concerning the scalar sector our results are shown in Table 1 (column three). We observe that for isovector states, there appears a candidate for the $`a_0(980)`$, the $`{}_{}{}^{3}P_{0}^{}`$ member of the lowest $`{}_{}{}^{3}P_{J}^{}`$ isovector multiplet. The other candidate, the $`a_0(1450)`$, is predicted to be the scalar member of a $`{}_{}{}^{3}P_{J}^{}`$ excited isovector multiplet. This reinforces the predictions of the naive quark model, where the $`LS`$ force makes lighter the $`J=0`$ states with respect to the $`J=2`$. The assignment of the $`a_0(1450)`$ as the scalar member of the lowest $`{}_{}{}^{3}P_{J}^{}`$ multiplet would contradict this idea, because the $`a_2(1312)`$ is well established as a $`q\overline{q}`$ pair. The same behavior is evident in the $`c\overline{c}`$ and the $`b\overline{b}`$ spectra, making impossible to describe the $`a_0(1450)`$ as a member of the lowest $`{}_{}{}^{3}P_{J}^{}`$ isovector multiplet without spoiling the description of heavy-quark multiplets. However, in spite of the correct description of the mass of the $`a_0(980)`$, the model predicts a pure light-quark content, what seems to contradict some experimental evidences. The $`a_0(1450)`$ is predicted to be also a pure light quark structure obtaining a mass somewhat higher than the experiment.

In the case of the isoscalar states, one finds a candidate for the $`f_0(600)`$ with a mass of 402 MeV, in the lower limit of the experimental error bar and with a strangeness content around $`8\%`$. The $`f_0(980)`$ and $`f_0(1500)`$ cannot be found for any combination of the parameters of the model. It seems that a different structure rather than a naive $`q\overline{q}`$ pair is needed to describe these states. The $`f_0(1500)`$ is a clear candidate for the lightest glueball amsl  and our results support this assumption. For the $`f_0(1370)`$ state (which may actually correspond to two different states blac ) we obtain two almost degenerate states around this energy, the lower one with a predominantly nonstrange content, and the other with a high $`s\overline{s}`$ content. Finally a state corresponding to the $`f_0(1710)`$ is obtained.

In the $`I=1/2`$ sector, as a consequence of the larger mass of the strange quark as compared to the light ones, our model always predicts a mass for the lowest $`0^{++}`$ state 200 MeV greater than the $`a_0(980)`$ mass. Therefore, being the $`a_0(980)`$ the member of the lowest isovector scalar multiplet, the $`\kappa (900)`$ cannot be explained as a $`q\overline{q}`$ pair. We find a candidate for the $`K_0^{}(1430)`$ although with a smaller mass.

## 4 The scalar mesons as tetraquark

Unlike the $`q\overline{q}`$ pairs tetraquark structures, suggested twenty years ago by Jaffe jaff , can couple to $`0^{++}`$ without orbital excitation and therefore could be serious candidates to explain the structure of the scalar mesons.

In this section we study tetraquark bound states by solving the Schrödinger equation using a variational method where the spatial trial wave function is a linear combination of gaussians. The technical details are given in Ref. vij2 . Due to the presence of the kaon-exchange there is a mixture among different configurations with the same isospin. In particular, in the isoscalar sector the configurations: $`[(qq)(\overline{q}\overline{q})]`$, $`[(qs)(\overline{q}\overline{s})]`$, and $`[(ss)(\overline{s}\overline{s})]`$ are mixed. The same happens in the isovector case for the configurations: $`[(qq)(\overline{q}\overline{q})]`$, and $`[(qs)(\overline{q}\overline{s})]`$, and in the $`I=1/2`$ case for the configurations: $`[(qq)(\overline{q}\overline{s})]`$, and $`[(qs)(\overline{s}\overline{s})]`$. In all cases $`q`$ stands for a $`u`$ or $`d`$ quark.

The results are shown in Table 1 (column four) where we present the lowest states for the three isospin sectors. As one can see, there appear two states, in the isoscalar and isovector sectors, with almost the same mass, although too high to be identified with the $`f_0(980)`$ and $`a_0(980)`$. In the $`I=1/2`$ sector, there appears a candidate to be identified with the $`\kappa (900)`$. It has been recently argued the possible importance of three-body forces arising from the confining interaction for those systems containing at least three quarks dimi . We have performed a calculation including a three-body confining term as the one reported in Ref. dimi  fixing its strength to reproduce the mass of the $`f_0(980)`$. The results are shown in Table 1 (column five) and as can be seen the degeneracy between the isoscalar and isovector states remains although their masses are now compatible with the experimental data. The lowest state of the isoscalar and $`I=1/2`$ sectors are almost not affected.

Using the same interaction and formalism we have calculated the $`D_{S_J}^{}(2317)`$ as a $`[(uc)(\overline{u}\overline{s})]`$ tetraquark. The result we obtain, M=2389 MeV, suggests that this state could also have a significant tetraquark component.

As a summary, we have found tetraquark bound states in the region of the light scalar mesons and in the $`D_{S_J}^{}(2317)`$. Our results suggest that some states, as it is the case of the $`a_0(980)`$ and $`f_0(600)`$, could present a significant mixture of $`q\overline{q}`$ and tetraquark structures, but it assigns a clear tetraquark structure to the $`f_0(980)`$ and the $`\kappa (900)`$. However, more accurate calculations including the exchange terms in the variational wave function which are negligible for the heavy-light tetraquarks and the explicit coupling to $`q\overline{q}`$ channels should be done before drawing any definitive conclusion.

This work has been partially funded by Ministerio de Ciencia y Tecnología under Contract No. BFM2001-3563 and by Junta de Castilla y León under Contract No. SA-109/01.


text:
# FUSE Observation of the Ultramassive White Dwarf PG 1658+441

## 1 Introduction

The DA white dwarf PG 1658+441 was discovered in the Palomar-Green survey (Green et al. 1986). Liebert et al. (1983) announced the presence of a magnetic field in this white dwarf and presented a detailed analysis of the energy distribution consistent with a 30,000 K pure hydrogen star. They measured a mean surface field of $`B_s=2.3`$ MG from the Zeeman splitting of the Balmer lines. The magnetic field of PG 1658+441 was studied in detail by Achilleos & Wickramasinghe (1989) using offset magnetic dipole models. They measured a field of $`B_d=3.0\pm 0.5`$ MG with a viewing angle between 40° and 50°. Schmidt et al. (1992) reported the first surface gravity measurement by fitting the Balmer line profiles with synthetic profiles including Zeeman splitting. They reached the surprising conclusion that PG 1658+441 has an abnormally high surface gravity of $`\mathrm{log}g=9.36\pm 0.07`$ corresponding to a stellar mass of $`1.31\pm 0.02`$ $`M_{}`$. The mass of PG 1658+441 appears much higher than the average mass of field white dwarfs ($`0.550.60M_{}`$). While searching for helium in ultramassive white dwarfs, Dupuis et al. (2002) analyzed the EUVE spectrum of PG 1658+441 and found evidence for a small quantity of atmospheric contaminants which was assumed to be helium.

The origin of PG 1658+441 and, more generally, of the population of hot ultramassive white dwarfs discovered in recent EUV/soft X-ray surveys is yet unknown. These objects appear to form a secondary peak in the white dwarf mass distribution near 1.25 $`M_{}`$ and a substantial fraction of these objects are also magnetic (Vennes 1999). As recently discussed by Dupuis et al. (2002), the problem posed with the mass of these stars is that they exceed the maximum mass of a C/O core ($`1.1M_{}`$) formed by the evolution of intermediate mass stars (Weidemann 2000). While it is not unconceivable to form white dwarfs with O/Ne core with mass exceeding 1.1 M through single star evolution channels (Nomoto 1984; Garcia-Berro et al. 1997), Schmidt et al. (1992) raised the issue that if PG1658+441 originates from a star with a mass in excess of 6–8 M, it should be relatively young and would be likely associated with a stellar cluster. Because PG1658+441 is not an obvious member of any cluster, Schmidt et al. (1992) proposed instead that PG1658+441 could be the product of the coalescence of a double-degenerate system. Such a system is expected to rotate rapidly and show photometric or magnetic field variability, but there is no evidence yet that PG1658+441 is variable (Liebert et al. 1983; Shtol’ et al. 1997). There is also the possibility that PG1658+441 originates from an older Ap star, and that a peculiar initial-mass to final-mass relation applies to the evolution of magnetic stars.

Why observe ultramassive white dwarfs with FUSE? Atmospheric parameters are usually obtained from a model atmosphere analysis of Balmer line profiles. With recent advance in FUV astronomy, spectroscopy of the hydrogen Lyman series has been made possible, and an independent diagnostic of atmospheric properties of hot white dwarfs is now available (see for example Finley et al. 1997, Barstow et al. 2001, and Barstow et al. 2003). The Lyman line series also allows for a test of line broadening theory in the high density conditions encountered in the atmosphere of ultramassive white dwarfs. The high spectral resolution ($`R=1700020000`$) achieved with FUSE allows for the detection of heavy element lines and helps better characterize the atmospheric composition of ultramassive white dwarfs. The detection of helium in the EUVE spectrum of GD 50 by Vennes et al. (1996) indicates that trace elements may be present in ultramassive white dwarf atmospheres. Although a study by Dupuis et al. (2002) shows that GD 50 is quite unique, the case is still open due to the limited sensitivity of EUVE. FUSE observations of this class of objects are needed to provide an independent determination of their stellar parameters.

We present an analysis of the FUSE spectrum of the ultramassive white dwarf PG 1658+441. In section 2, we describe the observations and data reduction. In section 3, we present the model atmosphere analysis and describe the abundance measurements of C, N, Si, and P. We also analyze the Lyman line profiles which clearly show the quasi-molecular satellite features of H$`\stackrel{+}{_2}`$ in the wings of Lyman-$`\beta `$ and Lyman-$`\gamma `$. This result is interesting in view of the high surface gravity and relatively high effective temperature of the star. We conclude by discussing the implications of our results on the origin of ultramassive white dwarfs and we summarize our findings.

## 2 Observations and Data Reduction

PG 1658+441 was observed as part of our FUSE Cycle 2 GI program B122. The spacecraft and the instruments are described in Moos et al. (2000, 2002) and an initial account of the calibration is given in Sahnow et al. (2000). FUSE provides spectra from 905 to 1187 Å with a velocity resolution of 15–20 km s-1 in four independent channels (SiC1, SiC2, LiF1, LiF2) each split in two segments (A and B). The sequence of 16 exposures was initiated on 17 July 2001 at 13:42:45 (UT) for a total exposure time of 28483 seconds. The target was centered in the LWRS aperture and the observations were made in TTAG mode. We have reprocessed the data with CALFUSE version 2.2.2. We found the target to be out of the aperture during 4 out of the 16 exposures in the LiF2 channel.

We have coadded the spectra of each channel from all exposures making sure we selected only photometric exposures. The coaddition is made using the program “cf\_combine” provided with CALFUSE. Figures 1 and 2 show the combined spectra binned by 24 pixels for each segment of the FUSE spectrometer. The Lyman lines series from the white dwarf photosphere clearly dominate the spectrum. A broad depression centered on 1150 Å was initially evident in the LiF1B spectrum, but it was not observed in the LiF2A spectrum which covers a similar wavelength range. This instrumental artifact, nicknamed the “worm”, is thought to be due to blockage of the spectral image by a grid wire above the surface of the micro-channel plate detector. We have derived a “worm” correction using the LiF2A segment as a reference spectrum. We smoothed the ratio of the LiF2A to LiF1B spectra using cubic splines and corrected the LiF1B spectrum by multiplying it by this smoothed ratio.

## 3 Analysis of the FUV spectrum of PG 1658+441

In order to analyze the spectrum of PG 1658+441, we used the model atmosphere code TLUSTY(v. 195) (Hubeny 1988; Hubeny & Lanz 1995) and the spectral synthesis code SYNSPEC version 45 (SYNSPEC45) which incorporates hydrogen Lyman line opacities (Lemke 1997), and SYNSPEC version 48 (SYNSPEC48) which incorporates new quasi-molecular satellite line opacities of H2 and H$`\stackrel{+}{_2}`$ computed by N. Allard for Lyman-$`\alpha `$, Lyman-$`\beta `$, and Lyman-$`\gamma `$. The quasi-molecular opacity tables of Lyman-$`\beta `$ were specifically recomputed for a temperature of 30,000 K, more appropriate for PG 1658+441. The general method is described by Allard et al. (1998) who computed theoretical line profiles for Lyman-$`\beta `$ including quasi-molecular satellites due to H$`\stackrel{+}{_2}`$ with a variable dipole.

### 3.1 Quasi-Molecular Satellites

First, we computed a pure hydrogen model atmosphere in non-LTE adopting published atmospheric parameters of PG 1658+441, $`T_{\mathrm{eff}}=30510\pm 200`$K and $`\mathrm{log}g=9.36\pm 0.07`$ (Schmidt et al. 1992). Next, we calculated a synthetic spectrum using SYNSPEC45 and covering the FUSE spectral range. This generic model is normalized to match the absolute flux level in LiF1A and LiF1B. The normalization corresponds to a $`V`$ magnitude of about 14.79 and agrees with the same normalization based on International Ultraviolet Explorer spectra (Dupuis et al. 2002).

Figures 1 and 2 show a comparison between the model spectrum (solid lines) and the FUSE spectra. The generic model seems in good agreement with the observed line profiles, even though we have not yet attempted a formal fit of the Lyman line profiles at this point (see section 3.2). However, a closer inspection of Figures 1 and 2 shows discrepancies in the Lyman-$`\beta `$ and Lyman-$`\gamma `$ red wings.

Broad absorption troughs are apparent at 1058 Å and 1076 Å in the red wing of Lyman-$`\beta `$ of the LiF1A, LiF2B, SiC1A, and SiC2B spectra. The troughs are attributed to H$`\stackrel{+}{_2}`$ quasi-molecular satellite features previously detected in cooler white dwarfs by the Hopkins Ultraviolet Telescope (HUT) (Koester et al. 1996), ORFEUS (Koester et al. 1998), and FUSE (Wolff et al. 2001; Hébrard et al. 2002). In addition, the features possibly present on the red wing of Lyman-$`\gamma `$ do correspond to quasi-molecular features detected in the HUT spectra of WD1031$``$114 and WD0644+375 by Koester et al. (1998) and in the FUSE spectrum of CD$``$38o10980 by Wolff et al. (2001) and included for the first time in model comparisons presented by Hébrard et al. (2003). PG 1658+441 is the hottest white dwarf so far to show these FUV spectral features.

The validity of the Balmer and Lyman line profiles is questionable because in high-gravity stars the Lyman line wings form at depths where the electronic density is between $`10^{17}10^{18}`$ cm-3, near the limit of Lemke’s tables (Lemke 1997). Moreover, the tables should not be applied directly to PG 1658+441 because they do not include the satellite lines from quasi-molecules of H$`\stackrel{+}{_2}`$, but the good agreement between the synthetic spectra and the observations (Figs. 1 and 2) away from the quasi-molecular satellite features suggests that they are a good approximation.

The emergence of quasi-molecular satellite lines of H$`\stackrel{+}{_2}`$ in PG 1658+441 is possibly due to its high surface gravity which implies larger particle densities in the line forming region. Because PG 1658+441 is hot, hydrogen should be more ionized but the high surface gravity implies larger pressures (or larger electronic densities) in the atmosphere which in turn favor recombination. Because of this effect, we find the number density of shydrogen atoms to be comparable with the one found in the atmospheres of DA in which the quasi-molecular satellites lines have been previously detected. These stars are typically cooler with $`T_{\mathrm{eff}}`$ up to about 20,000K and less massive with $`\mathrm{log}g`$ close to 8. Because the quasi-molecular opacity is proportional to the proton and H I densities, we predict that the H$`\stackrel{+}{_2}`$ quasi-molecular satellite lines become apparent at higher effective temperatures in WD with larger surface gravities.

Finally, adopting the same effective temperature and surface gravity ($`T_{\mathrm{eff}}=30,510`$ K, $`\mathrm{log}g=9.36`$), we computed a synthetic spectrum using SYNSPEC48 which includes quasi-molecular satellite line opacities. Figures 1 and 2 show a comparison between the data and the synthetic spectra which include the quasi-molecular satellite lines of Lyman-$`\beta `$ and Lyman-$`\gamma `$ (dashed lines). The predicted satellite lines do not agree well with the data as the H$`\stackrel{+}{_2}`$ satellites at 1058 Å and 1076 Å are somewhat too deep and shifted toward shorter wavelengths. On the other hand, the predicted Lyman-$`\gamma `$ satellites are too weak. There is also a subtle kink predicted in the Lyman-$`\beta `$ profile near 1037 Å which is possibly present in the data and better seen in Figure 3 (see next section). This preliminary comparison is certainly promising but we must investigate whether or not the agreement could be improved by doing a formal fit of the FUSE spectrum.

### 3.2 Atmospheric Parameters

To determine the atmospheric parameters more accurately, we have computed a fine grid of model atmospheres and synthetic spectra. The spectra were computed with SYNSPEC48 which includes the quasi-molecular opacities of H$`\stackrel{+}{_2}`$ for Lyman-$`\beta `$ and Lyman-$`\gamma `$. The grid, appropriate for determining the atmospheric parameters of PG 1658+441, includes effective temperatures from 25,000 to 35,000K in steps of 1000K and surface gravities from 7.5 to 9.5 ($`\mathrm{log}g`$) in steps of 0.1 dex.

Each of the spectroscopic channels showing the Lyman lines are fitted separately so that we obtain independent measurements of the stellar parameters. We used a $`\chi ^2`$–minimization technique to find the best-fit parameters, $`T_{\mathrm{eff}}`$ and $`\mathrm{log}g`$. We excluded regions contaminated by airglow, interstellar absorption lines, as well as detector edges where the spectra are of poor quality. We neglected Zeeman splitting which could produce shifts of the order of 1 Å or less in the core of Lyman-$`\beta `$ and higher Lyman lines. The approximation does not affect the results because the line cores are already excluded from the fit. We also binned the spectra by 20 pixels (which corresponds to about 0.13Å) to improve the signal-to-noise ratio. For the fit of the Lyman-$`\beta `$ line, we restrict the fit to wavelengths between 1005Å and 1050Å to exclude regions affected by the quasi-molecular satellites of Lyman-$`\beta `$ and Lyman-$`\gamma `$. Including these regions had the consequence of driving the fits toward spurious solutions. The spectrum normalization (or scale factor) is also left as a free parameter in the fit. An overall wavelength shift of $`0.4`$ Å was applied to the model spectra to align optimally the models with the data. Table 1 summarizes the results. The quoted error bars are computed from the projection of the 3$``$$`\sigma `$ confidence contours in the $`T_{\mathrm{eff}}`$$`\mathrm{log}g`$ plane. A representative example of a fit is shown on Figure 3 for the Lyman-$`\beta `$ line in the LiF1A segment which shows how well the line is fitted up to where the quasi-molecular satellites become important.

As expected, there is some scatter in the parameter values derived from the different channels and segments which gives us an assessment of the systematic errors. The measured parameter values are consistent within the quoted error bars although the surface gravity measured from SiC2B appears to be greater. This last measurement is less reliable, as reflected by the larger error, because the Lyman-$`\beta `$ profile from SiC2B is missing most of the blue wing (the wavelength coverage ends at 1017 Å). Minor differences observed between different channels are probably due to subtle flux calibration effects such as light blocking by grid wires or sensitivity degradation which are not perfectly accounted for in the data reduction. Nonetheless, the values obtained by taking the weighted mean of the six measurements are in agreement with those of Schmidt et al. (1992) although we seem to favor slightly lower values of the gravity and effective temperature.

Our results confirm the analysis of Schmidt et al. (1992) which was possibly affected by Zeeman splitting observed in the Balmer line cores. In general, they validate gravity measurements in low-field white dwarfs which are based on the modeling of spectral line wings computed at zero magnetic field. It is also interesting to note that Barstow et al. (2001) obtained a slightly lower value of the surface gravity based on the FUV HUT spectrum of the ultramassive white dwarf GD 50 compared to values based on optical spectra, but such differences do not appear worst than for less massive DA stars in their sample. Note that Friedrich et al. (1996) suggested that Stark broadening due to an electric field parallel to the magnetic field may be important in magnetic white dwarfs and that it may explain the broadness of PG 1658+441 Balmer lines. However, their best fit to an optical spectrum of PG 1658+441 is quite poor and may indicate that the effect is overestimated. We have not considered this effect.

### 3.3 Metal Abundance Upper Limits

A search for heavy elements in PG 1658+441 is motivated by the detection of helium in the atmosphere of GD 50 (Vennes et al. 1996) and by the possible detection of helium or other contaminants in the extreme-ultraviolet spectrum of PG 1658+441 (Dupuis et al. 2002). These observations suggest that other ultramassive white dwarfs may have some metals in their atmospheres. These metals may provide additional clues required to better understand the origin of these stars as they may be relics of the original composition of the ultramassive white dwarfs. Various scenarios for producing ultramassive white dwarfs predict hydrogen poor atmospheres. Therefore, it is important to understand the nature of ultramassive white dwarfs which appear to have hydrogen-rich atmospheres.

Thus, we have searched for the presence of metals in the atmosphere of PG 1658+441 but have found none of the usual suspects. Nevertheless, we have used the spectra to set upper limits on elements likely to be detected in FUSE spectra of white dwarfs, that is C, N, Si, and P. We computed non-LTE synthetic spectra for each element selected over a range of abundances varying from X/H = $`10^9`$ to $`10^5`$ in steps of 0.25 dex. The strongest predicted lines (as shown on Fig. 4) are used to compute a 3–$`\sigma `$ upper limit on the abundance. Note that the upper limits are measured at the rest wavelengths of the selected transitions, and that we corrected the wavelength scale for a predicted gravitational redshift of 211 km s-1. There is an additional uncertainty related to our lack of knowledge of the radial velocity of PG 1658+441. Figure 4 shows the spectra on a velocity scale with $`v=0`$ corresponding to the laboratory wavelength. The noise appears uniform over a range of $`\pm 500`$ km s-1. The linear Zeeman effect is not taken into account in the synthetic spectra, but it could shift the $`\sigma `$ components by $`\pm 360`$ km s-1 (assuming a surface average field of 2.3 MG). Limiting the analysis to the un-shifted component, this effect adds $`0.5`$ dex to the abundance upper limits. However, the spectra do not show any evidence for shifted components. The results are shown in Figure 4 and summarized in Table 2. The upper limits are rather stringent and shows that PG 1658+441’s atmosphere is metal poor and is consistent with having a pure hydrogen atmosphere.

## 4 Discussion and Summary

The fact that no trace elements are detected in the hydrogen atmosphere of PG 1658+441 is perhaps not so surprising. At such a high surface gravity we do not expect radiative acceleration to support detectable abundance of photospheric heavy elements (Chayer et al. 1995). On the other hand, the hydrogen-rich composition of PG 1658+441 is quite interesting. The surface of white dwarfs emerging from the high-end of the main-sequence progenitor mass range (8–10 M) is expected to be depleted of hydrogen and even probably of helium (Garcia-Berro et al. 1997). It may suggest that there was a sufficient amount of hydrogen left deeper in the envelope to form, over time, a hydrogen-dominated atmosphere. On the other hand, the progenitor of PG 1658+441 may have been an early-A or late-B magnetic star which developed a higher degenerate core mass because of the presence of a magnetic field.

Alternatively, the hydrogen-rich atmosphere of PG 1658+441 may suggest that hydrogen was accreted from the interstellar medium as suggested by Segretain et al. (1997). However, in the present case, accretion of interstellar hydrogen may be inhibited by the presence of a magnetic field. Using the formalism of Angel et al. (1981) and assuming that PG 1658+441 is accreting at the Eddington rate as argued by Koester (1976), we calculated that a mean equatorial field of about 4136 G would be sufficient to prevent the accretion flow from reaching the surface of the star (assuming a particle density of 1 cm-3 and a velocity of 50 km s-1). With a mean surface field of 2.3 MG (Schmidt et al. 1992), it is unlikely that PG 1658+441 would accrete significantly. Therefore, it is unlikely that the star entered the white dwarf cooling sequence with an atmosphere devoid of hydrogen.

The detection of quasi-molecular satellite absorption in the Lyman-$`\beta `$ and Lyman-$`\gamma `$ wings of PG 1658+441 shows that these features can be observed even in hot white dwarfs and that they should be included in the modeling of Lyman lines up to effective temperatures higher than assumed in the past. Our analysis indicates that the treatment of quasi-molecular satellite opacities could be improved, particularly in the range of ionic densities encountered in the atmosphere of ultramassive white dwarfs. Our current approach does indeed have some limitations. For instance, the quasi-molecular satellite profiles we use are based on an expansion in density (Allard et al. 1994) valid for low densities of hydrogen atoms ($`n10^{18}`$ cm-3). This approximation is appropriate for white dwarfs with cooler and lower surface gravities but is less acceptable for ultramassive white dwarfs. This might explain in part why the fits shown in previous investigations (Koester et al. 1998; Wolff et al. 2001; Hébrard et al. 2003, for examples) appear somewhat better. Another improvement would be to include a variable dipole moment in the calculations of the satellites of Lyman-$`\gamma `$ as suggested by Allard et al. (2003). This may enhance the strength of the Lyman-$`\gamma `$ satellites and help improve the fit shown on Fig. 2. At this stage, we have not included the quasi-molecular satellites of the Lyman lines in our model atmosphere calculations because those features are rather weak in comparison with the Stark wings at the relatively high effective temperature of PG 1658+441. In a future investigation, we plan to include the quasi-molecular satellites in the model atmosphere and interpolate in tables of profiles computed for densities appropriate for ultramassive white dwarfs.

In summary, we have presented an analysis of the FUSE spectrum of the ultramassive white dwarf PG 1658+441 and we find that:

* The FUSE spectra are well represented by a pure hydrogen non-LTE model spectrum computed with the atmospheric parameters determined from optical spectroscopy.
* Lyman-$`\beta `$ satellite absorption has been detected at 1058 Å and 1076 Å making PG 1658+441 the hottest star in which these features are detected.
* Upper limits on the abundance of C, N, Si, and P imply that PG 1658+441 has a highly pure hydrogen atmosphere.

This work is based on observations made with the NASA-CNES-CSA Far Ultraviolet Spectroscopic Explorer. FUSE is operated for NASA by the Johns Hopkins University under NASA contract NAS5-32985. This research is funded by NASA grants NAG5-11570 and NAG5-11844.


text:
# A New Monomeric Interpretation of Intrinsic Optical Bistability Observed in 𝑌𝑏^limit-from3+ - Doped Bromide Materials

## Abstract

We present a mechanism able to show intrinsic bistable behaviour involving single $`Yb^{3+}`$ ions embedded into bromide lattices, in which intrinsic optical bistability (IOB) has been observed. The mechanism is based on the experimentally found coupling between the $`Yb^{3+}`$ ion and the totally symmetric local mode of vibration of the $`[YbBr_6]^3`$ coordination unit. The model reproduces the IOB observed in $`CsCdBr_3:1\%Yb^{3+}`$ and allows to understand the experimentally found presence of the phenomenon in the other bromides, but its absence in $`Cs_3Lu_2Cl_9:Yb^{3+}.`$

A physical system without mirrors of an external cavity, exhibiting two stable outputs for a single input intensity of an exciting laser, is said to exhibit intrinsic optical bistability (IOB) . Hehlen et al.  have discovered that $`Yb^{3+}`$ \- doped $`Cs_3X_2Br_9`$ \[$`X=Y,Lu`$\] and $`CsCdBr_3`$ show IOB when exciting into the $`{}_{}{}^{2}F_{7/2}^{}(0)^2F_{5/2}(0^{})`$ and $`{}_{}{}^{2}F_{7/2}^{}(0)^2F_{5/2}(2^{})`$ electronic transitions of the $`Yb^{3+}`$ ion with a narrow band laser of suitable near-infrared (NIR) frequency. On the contrary, no IOB was observed under excitation into the $`{}_{}{}^{2}F_{7/2}^{}(0)^2F_{5/2}(1^{})`$ transition . When doped in these crystalline bromide materials some $`Yb^{3+}`$ ions form strongly interacting $`Yb^{3+}Yb^{3+}`$ dimers . Due to this interaction, excitation of the $`Yb^{3+}`$ $`{}_{}{}^{2}F_{7/2}^{}^2F_{5/2}`$ transition (around $`10000cm^1`$) results in both, $`{}_{}{}^{2}F_{5/2}^{}^2F_{7/2}`$ emission of single $`Yb^{3+}`$ ions in the NIR spectral range, as well as $`Yb^{3+}Yb^{3+}`$ dimer emission in the visible (VIS) at twice the resonance frequency of the single $`Yb^{3+}`$ ion  through a cooperative luminescence mechanism . The IOB phenomenon consists of a temperature and excitation power dependent hysteretic behaviour of the transmitted, NIR and upconverted VIS emission intensity (with the system switching output intensity at two critical incident excitation powers). To date, it is discussed if the phenomenon is due to the presence of strongly interacting $`Yb^{3+}Yb^{3+}`$ dimers or if it is a single $`Yb^{3+}`$ ion property. The theories relying on the former interpretation are usually referred to as “dimeric”while those supporting the latter idea as “monomeric”. Dimeric theories substantially rely on the existence of intra-pair electronic coupling of $`Yb^{3+}Yb^{3+}`$ dimers . According to these theories, IOB should be observed *a fortiori* in $`Yb^{3+}`$ \- doped $`Cs_3Lu_2Cl_9`$, as this material is isostructural to $`Yb^{3+}`$ \- doped $`Cs_3X_2Br_9`$ \[$`X=Y,Lu`$\]  and exhibits a stronger electronic exchange interaction for $`Yb^{3+}Yb^{3+}`$ dimers . However, no experimental evidence for IOB was found in $`Yb^{3+}`$-doped $`Cs_3Lu_2Cl_9`$ . Furthermore, in $`Cs_3X_2Br_9:Yb^{3+}`$ ($`X=Y,Lu`$) only a few $`Yb^{3+}`$ ions form dimers at low concentration of impurity ions. Nevertheless, both bromide lattices exhibit IOB even at a low $`Yb^{3+}`$ \- doping level . A monomeric mechanism to account for IOB in these systems has been proposed by Gamelin et al.  who presented a phenomenological model, according to which IOB depends on a single - $`Yb^{3+}`$ ion property in these materials. In this model, IOB is attributed to the interplay of a laser heating effect and a strongly increasing, non-linear dependence of the material’s absorbance with the internal temperature. With its phenomenological nature, this theory succeeds in reproducing many experimental trends, but does not provide the microscopic mechanism of IOB. In this letter, we present the first microscopic, monomeric mechanism that allows a single active ion to exhibit all essential aspects observed in IOB. To date, the microscopic theories presented in literature describing how a collection of two-level atoms can exhibit IOB are based on electronic coupling between two or more near absorbers . In these cases, the electronic interaction between the coupled absorbers provides the feedback loop required to obtain intrinsic bistability. Together with the intrinsic two-level atom non linearity and the use of single-atom density matrix approach (an indispensable step in order to obtain non linear equations), this feedback loop gives rise to inversion dependent atomic resonance frequency renormalization, and thus possible IOB. With a few exceptions, the role played by the host lattice, and in particular of the first coordination sphere of the $`Yb^{3+}`$, in the mechanism responsible for IOB in $`Yb^{3+}`$-doped bromide host materials has generally been regarded as secondary. However, Hehlen et al. measured different internal temperatures for the system on the two branches of the hysteresis loop . This suggests that with the first switch occurring, the lattice undergoes heating. Additionally, Lüthi  has found that $`Cs_3Lu_2Br_9:10\%Yb^{3+}`$ shows bistable behaviour even when excited with a chopped laser and off-times as long as $`75msec`$. This means that the system prepared on the upper branch of the bistability region does not decay to the lower branch for times as long as $`75msec`$ after turning off the exciting laser beam. If IOB were due to a cooperative mechanism based on electronic $`Yb^{3+}Yb^{3+}`$ coupling, and the feedback loop provided by the same two ytterbium ions, it would be hard to justify this behavior. After $`75msec`$, all ions in their excited state are expected to have decayed to the ground state, and there seems to be no physical reason for the system to save memory of its past history. The fact that the nuclear relaxation times of the lattice are much longer than the electronic ones, suggests that the lattice nuclear degrees of freedom play a relevant role in the bistability mechanism. Supported by the findings of a Raman spectroscopy study on $`CsCdBr_3:1\%Yb^{3+}`$  showing strong coupling between the electronic $`{}_{}{}^{2}F_{7/2}^{}(0)^2F_{5/2}(2^{})`$ transition of $`Yb^{3+}`$ and the totally symmetric vibrational local mode of the $`[YbBr_6]^3`$ coordination unit (into which $`Yb^{3+}`$ is embedded in all bromide lattices showing IOB), we propose that this vibrational coupling provides the indispensable feedback loop to observe IOB (involving only a single $`Yb^{3+}`$ ion).

Our physical system consists of the electronic structure of a single $`Yb^{3+}`$ ion and the totally symmetric vibrational local mode $`A_{1g}`$ of the $`[YbBr_6]^3`$ octahedron. We treat the former as a two-level system ($`|a`$ and $`|b`$ ground and excited states, respectively, $`\mathrm{}\omega _0`$ energy difference, $`|ij|=\sigma _{ij}`$ $`i,j=a,b`$) and the latter as a quantum harmonic oscillator ($`a`$ and $`a^+`$ annihilation and creation operators, respectively, and angular frequency $`\omega _k`$). We write the free Hamiltonian $`H_0`$ of the system as

$$H_0=H_{Yb}+H_k=\mathrm{}\omega _0\sigma _{bb}+\mathrm{}\omega _ka^+a$$
(1)
For the interaction Hamiltonian $`V_{Ybk}`$ between the ion and the $`A_{1g}`$ local mode we choose a bilinear operatorial form (in order for the oscillation centre to depend on the electronic state in which the ion is )

$$V_{Ybk}=ϵ(\sigma _{bb}\sigma _{aa})(a+a^+)$$
(2)
When a coherent classical electromagnetic field is applied to the ion, using the rotating wave and the electric dipole approximations, the coupling between the $`Yb^{3+}`$ ion and the external field is described by the observable $`V_{YbL}=\frac{\mathrm{}\mathrm{\Omega }}{2}(\sigma _{ba}e^{i\omega t}+\sigma _{ab}e^{i\omega t})`$, where $`\mathrm{\Omega }`$ and $`\omega `$ stand for the Rabi and the laser frequency. Denoting the total Hamiltonian operator of the system with $`H=H_0+V_{Ybk}+V_{YbL}`$, and taking into account both the ion spontaneous emission and the dissipative interaction of the ion and the local mode with the bulk modes of the lattice (described by a suitable superoperator $`S`$), the dynamics of the system obeys the generalized Von-Neumann equation

$$\dot{F}=\frac{i}{\mathrm{}}[H,F]+SF$$
(3)
where $`F`$ is the density operator of the whole system. $`F(t)`$ can be written as $`F(t)=Tr_kF(t)Tr_{Yb}F(t)+F^{Ybk}(t)`$, where $`Tr_kF`$ and $`Tr_{Yb}F`$ are the reduced density operators of the ion and the local mode, respectively, while $`F^{Ybk}(t)`$ represents the correlations existing at the time $`t`$ between $`Yb`$ and $`k`$. Assuming $`V_{Ybk}H_0`$, we approximate $`F(t)`$ with its factorized part

$$F(t)Tr_kF(t)Tr_{Yb}F(t)$$
(4)
This is the main hypothesis of our model. In this way, analogously to other theories , we *assume* that single-ion density matrix approach is allowed. A necessary condition for this approach to be valid is that $`V_{Ybk}`$ is much smaller than the free Hamiltonian $`H_0`$ (later on we will show this is the case). It must be remarked that the feasibility of this treatment - indispensable to have non linear equations and thus possible IOB - poses a delicate and lively debated  theoretical problem questioning even the existence of a microscopic explanation of some IOB phenomena . We justify assumption (4) with the ability of our model to reproduce the observed IOB phenomena. Inserting eq.(4) in eq.(3) and doing the trace over $`k`$ and $`Yb`$, we obtain two coupled master equations in $`Tr_kF(t)`$ and $`Tr_{Yb}F(t)`$. Performing the variable change $`\rho _{ab}^{}=\rho _{ab}e^{i\omega t}`$, it can be proved that the ion dynamic equations in $`\dot{\rho }_{bb}`$ and $`\dot{\rho ^{}}_{ab}`$ differ from the usual optical Bloch equations in the rotating frame for the presence of an effective resonance frequency $`\omega _0^{eff}(t)=\omega _0+2\frac{ϵ}{\mathrm{}}a+a^+(t)`$, linearly dependent on the mean value of the local mode operator $`\left(a+a^+\right)`$. This proves that the dynamic response of the ion to the applied field depends on the state of the local mode (in its turn depending on the ion state). It is easy to show that for stationary conditions (when setting all temporal derivatives to zero) the resonance frequency can be written as

$$\omega _0^{eff}=\omega _0\beta (\rho _{bb}\rho _{aa})$$
(5)
where

$$\beta =4\frac{ϵ^2}{\mathrm{}^2}\frac{\omega _k}{\omega _k^2+\gamma ^2}$$
(6)
with $`\gamma `$ being the relaxation rate of the local mode. We thus find that the overall effect of the coupling between the electronic transition of the $`Yb^{3+}`$ ion and the $`A_{1g}`$ local mode of the $`[YbBr_6]^3`$ unit, consists of an ion resonance frequency renormalization linearly dependent on inversion ($`\rho _{bb}\rho _{aa}`$). Introducing the absorption parameter $`s=2\rho _{bb}`$ ($`0s1`$), it is possible to prove that, under stationary conditions, $`s`$ and the incident excitation power $`P`$ must obey the equation

$$c_3s^3+c_2s^2+\left(c_1+bP\right)sbP=0$$
(7)
with

$`c_3=4\left({\displaystyle \frac{\beta }{\mathrm{\Gamma }_h}}\right)^2c_2=8{\displaystyle \frac{\beta }{\mathrm{\Gamma }_h}}{\displaystyle \frac{\delta +\beta }{\mathrm{\Gamma }_h}}c_1=4\left({\displaystyle \frac{\delta +\beta }{\mathrm{\Gamma }_h}}\right)^2+1`$
$`b={\displaystyle \frac{3e^2}{\pi mcϵ_0}}{\displaystyle \frac{F_{ab}}{\mathrm{}\omega _0}}cos^2\phi {\displaystyle \frac{1}{\mathrm{\Gamma }_1\mathrm{\Gamma }_h}}`$ (8)
where $`\mathrm{\Gamma }_h`$ is the homogeneous linewidth of the electronic transition, $`\delta =\omega _0\omega `$ is the detuning of the laser frequency, $`F_{ab}=(2m\omega _0/3\mathrm{})|a|𝐫|b|^2`$ stands for the oscillator strength of the transition $`|a|b`$, $`\phi `$ is the angle between the laser polarization versor and $`a|𝐫|b`$, and $`\mathrm{\Gamma }_1`$ is the ion longitudinal decay time. Eq.(7) is a third degree equation in the unknown $`s`$ dependent on many parameters. It can be shown that the condition for the existence of a $`P`$-range inside which, for any value of $`P`$, two stable solutions for $`s`$ are found (IOB) is

$$\frac{\beta }{\mathrm{\Gamma }_h}\frac{3\sqrt{3}}{2}$$
(9)
This provides a threshold value of the local mode - $`Yb^{3+}`$ ion interaction strength required to observe intrinsic bistability compatibly with our microscopic feedback mechanism. Moreover, eq.(7) shows that, for an increasing value of the bistability parameter $`\beta /\mathrm{\Gamma }_h`$, the incident power bistability range moves toward higher values of $`P`$ and, in addition, becomes wider. We now use our theory to reproduce the two incident powers at which switching occurs measured in $`CsCdBr_3:1\%Yb^{3+}`$ . Guillot-Noël et al. successfully performed an analogous simulation to test their dimeric theory . In our study we use: $`\delta =0`$ (perfect resonance), $`\omega _0=10602.8cm^1`$ , $`\mathrm{\Gamma }_h=0.6cm^1`$ , $`1/\mathrm{\Gamma }_1=7.810^4s`$ , $`F_{ab}=10^6`$ and, additionally, we take $`cos^2\phi =1/2`$, assuming that the incident laser beam is unpolarized . Using the bistability parameter to optimise our simulation, good agreement is achieved between the theoretical prediction and the experimentally measured values  for the two incident powers at which switching occurs, as shown in Fig. 1.

Using the optimised value of $`\beta /\mathrm{\Gamma }_h=7.5`$, the ion-local mode coupling constant $`ϵ`$ \[eq.(2)\] can be estimated. From eq.(6) and taking into account the fact that the lattice nuclear relaxation times are typically very long compared to $`1/\omega _k`$, we get $`ϵ\frac{\mathrm{}}{2}\sqrt{\beta \omega _k}`$. With $`\omega _k=191cm^1`$ , we obtain $`ϵ=14.66cm^1`$ for $`CsCdBr_3:1\%Yb^{3+}`$. This value confirms our *a priori* assumption that $`V_{Ybk}`$ is much weaker than $`H_0`$. With this value for $`ϵ`$, the displacement of the oscillation centre for the local mode $`k`$, for an $`Yb^{3+}`$ ion undergoing the transition from $`|a`$ to $`|b`$, as implied by the operatorial form assumed for $`V_{Ybk}`$ \[eq.(2)\], turns out to be in the order of $`1\%`$ . This result agrees with typical values found for lattice point defects due to incorporation of rare earth impurity ions , and confirms that the coupling constant $`ϵ`$ we get from our model is of the correct order of magnitude.

Another relevant test for our theory is its ability to predict the presence of IOB in $`Cs_3X_2Br_9:Yb^{3+}`$ \[$`X=Y,Lu`$\], but the absence of the phenomenon in $`Cs_3Lu_2Cl_9:Yb^{3+}`$ . In the bromide compounds, the angular frequency $`\omega _k`$ is $`190cm^1`$ , similar to the one found in $`Cs_3Er_2Br_9`$ . From comparison of the absorption and emission spectra, we estimate $`\mathrm{\Gamma }_h`$ and $`ϵ`$ to be comparable to the ones measured in $`CsCdBr_3:Yb^{3+}`$ . These assumptions made, we expect the bistability parameter $`\beta /\mathrm{\Gamma }_h`$ for $`Cs_3X_2Br_9:Yb^{3+}`$ \[$`X=Y,Lu`$\] to have a comparable value to the one calculated for $`CsCdBr_3:1\%Yb^{3+}`$, and can understand why IOB is found in the former compounds. For the chloride, $`\omega _k`$ is $`285cm^1`$ . From inspection of the absorption and emission spectra, we estimate $`\mathrm{\Gamma }_h`$ to fall in the range $`0.30.7cm^1`$ , and the $`Yb^{3+}`$ \- $`A_{1g}`$ local mode coupling constant $`ϵ`$ to be about 5 times smaller than in bromide compounds. Based on these values, the bistability parameter $`\beta /\mathrm{\Gamma }_h`$ is estimated to be lower than the critical threshold of $`3\sqrt{3}/2`$ required for IOB. We finally notice that absorption spectra of $`Cs_3Lu_2Br_9:Yb^{3+}`$ show that the $`Yb^{3+}`$ ion $`{}_{}{}^{2}F_{7/2}^{}(0)^2F_{5/2}(1^{})`$ transition exhibits a very weak coupling to the $`A_{1g}`$ local mode compared to $`00^{}`$ and $`02^{}`$. This explains why no IOB is observed under excitation into this peculiar transition.

In summary, we have presented the first microscopic mechanism according to which a single $`Yb^{3+}`$ ion is able to exhibit IOB when doped in $`Cs_3X_2Br_9`$ \[$`X=Y,Lu`$\] or $`CsCdBr_3`$. The mechanism couples the electronic two-level system of the $`Yb^{3+}`$ ion with the totally symmetric local mode of the vibration of the $`[YbBr_6]^3`$ coordination unit. The theory has been used to reproduce basic aspects of the IOB observed in $`CsCdBr_3:1\%Yb^{3+}`$. We emphasize that our theory not only allows to understand, why IOB has been observed in $`CsCdBr_3:1\%Yb^{3+}`$ and $`Cs_3X_2Br_9`$ \[$`X=Y,Lu`$\], but also, why it is absent in the chloride compound $`Cs_3Lu_2Cl_9:Yb^{3+}`$. These results provide new insight into the microscopic mechanism underlying the IOB phenomenon and guidelines for expanding the material basis in the search for IOB systems, and their potential applications.

We gratefully acknowledge inspiring discussions with O. Guillot-Noël of ENSCP, Paris and M.P. Hehlen of M. Lujan Jr. of Neutron Scattering Centre, Los Alamos National Laboratory, Los Alamos, New Mexico.

*Note*: Octahedral coordination of the $`[YbBr_6]^3`$ unit, instead of the actual $`C_3`$ symmetry, is assumed throughout the text with no impact on our conclusions.