text:
# Temperature Dependence of the Spin Polarization in the Fractional Quantum Hall Effects

\[

## Abstract

Using a Hamiltonian formulation of Composite Fermions that I recently developed with R. Shankar, I compute the dependence of the spin polarization on the temperature for the translationally invariant fractional quantum Hall states at $`\nu =1/3`$ and $`\nu =2/5`$. I compare my results to experiments at $`\nu =1/3`$, and find reasonably good agreement.

\]

The fractional quantum Hall (FQH) effect has introduced us to new, highly correlated, incompressible states of electrons in high magnetic fields. A unified understanding of all fractions $`\nu =p/(2sp+1)`$ was achieved by the Composite Fermion picture of Jain, in which the electrons are dressed by $`2s`$ units of statistical flux to form Composite Fermions (CFs). At a mean field level, the CFs see a reduced field $`B^{}=B/(2sp+1)`$, in which they fill $`p`$ CF-Landau levels (CF-LLs), and exhibit the integer quantum Hall effect.

Due to the small $`g`$ factor of electrons in GaAs, spins may not be fully polarized in FQH states. Transitions between singlet, partially polarized, and fully polarized states (based on gap measurements) have been observed for a number of fillings, which can be understood in terms of CF’s with a spin. The transitions happen when an unoccupied CF-LL of one spin crosses the occupied CF-LL of the opposite spin.

While these low temperature measurements are in satisfactory agreement with the ground states predicted in the CF picture, in order to understand the temperature dependence of the polarization $`P(T)`$ one has to consider all excited states as well. Detailed measurements of $`P(T)`$ for the $`\nu =1/3`$ state have recently appeared in the literature. It is well-known that the $`\nu =1/3`$ state is spontaneously polarized at $`T=0`$, even when the Zeeman coupling $`E_Z=g\mu B_{tot}`$ is zero. In this it is analogous to the $`\nu =1`$ state, which has been extensively studied theoretically and experimentally. There are, however, significant differences between the two cases at finite $`T`$. The $`P(T)`$ curve for $`\nu =1/3`$ has a different shape, and has empirically been fitted to a noninteracting form $`P(T)=\mathrm{tanh}\left(\mathrm{\Delta }/4k_BT\right)`$ where $`\mathrm{\Delta }`$ is found to approximately twice $`E_Z`$.

In a recent paper, MacDonald and Palacios identified a key qualitative feature that makes $`\nu =1/3`$ very different from $`\nu =1`$. In the $`\nu =1`$ case the particle-hole excitations are very high in energy compared to $`E_Z`$, and are frozen out at all low temperatures of interest. Consequently, the $`T`$ dependence of $`P`$ comes mainly from spin wave excitations and their interactions. This is the reason why long-wavelength effective theories such as the continuum quantum ferromagnet approach are successful. However, for $`\nu =1/3`$, particle-hole excitations are on the same scale of energy as $`E_Z`$, and cannot be ignored at any $`T`$. MacDonald and Palacios use a simplified model to illustrate this feature, but the model is not sufficiently detailed to enable a calculation of $`P(T)`$ for a realistic sample. Also, the model cannot be readily extended to non-Laughlin fractions.

The goal of this paper is to describe a general analytical method for approximately computing $`P(T)`$ for an arbitrary principal fraction for realistic samples. In the last few years, R. Shankar and I have developed a Hamiltonian formalism which allows us to carry out approximate calculations for any physical quantity in the fractional Hall regime. Our central result is a formula for the LLL-projected electronic charge density at small $`q`$:

$$\rho _e(q)=\frac{\underset{j}{}e^{iqx_j}}{2ps+1}il^2(\underset{j}{}(q\times \mathrm{\Pi }_j)e^{iqx_j})$$
(1)
where $`\stackrel{}{x}_j`$ is a CF coordinate, $`l=1/\sqrt{eB}`$ is the magnetic length, and $`\stackrel{}{\mathrm{\Pi }}_j=\stackrel{}{P}_j+e\stackrel{}{A}^{}(r_j)`$ is the velocity operator of the CFs. The low-energy Hamiltonian is $`H=\frac{1}{2}\frac{d^2q}{(2\pi )^2}v(q)\overline{\rho }(q)\overline{\rho }(q)`$ where $`v(q)`$ is the electron-electron interaction. To include the effects of finite sample thickness, and to stay within the limitations of our small-$`q`$ approach, we work with a modified Coulomb interaction of the form $`v(q)=e^{\lambda q}2\pi e^2/q`$, where the length $`\lambda `$ is connected to the thickness. A notable feature of the formalism is that energy dispersions arise entirely from interaction effects. The Hartree-Fock (HF) approximation has been applied to the above Hamiltonian, and reasonable success has been obtained in computing gaps (and scaling relations between gaps) for the principal fractions, and magnetoexciton dispersions. Most recently Shankar has computed $`P(T)`$ for the compressible half-filled LL. The reason HF works so well is that our Hamiltonian is expressed directly in terms of CF coordinates. The CF’s have the right (fractional) charge and dipole moment in our formalism, and corrections to HF are expected to be small for most physical quantities.

We will compute $`P(T)`$ for $`1/3`$ and $`2/5`$ by simply carrying out the Hartree-Fock (HF) approximation for CFs at finite $`T`$. We will use a CF-LL cutoff to ensure that the correct number of electronic states occur in the Hilbert space. This implies that one must keep 3 CF-LLs for $`\nu =1/3`$ and 5 CF-LLs for $`\nu =2/5`$. While this restriction is relatively unimportant at very low $`T`$, it becomes increasingly important as the temperature increases, and occupations of the excited CF-LLs become significant.

Let us proceed to the results. We first consider the 10W sample of Khandelwal et al. The sample parameters are $`B_{}=9.61T`$, $`B_{tot}=12T`$, and each quantum well has a thickness of 260$`\AA `$. This implies that the Coulomb energy scale is $`E_C=e^2/\epsilon l_0160K`$ and the Zeeman energy is $`E_Z=0.0175E_C`$. The nominal thickness of the sample is $`\lambda 3l_0`$.

Figure 1 shows the HF prediction from our theory for $`\lambda =3l_0`$ compared to the experimental data (dashed line). The agreement is very gratifying. However, it must be regarded as fortuitous, since the simple model potential that we have assumed is unlikely to reproduce the complicated effects of density distributions in the actual sample with finite thickness. In other words, $`\lambda `$ is related to the sample thickness only in a complicated and indirect way. Another way to determine the value of the effective $`\lambda `$ is to go to the calculation of Shankar for $`\nu =1/2`$, in which he found a reasonable fit to $`P(T)`$ assuming $`\lambda 1.75l_0`$ for the same 10W sample. Accounting for the fact that the magnetic length changes when one changes filling at constant density, we estimate $`\lambda 2.2l_0`$ for $`\nu =1/3`$. Figure 1 also shows the prediction for this value (solid line). It can be seen that the predicted curve lies above the data over a range of intermediate $`T`$. This is only to be expected since the simple HF does not include the effects of spin wave excitations, or of the modification of single-particle energies due to the interaction of CFs with spin waves. The agreement between theory and experiment shown in Figure 1 is reasonably good, even for $`\lambda =2.2`$. It can also be seen that changes in $`\lambda `$ do not make huge changes in $`P(T)`$.

Figure 2 shows the same type of comparison for the data of Melinte et alfor their M242 sample. Here the sample parameters are $`B_{tot}=B_{}=17T`$, and the nominal well thickness is 250$`\AA `$. This implies that $`E_C210K`$, $`E_Z=0.0186E_C4K`$, and a nominal value for $`\lambda =4.0l_0`$. Once again this value for $`\lambda `$ is very likely an overestimate, so I have also calculated the prediction for $`\lambda =2.0l_0`$. There is a lot of scatter in the data at low $`T`$, due to the very long times needed to measure the Knight shift, and the error bars are also large at low $`T`$. This gives us some latitude in defining what we mean by the Knight shift corresponding to $`P=1`$. In any reasonable theory one expects to find that $`P=1`$ for $`TE_Z`$, and expects to see this saturated value of $`P`$ up to about $`T=0.5E_Z`$ or so.

Based on these considerations I have used two values of $`K_{s,P=1}`$, $`21kHz`$ and $`19kHz`$ both of which lie within the error bars of the low $`T`$ data. One possibility that can explain this spread is that spin-reversed quasiparticles are present in the ground state due to disorder, which can bring down the “saturated” value of the Knight shift. The $`21kHz`$ value was used by Melinte et al in a phenomenological $`\mathrm{tanh}(\mathrm{\Delta }/k_BT)`$ fit to obtain $`\mathrm{\Delta }=1.7E_Z`$. I believe that the fit for $`K_{s,P=1}=19kHz`$ is closer to the truth, since then the experimental saturation region is about $`0.5E_Z`$. The agreement between theory and experiment for this value of $`K_{s,P=1}`$ are much better than for $`K_{s,P=1}=21kHz`$. Overall the agreement is somewhat worse than for the Khandelwal et al data, but still adequate, considering the simple nature of the approximation.

Why is HF so good in this case while it was so poor for $`\nu =1`$? To answer this question let us turn to the spin wave dispersions. These can be computed in the manner described in my magnetoexciton calculation, and are shown in Figure 3 for $`\lambda =2.2l_0`$ and $`3.0l_0`$ for $`E_Z=0.0175E_C`$ and $`T=0`$. The $`q0`$ limit is required to be $`E_Z`$ by Larmor’s theorem, while the $`q\mathrm{}`$ limit is the spin-reversed particle-hole gap $`\mathrm{\Delta }_{SR}`$. Figure 3 explicitly illustrates the feature that the spin-flip particle-hole excitations are at the same energy scale as $`E_Z`$. This in turn points to the need for a more sophisticated theory of these ferromagnets which includes particle-hole excitations, spin waves, and their interactions with themselves and each other. Such a complete theory does not yet exist.

Figure 4 shows the evolution of $`\mathrm{\Delta }_{SR}`$ with $`T`$. Recall that all the energy splittings in our theory come from interactions, and as the occupations of the states change with $`T`$ so do their energies. It is clear that as $`T`$ becomes large $`\mathrm{\Delta }_{SR}`$ tends rapidly towards $`E_Z`$. This implies that the finite $`T`$ spin wave dispersion becomes increasingly flat at $`T`$ increases. Recall that a noninteracting theory would have a completely flat dispersion, that is, $`\omega (q)=E_Z`$.

Therefore, as the temperature increases, our theory becomes more weakly interacting, and our HF becomes more accurate. This trend can be expected to continue until a temperature scale when CFs cease to exist. There are no obvious signs of such a scale in the data.

Our theory is very general, and can be applied to any fractional Hall state. To ilustrate this Figure 5 shows the $`P(T)`$ curves for $`\nu =2/5`$ for $`\lambda =1.5l_0`$ for a range of Zeeman couplings. Note the nonmonotonicity of the curves that start from the singlet ground state at $`T=0`$. There is a transition to the fully polarized state around $`E_Z=0.01E_C`$. Note also that I have allowed only translationally invariant HF states, which ignores possible partially polarized states that I have recently proposed to explain intriguing observations by Kukushkin et al of a state with half the maximal polarization for $`\nu =2/5`$, which is not allowed as a translationally invariant CF state. I plan to explore the temperature dependence of the polarization, and other properties of this inhomogeneous state more thoroughly in a future publication.

Finally, let us compare our results to the only other method that can compute $`P(T)`$ for arbitrary fractions, which is exact diagonalization (keeping all the excited states) and subsequent calculation of thermodynamic quantities. Due to computational limitations, this method is restricted to fairly small systems. For example, the largest system studied by Chakraborty and Pietilainen for $`\nu =1/3`$ has 5 electrons, and for $`\nu =2/5`$ has 4 electrons. For $`\nu =1/3`$ the exact diagonalization result lies above our predictions (and the experiment) for $`T>4K`$. In fact, at $`T=0.09E_C14K`$, the exact diagonalization prediction seems to be almost a factor of two above our prediction, which essentially coincides with the experiment (Figure 1). This discrepancy might be the result of finite thickness or finite size corrections. However, at low $`T`$ the exact diagonalization result follows the data more closely than our HF approximation (in all the above comparisons I have used the $`g=0.5`$ line in Figure 2 of ref and compared to the 10W sample of Khandelwal et al. This sample has the closest parameters to those used in ref). For $`\nu =2/5`$, our results reproduce the nonmonotonicity of $`P(T)`$ for those values of $`E_Z`$ where the singlet state is the ground state, and the peaks in $`P(T)`$ occur at roughly the same $`T`$ in our results and the exact diagonalization results. However, the same overall pattern holds for $`\nu =2/5`$, namely, the results of Chakraborty and Pietilainen are below ours for low $`T`$, but are higher for $`T>0.02E_C`$, where they once again see a $`1/T`$ tail with a large coefficient. It would be interesting to explore the finite size systematics to see if the large $`T`$ tail is suppressed for larger sizes.

In summary, I have presented an approximate analytical method for computing the temperature dependence of the polarization for an arbitrary fractional quantum Hall state. An important open problem is the development of a formalism that can successfully compute corrections to HF by including particle-hole excitations, spin waves, and their interactions in a self-consistent way at finite $`T`$. While the problem is not pressing for $`\nu =1`$, where large-$`N`$ treatments give quite good agreement with experiment, it is sorely needed for fractional Hall ferromagnets. It will be interesting to see whether such corrections to HF can improve the agreement between our predictions and experiment.

It is a pleasure to thank R. Shankar, J. K. Jain, S. E. Barrett, and N. Freytag for helpful conversations, and the latter two for sharing their raw data.


format: text

text:
# Spectroscopy of Matter Near Criticality

\[

## Abstract

We propose that the finite-frequency susceptibility of matter near a class of zero-temperature phase transition exhibits distinctive excitonic structure similar to meson resonances.

\]

In this paper we predict a new spectroscopic effect that should occur very generally at quantum phase transitions described by $`O(n)`$ $`\sigma `$-models

$$=|_\nu \stackrel{}{\sigma }|^2+\mu |\stackrel{}{\sigma }|^2\lambda |\stackrel{}{\sigma }|^4.$$
(1)
These are thought to describe any continuous phase transition from a fully gapped quantum disordered insulator to a state with continuous broken symmetry, such as an antiferromagnet or a superconductor. The effect we predict, shown in Fig. 1, is a sequence of resonances in the susceptibility

$$\chi _q(\omega )=\stackrel{}{\sigma }(\stackrel{}{r},t)\stackrel{}{\sigma }(0,0)e^{i(qr\omega t)}𝑑r𝑑t,$$
(2)
when the system is slightly detuned from criticality. On the insulating side these may be understood physically as bound states of excitons. On the ordered side they may be understood as bound states of Goldstones.

This effect is important because it is formally similar to a meson spectrum in particle physics, and thus relevant to the larger issue of whether behavior like that of the empty vacuum of space might occur as an emergent phenomenon in ordinary matter. It has been known since the 1970s that matter undergoing a continuous zero-temperature phase transition is a serious candidate for this, and there is widespread agreement that the low-energy properties of such matter should be described by simple relativistic field theories . However, little is known beyond these basic facts. First-principles theory is impossibly difficult except in a few limiting cases, and all well-characterized quantum phase transitions presently known in the laboratory are either first-order, such as the solidification of 4He under pressure, or disorder-dominated. Thus the entire subject of matter undergoing a continuous quantum phase transition in the absence of disorder is largely unexplored, notwithstanding its larger implications .

Our prediction rests on a physical analogy between critical systems and 1-dimensional spin chains. The standard O(3) $`\sigma `$-model has a quantum phase transition as a function of $`\mu `$. The susceptibility at the phase transition has been calculated by $`ϵ`$-expansion by various groups  and found to have the relativistic form $`q^{\eta 2}`$, or

$$\chi _q(\omega )\frac{1}{[v^2q^2\omega ^2]^{1\eta /2}},$$
(3)
where $`v`$ is the speed characterizing the Goldstone mode on the ordered side of the transition and $`\eta 0.031\pm 0.022`$ is the anomalous exponent . Half-integral spin chains not at a phase transition, on the other hand, have been calculated by exact diagonalization methods to have the strikingly similar susceptibility

$$\chi _q(\omega )\frac{1}{[v^2(q\pi )^2\omega ^2]^{1/2}},$$
(4)
where $`v`$ is a parameter characteristic of antiferromagnetic exchange. As shown in Fig. 2, in either case one sees a diverging branch cut with a threshold on a relativistic light line. Our idea is that that these spectra are similar because they are the physically the same thing. The spin chain spectra are broad because the lowered dimensionality effectively makes the system critical by preventing antiferromagnetism and destabilizing the system toward spin-peirels order. The higher-dimensional critical points have a ground state and excitation spectrum similar to that of the spin chains. However, it has been known since the work of Fadeev and Takhtajan  that the breadth of the spin-chain spectra is due ultimately to decay of the injected spin wave into two or more spinons \- neutral, spin-1/2 particles not postulated in the underlying equations of motion that may be separated to infinity, do not bind, and do not decay. Thus our idea is that critical points in higher dimension possess new, previously unidentified elementary excitations with integrity analogous to spinons. The pairwise binding of these away from criticality is the basis of the prediction in Fig. 1.

Spinons are particularly easy to write down and explain in the case of the Haldane-Shastry Hamiltonian

$$_{HS}=J(\frac{2\pi }{N})^2\underset{\alpha <\beta }{\overset{N}{}}\frac{\stackrel{}{S}_\alpha \stackrel{}{S}_\beta }{|z_\alpha z_\beta |^2}.$$
(5)
This choice of Hamiltonian is convenient, but not essential, as spinons are believed to be a universal feature of such models, and indeed were originally discovered in the Bethe solution of the Heisenberg chain . Here $`z_\alpha =\mathrm{exp}(i2\pi \alpha /N)`$ is a lattice site on the unit circle expressed as a complex number, and $`\stackrel{}{S}_\alpha `$ is a spin-1/2 spin operator for that site. If the number of sites $`N`$ is even, then the ground state of this Hamiltonian is the spin singlet

$$\mathrm{\Psi }(z_1,\mathrm{},z_{N/2})=\underset{j<k}{\overset{N/2}{}}(z_jz_k)^2\underset{j}{\overset{N/2}{}}z_j,$$
(6)
where $`z_j`$ denotes the location of the jth $``$ site, all other being down. The energy eigenvalue of this state is $`J\pi ^2(N+5/N)/24`$. If $`N`$ is odd, on the other hand, then the ground state is a doublet, the $``$ component of which is one of the states

$$\mathrm{\Psi }_m(z_1,\mathrm{},z_M)=\underset{\alpha }{}(z_\alpha ^{})^m\mathrm{\Phi }_\alpha (z_1,\mathrm{},z_M)$$
(7)
$$\mathrm{\Phi }_\alpha (z_1,\mathrm{},z_M)=\underset{j}{\overset{M}{}}(z_\alpha z_j)\underset{j<k}{\overset{M}{}}(z_jz_k)^2\underset{j}{\overset{M}{}}z_j,$$
(8)
where $`M=(N1)/2`$. The states $`\mathrm{\Psi }_m`$ are propagating spinons of momentum

$$q_m=N\pi /22\pi (m+1/4)/N(mod2\pi ).$$
(9)
They are eigenstates of the Hamiltonian with eigenvalue

$$E_m=J(\frac{\pi ^2}{24})(N\frac{1}{N})+\frac{J}{2}(\frac{2\pi }{N})^2m(Mm).$$
(10)
Thus the spinon may be thought of as a propagating spin-1/2 defect caused by the addition of an extra site to an even-N chain. Combining the momentum and energy expressions, one obtains the dispersion relation

$$E_q=\frac{J}{2}\left[(\frac{\pi }{2})^2q^2\right](mod\pi ).$$
(11)
This tells us that spinons are relativistic particles with Dirac points at $`q=\pm \pi /2`$. The states $`\mathrm{\Phi }_\alpha `$, which are mathematically equivalent to fractional quantum Hall quasiparticles, describe a $``$ spinon at site $`\alpha `$.

Spinons attract strongly. The spin susceptibility

$$\mathrm{Im}\chi _q(\omega )=\underset{x}{}|<x|S_{q\pi }^z|0>|^2\delta (\omega E_x)$$
(12)
of the even-N Haldane-Shastry model, where $`|0>`$ denotes the ground state, $`|x>`$ and $`E_x`$ denote an excited state and corresponding excitation energy, and

$$S_q^z=\underset{\alpha }{}z_\alpha ^{qN/2\pi }S_\alpha ^z,$$
(13)
is given exactly by

$$\mathrm{Im}\chi _q(\omega )=\frac{J}{8\pi }_{\pi /2}^{\pi /2}_{\pi /2}^{\pi /2}\frac{|kk^{}|}{\sqrt{E_kE_k^{}}}$$

$$\times \delta (k+k^{}+q)\delta (E_k+E_k^{}\omega )dkdk^{}$$
(14)
for $`\pi <q<\pi `$ mod $`(2\pi )`$. The physical meaning of this formula is that spinon pairs of total momentum $`q`$ excited out of the vacuum escape to infinity but also attract each other, in that the wavefunction amplitude for the two particles to coincide is enhanced over the flat 2-spinon joint density of states by a matrix element that strongly favors low-energy decays. This effect is crucial for producing the divergence of the branch cut, for the joint density of states is flat.

The strong interaction among spinons cause their field theory to be ambiguous. Writing the Dirac sea of electrons as

$$|\mathrm{\Phi }>=\underset{|k|<\pi /2,\sigma }{}c_{k\sigma }^{}|0>(mod2\pi ),$$
(15)
where $`c_{k\sigma }=N^{1/2}_\alpha z_\alpha ^{kN/2\pi }c_{\alpha \sigma }`$, we find that the even-N ground state of Eq. (6) is exactly $`|\mathrm{\Psi }>=P_G|\mathrm{\Phi }>`$, where

$$P_G=\underset{\alpha }{}(1c_\alpha ^{}c_\alpha ^{}c_\alpha c_\alpha )$$
(16)
is the Gutzwiller projector, and that the spinon-pair eigenstates are superpositions of the wavefunctions

$$P_Gc_{\pi k}^{}c_k^{}|\mathrm{\Phi }>=P_Gc_{\pi k^{}}^{}c_k|\mathrm{\Phi }>.$$
(17)
Thus spinons have a natural fermi representation in terms of conventional particle and hole excitations of a Dirac sea or metal. However, the representation is degenerate, in that an $``$ spinon has allowed momentum only in the range $`\pi /2<k<\pi /2`$ may be written either as an $``$ particle in this range or a $``$ hole at $`\pi k`$. They are thus like Kogut-Susskind fermions, except without the “doubling problem” .

In order to talk concretely about spinon-like excitations potentially present at higher dimensional critical points let us now identify a specific transition that the $`\sigma `$-model might describe. We propose the Hamiltonian

$$=\underset{jk\sigma }{}t_{jk}c_{j\sigma }^{}c_{k\sigma }+U\underset{j\sigma \sigma ^{}}{}c_{j\sigma }^{}c_{j\sigma ^{}}^{}c_{j\sigma ^{}}c_{j\sigma }$$
(18)
$$t_{jk}=|t_{jk}|\mathrm{exp}(i_j^k\stackrel{}{A}𝑑\stackrel{}{s})_{\mathrm{plaquette}}\stackrel{}{A}𝑑\stackrel{}{s}=\pi $$
(19)
borrowed from cuprate superconductor theory . This describes electrons moving on a square planar lattice in the presence of a magnetic flux $`\pi `$ per plaquette, near-neighbor and second-neighbor integrals $`t`$ and $`t^{}`$, and an on-site coulomb repulsion $`U`$. When $`U`$ is small this system is a conventional insulator with band structure

$$E_k=\pm 2t\sqrt{\mathrm{cos}^2(k_x)+\mathrm{cos}^2(k_y)+m^2},$$
(20)
in Landau gauge $`(\stackrel{}{A}=\pi y\widehat{x})`$ with $`m=m_0\mathrm{sin}(k_x)\mathrm{sin}(k_y)`$, $`m_0=2t^{}/t`$. When $`U`$ is large it is an ordered antiferromagnet with Heisenberg exchange. The magnetic field in this model is simply a device for constructing a band insulator with one electron per unit cell, and thus enabling the transition to be to a half-integral antiferromagnet.

The parameter $`\stackrel{}{\sigma }`$ in Eq. (1) corresponds to conventional spin, and $`\mu `$ corresponds roughly to $`U`$. When $`U`$ is small, the ladder approximation for the the spin susceptibility is appropriate and gives

$$\chi _q(\omega )=\frac{\chi _q^0(\omega )}{1+U\chi _q^0(\omega )},$$
(21)
$$\chi _q^0(\omega )=\frac{1}{(2\pi )^3}\mathrm{Tr}[G_k(E)G_{k+q}(E+\omega )]𝑑E𝑑k$$
(22)
$$G_k^1(E)=E2t\left[\mathrm{cos}(k_x)\alpha _x+\mathrm{cos}(k_y)\alpha _y+m\beta \right]$$
(23)
$$\alpha _x=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\alpha _y=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\beta =\left[\begin{array}{cc}0& i\\ i& 0\end{array}\right].$$
(24)
The imaginary part of $`\chi _q(\omega )`$ consists of a broad continuum with an edge at the pair-creation threshold and a $`\delta `$-function representing a spin-1 exciton bound down below this threshold by an amount depending on $`U`$. This is illustrated in Fig. 3. This exciton at momentum $`Q=(\pi ,\pi )`$ has integrity, i.e. does not decay, even when this calculation is corrected to all orders, because it is the lowest-energy excitation of the system. The energy gap at $`(\pi ,\pi )`$ decreases as $`U`$ is increased and collapses to zero at the critical value $`U_c`$. As it does so the exciton dispersion becomes relativistic. At the transition the Hartree-Fock equations become unstable to spin density wave formation, and the same ladder sum computed with the broken-symmetry electron propagator then gives three relativistic $`\delta `$-functions - two light-like spin waves and one massive amplitude mode. The mass of the latter converges to zero at the transition. Both the integrity of the spin waves and their linear dispersion survive correction to all orders by virtue of Goldstone’s theorem .

Let us now guess, by analogy with Eq. (17) that the spinons that emerge at the transition are Gutzwiller projections of conventional particles and holes of this band structure with $`m_0`$ set to zero. Note that this is not equivalent to assuming that the true electron gap collapses. The system continues to be an insulator across the transition and has no low-energy charged excitations of any kind. This makes the spinons properly relativistic, forces their momenta to add up properly to that of the corresponding spin wave, restricts their momenta to half the Brillouin zone, and causes an $``$ particle at momentum k and a $``$ hole a momentum $`(\pi ,\pi )k`$ to describe the same spinon.

Let us now simulate the attraction between spinons by re-evaluating Eq. (21) with $`m_0`$ set to zero and $`U`$ tuned to criticality. This is done in Fig. 4. It may be seen that one properly obtains a continuum with a divergent branch cut on the relativistic light line. As is the case in 1 dimension, the spectrum is a continuum because the spinon pairs always escape to infinity. The threshold is sharp because lower-energy decays are kinematically forbidden. The edge is enhanced because the attractive force between spinons enhances the wavefunction when the particles come together with small relative momentum. But the enhancement is divergent in this case because of criticality. At the phase transition an arbitrarily small perturbation must be able to push the system into the ordered state, i.e. cause the spinons to bind. This forces the susceptibility to be infinite. The exponent of the divergence is 1 in this calculation rather than the correct $`2\eta `$, but this simply means that spinons interact by means of a potential with a long-range tail, an effect we have verified numerically.

Let us now simulate the effect of pushing to the disordered side of the transition by increasing $`m_0`$ slightly keeping $`U`$ fixed. As shown in Fig. 4, this opens a gap in the continuum and allows for a single bound state, which is the exciton in Fig. 3. The presence of this continuum and its matching the critical susceptibility at high frequencies is required by the principle that a fast measurement should not be able to distinguish a critical system from an almost-critical one. The gap in the continuum at low frequencies is required for the bound state to be sharp. The sharp bound state is required by the integrity of the exciton as the lowest-energy excitation of the disordered phase.

However, a more thoughtful analysis leads to the spectrum of Fig. 1. The continuum threshold in the bottom of Fig. 4 cannot be physically right because spinons are not legitimate elementary excitations of a band insulator. They must cease to exist immediately when the system is pushed away from criticality. The simplest way for this to occur is if the long-range tail of the interaction diverges so as to cause confinement. The result would be a bound-state spectrum more like that of a meson than a traditional exciton. The spectrum in Fig. 1 is calculated assuming that the interaction between the spinons is proportional to $`\mathrm{ln}(r)`$ rather than a point interaction. We specifically compute

$$\mathrm{Im}\chi _q(\omega )\underset{n}{}|\psi _n(0)|^2\delta (\omega E_n)$$
(25)
$$\frac{1}{2}^2\psi _n+\mathrm{ln}(r)\psi _n=E_n\psi _n.$$
(26)
The other obvious possibility is a string potential , expected on the ordered side of the transition, but the result in that case is not much different. In either case one sees a series of bound state resonances that merge into the continuum as the energy scale is raised. Whether these have integrity or decay is determined simply by whether they are above or below the pair-production threshold for the first resonance. These resonances maybe understood physically as bound states of excitons, but a more apt analogy would be with charmonium .

Thus if these resonances are observed experimentally it will suggest strongly that a relativistic gauge theory of some kind, most likely SU(2) in this case , can emerge at critical points.

This work was supported primarily by the NSF under grant No. DMR-9813899. Additional support was provided by NASA Collaborative Agreement NCC 2-794 and by NEDO.


text:
# TGRS OBSERVATIONS OF GAMMA-RAY LINES FROM NOVAE. I. LIMITS ON THE POSITRON ANNIHILATION LINE IN FIVE INDIVIDUAL NOVAE

## 1 Introduction

In the standard thermonuclear runaway model of classical nova outbursts the energy source is the explosive burning of H in a degenerate layer on the surface of a white dwarf, composed of H-rich material accreted from a companion contaminated by material taken up from the white dwarf by diffusion. The nuclear burning time-scale is rapid compared to the lifetimes of certain key $`\beta `$-unstable isotopes in the reaction chain, and substantial abundances of these nuclei are therefore built up. The isotopes involved may be 7Be, 13N, 14,15O, 17,18F, 22Na or 26Al, depending on the temperature and the chemical composition of the contaminating material from the white dwarf (CO-rich or ONeMg-rich). Convection in the accreted layer carries the longer-lived ($`>100`$ s) of these species to the surface before they decay. The $`\gamma `$-ray lines resulting from decays at the surface are in principle observable, and would provide a rather direct test of the model.

The $`\beta `$-unstable nuclei involved are all proton-rich and therefore decay predominantly by emission of positrons. Thus in addition to $`\gamma `$-ray lines from de-excitation of the daughter nuclei, we also expect a line at 511 keV from the annihilation of the positrons with the ambient material (Clayton & Hoyle 1974, Leising & Clayton 1987). Several previous experiments have searched for the decay lines from the relatively long-lived isotopes 22Na and 7Be, either from individual nearby novae (Leising et al. 1988, Iyudin et al. 1995) or from the accumulated production of many novae in a wide field of view (Leising et al. 1988, Harris, Leising & Share 1991, Harris et al. 1997). Here we report first results of the first search for the 511 keV line, which reflects primarily the production of shorter-lived isotopes such as 13N and 18F.

The Transient Gamma Ray Spectrometer experiment (TGRS) is very well suited to the search for the 511 keV line, for several reasons. First, it is located on board a spacecraft whose orbit is so elliptical that it has spent virtually all of its mission (1994 November–present) in interplanetary space (the WIND mission). In this environment the $`\gamma `$-ray background level is relatively low, and interruptions due to passages of Earth’s trapped radiation belts are minimal (typically lasting a few hours at intervals of several months). This is important because the bursts of 511 keV radiation are predicted to occur only during a $`6`$ hr period at an uncertain interval before the nova event (Leising & Clayton 1987, Hernanz et al. 1997, and §2). Second, the TGRS instrument is attached to the south-facing surface of the rotating cylindrical WIND body, which points permanently toward the south ecliptic pole. The detector is unshielded, and TGRS therefore has an unobstructed view of the entire southern (ecliptic) hemisphere. Third, and most importantly, the TGRS Ge detector has sufficient spectral resolution to detect a 511 keV line which is slightly Doppler-shifted away from the background 511 keV line which is always present in space experiments, which arises from annihilation following $`\beta ^+`$ decays of cosmic-ray spallation product nuclei. This background line is always very close to the rest energy, whereas the most recent nova models predict the source line to be blueshifted by 2–5 keV (Hernanz 1997, private communication); by comparison, the TGRS energy resolution at 511 keV has varied from 3–4 keV FWHM during the mission (Harris et al. 1998). No previous experiment has possessed the same advantages of long and continuous temporal coverage, broad spatial coverage, and fine spectral resolution.

Our analysis procedure (described in §3) relies heavily on the theoretically predicted properties of the 511 keV line. There has not been much incentive for theoretical work upon this line because of the difficulty of resolving the nova line from the background line in the previous generation of low-resolution experiments. The early work of Clayton & Hoyle (1974) and Leising & Clayton (1987) assumed white dwarf models which were CO-rich, and treated the explosion parametrically. Much fuller hydrodynamic models which can be applied at all stages from accretion through explosion and nucleosynthesis have recently been developed (Starrfield et al. 1992, Hernanz et al. 1996, Jose et al. 1997) and applied to the positron annihilation problem (Hernanz et al. 1997, Gómez-Gomar et al. 1998). In this series of models the largest 511 keV line fluxes come from the more massive ONeMg-rich white dwarfs; indeed, fluxes as high as $`10^2`$ photon cm-2 s-2 are expected in favorable cases. For the purpose of analysis we shall in general assume the line properties of the Hernanz et al. (1997) ONe2 model, which generates this level of flux over a $`6`$ hr period (§3).222 This model is referred to as A2 in Gómez-Gomar et al. (1998). In §2 we describe the difficulties which are involved in a search for the 511 keV line around the times of those novae which are known to have been in TGRS’s field of view during the mission. Our results for the individual novae are presented in §4, and we will show in §5 that for the ONe2 and other nova models TGRS sensitivity is good enough to justify a large-scale search for the line in the full 1995–1997 TGRS data.

## 2 Sample of Novae

There were five known classical novae in the southern hemisphere between 1995 January and 1997 June, whose relevant properties are given in Table 1. It should be noted that in many cases the observational material on which they are based is extremely incomplete and their properties are often very poorly known. Columns 1 and 2 of Table 1 identify the nova. Columns 3 and 4 indicate the day of discovery and the apparent magnitude $`m_V(t=0)`$ at discovery. It is important to note that many novae are probably discovered some time after visual maximum $`m_V`$, which occurs $``$1–4 days after the explosion and the 511 keV line emission. To bracket the time during which the explosion must have occurred, we give the date of the last prediscovery observation (or upper limit) in Column 5.

Observations during the early decline have been used to estimate the speed class $`t_2`$, defined as the time necessary for $`m_V`$ to fall by two magnitudes from the discovery value $`m_V(t=0)`$ (Column 6, Table 1). The typical situation is illustrated in Fig. 1. Two facts are apparent from the Table and the figure. First, the sample of southern novae is likely to be very incomplete. Gaps of $`>20`$ d in the coverage of the known novae are common; therefore, very fast novae with $`t_2`$ of the order of a few days could rise and then fall below the detectable level during such a gap. Warner (1989) estimates that this incompleteness could affect even novae so bright as $`m_V3`$. Second, the distance to any nova can only be determined by comparing the apparent magnitude at some known epoch (usually at visual peak) with an absolute magnitude $`M_V`$. Empirical formulae are available which relate the absolute magnitude $`M_V`$ to the speed class $`t_2`$ (Warner 1995). However, the peak $`m_V`$ is unknown since the time of the peak is unknown, and it can only be estimated by extrapolating back up the prediscovery light curve (dashed line in Fig. 1). Thus the estimated distance $`r`$ to the nova is a function $`r(t)`$ of the time before discovery at which the outburst is assumed to occur.

The relationships involved are the well-known distance modulus formula including extinction $`A_V`$

$$m_VM_V=5\mathrm{log}r5+A_V$$
(1)
and the extrapolation up the prediscovery light curve

$$m_V(t)=m_V(t=0)+2\frac{t}{t_2}$$
(2)
which can be combined to yield a transcendental relation for the estimated distance $`r(t)`$

$$\mathrm{log}r(t)+\frac{A_V(r)}{5}=1+\frac{m_V(t)M_V}{5}$$
(3)
if $`M_V`$ and $`A_V(r)`$ are given. We solved these equations numerically for each nova during each prediscovery interval using an empirical relation between $`M_V`$ and $`t_2`$

$`M_V`$ $`=`$ $`2.41\mathrm{log}t_210.7\text{for 5 d }<t_2<50\text{ d}`$ (4)
$`=`$ $`9\text{for }t_25\text{ d}`$ (5)
$`=`$ $`6.6\text{for }t_2>50\text{ d}`$ (6)
(Warner 1995) and the prescription for extinction as a function of celestial coordinates and of distance given by Kil’pio & Malkov (1997). The range of distances given in Column 7 of Table 1 for each nova was thus obtained. In all cases the smallest value holds true if the outburst occurred just after the last prediscovery measurement (Column 5) and the largest value corresponds to an outburst occurring at the time of discovery (Column 3). We neglect uncertainties due to the scatter in the empirical relation Eqs. (4–6), which are substantial, but still are less important than the uncertainty in the outburst epoch for relatively fast novae such as most of those in Table 1.

Very little information appears to be available on the abundances of the novae in Table 1. Woodward et al. (1998) have shown that N Cru 1996 was definitely of the neon subclass, involving an ONeMg white dwarf. More tentatively, enhancement of N but not of O in N Cir 1995 (Greeley, Blair & Long 1995) suggests that it might have occurred on a very massive ONeMg white dwarf (Politano et al. 1995). Since theoretical models (Hernanz et al. 1997, Gómez-Gomar et al. 1998) suggest that the 511 keV line is much more likely to be detectable from neon novae, both these novae might be promising candidates for detection if they are at the nearer ends of the ranges in Table 1 Column 7.

## 3 Observations and Analysis

### 3.1 Observations

Since launch in 1994 November the WIND spacecraft has been raised into a complicated and extremely elliptical orbit, most of which is spent outside Earth’s magnetosphere. Perigee passes taking the spacecraft within the magnetosphere (and therefore exposing it to high charged particle fluxes in the trapped radiation belts) are brief (typically $`1`$ d). The overall level of background $`\gamma `$-radiation is therefore relatively low, being mainly due to the cosmic diffuse background and the effects of irradiation by the Galactic cosmic rays.

The TGRS detector is a radiatively cooled 35 cm2 Ge crystal sensitive to energies between 20 keV–8 MeV. It is mounted on the top surface of the WIND spacecraft body, which permanently faces the south ecliptic pole; the spacecraft is stabilized by rotation in the ecliptic plane with a 3 s period (Owens et al. 1995). Although the detector is unshielded, some capability for localizing sources exists along the ecliptic plane, since a passive Pb occulter is fixed concentric with the detector and carried around it by the spacecraft’s rotation, thereby modulating signals from the rotation plane with a 3 s period (see figures in Teegarden et al. 1996).

The database of regular TGRS background spectra consists of count rates between 20 keV and 8 MeV which have been accumulated continuously in 1 keV energy bins during 24 min intervals since launch, except for the epochs of perigee passages. The data stream is also interrupted by memory readouts of triggered data from solar flares and $`\gamma `$-ray bursts for periods of $`2`$ hr (Owens et al. 1995). We summed these spectra into 6 hr intervals, corresponding to the expected time-scale of the signal from a classical nova (Hernanz et al. 1997). As noted in §2 above, it is frequently unclear how much time elapsed between the explosion and the visual discovery for any given nova. We were therefore usually obliged to take 6-hr spectra for the entire period prior to discovery, as far back as the latest available prediscovery observation (dashed line, Fig. 1). Of course, the implied apparent magnitude of the nova brightens as one proceeds back along the light curve, and in some cases this in turn implies that a very bright nova must somehow have escaped detection. We therefore did not consider epochs prior to discovery for which $`m_V`$, extrapolated according to speed class and corrected for extinction, would have fallen below 3.0 (cf. Warner 1989).

Two of the novae listed in Table 1 lay in the approximate direction of the Galactic center, near where the Galactic plane crosses the ecliptic, and so were occulted by the Pb occulter (N Sgr 1996 and N Sco 1997). The TGRS occulted data are returned in four broad 64-channel energy windows, only one of which is binned at 1-keV resolution; this window covers the region 479–543 keV and therefore contains the blueshifted nova line predicted by Leising & Clayton (1987) and Hernanz et al. (1997). For N Sgr 1996 and N Sco 1997 we therefore summed the occulted spectra over the same prediscovery time intervals as described above for background spectra. We thus expected to have an independent check (by means of location) on any positive signal discovered in the background spectra.

### 3.2 Analysis of Background Spectra

A characteristic TGRS background spectrum in the range 490–530 keV during a 6-hr interval is shown in Fig. 2. The overall features of such spectra, which are very stable from one spectrum to the next, are the following: first, a strong line due to positron annihilation in the spacecraft and the passive material surrounding the detector. This line is expected to lie at the exact rest energy 511.00 keV333 The instrument gain has varied slightly (by $`0.5`$%) during the mission. For the purpose of measuring line positions close to 511 keV, we can use the position of this background annihilation line as a reference energy to correct this gain shift. The procedure involved — fitting a Gaussian line and equating the position of its peak to 511.00 keV — is described by Harris et al. (1998).. The line shows a clear asymmetry in the form of a broad red wing of instrumental origin. This asymmetry is probably caused by charge collection losses; the increasing effects of radiation damage have caused it to become more pronounced as the mission has proceeded. Second, there is an underlying continuum at energies above and below the background line energy, which over a short range of spectrum may be approximated by a power law. Third, there is a discontinuity in the continuum at the line energy which is due to incomplete absorption of the line photons in the detector.

The principle behind our analysis was to fit the selected background spectra (see §3.1) with such components as those described above, plus lines fixed at the energies and widths predicted by theory from neon novae. Significant fitted amplitudes in these lines would be taken as detections. The predicted lines are described in Table 2; they correspond to the same annihilation line at two different epochs, 6 and 12 hours after the explosion \[Hernanz et al. 1997, Hernanz 1997 (private communication), Gómez-Gomar et al. 1998\]. We will refer to these as the ”6-hr” and ”12-hr” lines respectively. Note that, while the 6-hr line is quite well separated from the background line, the 12-hr line is only blueshifted from it by 2 keV and will thus be badly blended, even with TGRS’s superior energy resolution444As measured by Harris et al. (1998) the FWHM energy resolution of TGRS at energy 511 keV degraded during the mission from 3.2 keV to $`>4`$ keV due to radiation damage from cosmic-ray impacts; this loss of resolution appears to be due to the growth of red wings on spectral lines as described below.. We thus expect worse statistical and also systematic errors on the 12-hr line. The detection and measurement of the nova line signal (if any) thus depend almost entirely on our results for the 6-hr line.

Fitting the full spectrum shown in Fig. 2 with the full complement of model components (background 511 keV Gaussian line; nova 6-hr or 12-hr Gaussian line555 The fits are not strongly sensitive to the precise line shape, so long as the FWHM is fixed at 8 keV (Table 2). For example, Leising and Clayton (1987) argued in favor of the broad flat-topped line shape expected from the near side of an expanding shell. We tested this line shape and found little difference from results obtained with the Gaussian shape which Hernanz (1997, private communication) used as an approximation.; power law; step function at 511 keV) does not lead to acceptable values of $`\chi ^2`$ (e.g. $`500`$ for 35 degrees of freedom in Fig. 2). The main reason for this is the asymmetric shape on the red side of the background 511 keV line which is apparent in Fig. 2, whose origin is unclear. It became more pronounced as the mission proceeded, and was therefore probably due to accumulated radiation damage; an investigation of this effect is in progress (Kurczynski et al. 1999).

Although the anomalous red component of the background 511 keV line is not perfectly understood, this is not of crucial importance to our search for excess emission on the blue side of the profile. We have found that the best fits are obtained when a simpler model (background 511 keV line; nova 6-hr or 12-hr line; constant approximation to power law) is fitted to energies $`511.0`$ keV only (Fig. 3a, b). The values of $`\chi ^2`$ for such fits typically lie between 15–25 for 17 degrees of freedom.

In detail, our method of fitting each spectrum was to vary the model parameters until the count rate was best reproduced according to the minimum of $`\chi ^2`$; errors were calculated by the method of Lampton, Margon & Bowyer (1976) for mapping the range of parameter space where the minimum $`\chi ^2`$ value is exceeded by 1. The line amplitudes from the count rate were normalized by the detector effective area (a function of the direction to the nova), which is known from Monte Carlo simulations (Seifert et al. 1997). In Fig. 3b (dotted line) we use this detector response to show that a neon nova with the theoretically expected properties would be easily detectable at a distance 1 kpc.

A typical series of fitted amplitudes for the 6-hr line, during times when N Sco 1997 may have erupted, is shown in Fig. 4a. Three features of this result call for attention. First, the points are distributed about a mean value in a way consistent with random scatter ($`\chi ^2`$ per degree of freedom = 0.71). Second, (consistent with this) there is no individual point of high significance relative to the mean such as would be expected from a real nova line. Third, the mean value is not zero — there is a systematic positive nova line amplitude.

The systematic positive amplitude is probably due to a systematic departure of the background 511 keV line from Gaussian shape on its blue wing (for example, if the anomalous red component seen in Fig. 2 has a Gaussian shape, its blue tail would effectively contribute such a constant value to the nova line). We have found in all cases that the systematic positive offset is very stable. It has increased slowly during the mission; the time-scale of the increase ($``$ months to years) is much longer than the time scales ($`10`$ days) over which we search for lines from a given nova. In this respect the behavior of the positive offset is very similar to that of other measures of degradation of detector performance, such as energy resolution (Harris et al. 1998) and the amplitude of the anomalous red component seen in Fig. 2.

Apart from this stable and well-behaved offset (which can be removed from the nova line flux by simple subtraction) there appear to be few or no systematics in our nova line measurements made from background spectra. We therefore believe that single significant positive measurements (following subtraction of the systematic offset) would be valid evidence for short-lived ($`6`$ hr) line emission. We now turn to two internal checks which could be made on such a detection. In §3.3 we examine the usefulness of data modulated by the TGRS occulter for some celestial directions. In §3.4 we show that small improvements in sensitivity can be obtained from measuring the 511 keV line at epochs subsequent to its peak flux at $`t=6`$ hr (what we have called the ”12-hr line”).

### 3.3 Analysis of Occulted Spectra

The Pb occulter referred to in §3.1 subtends an angle of $`90^{}`$ along the ecliptic and approximately $`16^{}`$ across it, and is 1 cm in thickness. The spectra in 64 1-keV energy channels around 511 keV, when modulated with a 3 s period by the occulter, are binned into 128 angular channels (”sectors”) of $`2.8125^{}`$ in ecliptic longitude. We obtained spectra for N Sgr 1996 and N Sco 1997 by fitting occultation dips centered on the respective ecliptic longitudes ($`\lambda =275.66^{}`$ and $`268.73^{}`$) to each energy channel (cf. Teegarden et al. 1996, Harris et al. 1998). A typical spectrum for a 6 hr period is shown in Fig. 5.

All background features varying on time-scales $`>3`$ s are subtracted out of such a spectrum, notably the 511 keV background line which is so prominent in Fig. 2. The subtraction is very effective and the residual spectrum in Fig. 5 is almost featureless. For occultations having phase close to that of the Galactic center ($`\lambda =266.8^{}`$), the strongest feature in the residual spectrum is expected to be the diffuse Galactic 511 keV annihilation line (see e.g. Purcell et al. 1997). However the measurement by Harris et al. (1998) with the TGRS occulter during 90 d intervals indicates that during 6 hr intervals this line will be detected at a level of only $`0.35\sigma `$.

We therefore fitted the occulted spectra between 490–530 keV by a model containing two lines, one at exactly 511 keV due to the diffuse Galactic source, and one at either the 6-hr or 12-hr line energy (Table 2). The gain-corrected 511 keV line positions were taken from the corresponding background spectra (see footnote, §3.2). The 511 keV line widths were fixed at the values measured by Harris et al. (1998) during 90 d intervals containing the nova epoch.

Results from a single 6 hr period for the 6-hr line in N Sco 1997 are shown in Fig. 5 (solid, dashed and dot-dashed lines); in Fig. 6 we show the time series of all such fits for this nova. Although there do not appear to be any large systematic effects, it is immediately apparent that some of the 6 hr intervals exhibit anomalously large errors. These situations occur in a model which possesses multiple minima of $`\chi ^2`$ lying close to each other in the model’s parameter space. The region defined by $`\chi ^2+1`$ is then much larger than would naively be expected from the behavior of $`\chi ^2`$ in the region around the true minimum. In the present case, which is a two-line model where the lines are not well resolved, such situations tend to arise when the spectrum is best fitted by lines with amplitudes of opposite signs. We found this to be true of the badly-behaved fits in Fig. 6.

### 3.4 Combination of 6-hr Line and 12-hr Line Results

The time series of background spectrum fit results for N Sco 1997 are shown in Fig. 4a for the 6-hr line and in Fig. 4b for the 12-hr line. We see that the statistical errors on the 12-hr line are generally much larger, as expected, due to worse blending with the 511 keV background line. The systematic positive offset described in §3.2 for the 6-hr line is also worse in the 12-hr line flux measurement.

The 6-hr and 12-hr line fits are independent and can be combined so as to improve the significance of a detection in one of them. In our simple approximation, they are combined with a weight proportional to their contribution to the Hernanz et al. (1997) light-curve (c.f. Table 2) and proportional to the inverse square of the statistical error, in the usual way. Each 6-hr line measurement (Fig. 4a) is combined with the 12-hr line measurement from the following 6-hr interval (Fig. 4b). The results for N Sco 1997 are shown in Fig. 4c for background spectrum fits666In this figure, and from here onward, the positive offset in these measurements is subtracted off.

It is clear from Fig. 4b that the quality of the 12-hr line measurements is much worse than that of the 6-hr lines, as expected. Thus the 12-hr line measurements yield only a small gain in the sensitivity of the line search, when combined with the 6-hr line measurements which contain most of the information. This may be illustrated by the similarity between Figs. 4a and 4c; clearly very little information was added by including the 12-hr line.

## 4 Results

Our results for N Cir 1995, N Cen 1995, N Sgr 1996, and N Cru 1996 are shown in Figs. 7–10. For N Sgr 1996 there are both background and occulted spectra (Figs. 9a and 9b), for the others only background spectra. The results are given as the nova line flux at all epochs when the thermonuclear outburst could have occurred; note that (as derived in §2) the theoretical value with which they are compared is time-dependent (dotted line). The positive systematic offsets described in §3.2 have been subtracted. In the absence of any single strong positive 6-hr line detections, we have combined the 6-hr and 12-hr line results as described in §3.4.

Once the positive offset has been removed, the distribution of the measured line intensities is quite free from systematics, as previously deduced from the N Sco 1997 results (§3.2 and Fig. 4c). The distribution of the complete ensemble of measurements (264 points) agrees very well with a random distribution about zero. We therefore conclude that none of the novae in our sample emitted the expected blueshifted 511 keV line, down to a limiting flux which varied for each nova. We express this limit conservatively as the $`3\sigma `$ upper limit from the statistical errors on each 6-hr point; it is represented by dashed lines in Figs. 4c and 7–10, and the mean values for each nova are given in Table 3. The mean values vary slightly from nova to nova partly because of the slow degradation of detector performance (§3.2), but mainly due to the variation of effective area with angle of incidence of the nova in the detector.

The typical $`3\sigma `$ sensitivities in Table 3 may be compared with the prediction of $`1.6\times 10^2`$ photon cm-2 s-1 for a neon nova at 1 kpc (Table 2). From this we derive our most important conclusion, namely that the TGRS measurements are sensitive enough to detect neon novae to a distance of about 2.5 kpc at $`3\sigma `$ significance, for favorable estimates of the nova line emission (Hernanz et al. 1997, model ONe2).

There remains a small possibility that a 511 keV line event associated with one of the novae has been missed due to incomplete coverage. It can be seen from Figs. 4 and 6–10 that losses of entire 6-hr intervals are rare, generally involving perigee passes of $`1`$ d (Figs. 7, 8, 10). Briefer losses of data within a few 6-hr intervals cause larger than normal error bars, and arise when TGRS’s operation mode changes following a $`\gamma `$-ray burst trigger, or from telemetry errors. Overall, we estimate that the instrument’s temporal coverage at the characteristic sensitivities in Table 3 (2–3 $`\times 10^3`$ photon cm-2 s-1 in 6 hr at $`3\sigma `$) is 80–85%, and that $`97`$% of all 6 hr intervals contain at least some live time.

## 5 Discussion

### 5.1 Results for Individual Novae

As emphasized in §2, it is very difficult to use our results for the individual novae to constrain the ONe2 model because of our lack of knowledge of two related nova parameters, the epoch of explosion and the distance. In Figs. 4c and 7–10 we use the limited information given in Table 2 to plot the predicted fluxes $`\varphi _{pred}`$ as dotted lines, which depend on the time relative to discovery according to

$$\varphi _{pred}(t)=0.016/r(t)^2$$
(7)
where the distance $`r(t)`$ in kpc is obtained from Equation (3). Note that $`\varphi _{pred}`$ is not a light curve; as shown in Table 1, there is usually a wide range of allowed explosion times, and $`\varphi _{pred}`$ is the predicted flux at any one time of explosion. In two cases (N Sco 1997 and N Sgr 1996; Figs. 4c and 9a) the predicted flux is always below our $`3\sigma `$ upper limit, for any of the allowed explosion times; these novae must be too distant for our results to be used to constrain even the most optimistic nova model (Hernanz et al. 1997, model ONe2). The other three novae (N Cir 1995, N Cen 1995 and N Cru 1996) are predicted to yield a detectable $`\varphi _{pred}`$ only if the explosion occurred early in the allowed range of times \[prior to $`12`$ d before discovery for N Cir (Fig. 7), 7 d before discovery for N Cen (Fig. 8), and 9 d before discovery for N Cru (Fig. 10)\]. In other words, the allowed range of distances for these novae includes the TGRS detection limit $`2.5`$ kpc. The ONe2 model can be tested using our results for these novae if further investigations yield either the appropriate constraints on explosion time, or an independent measurement of the distance to the nova.

### 5.2 Implications for Nova Searches

The data analyzed in this paper covers about 70 d, representing only a small fraction of the whole TGRS database. An obvious extension of this work would be to analyze the entire database in order to detect, or place upper limits upon, the occurrence of all neon novae lying within the 2.5 kpc radius estimated in §4. There are good reasons for believing that many novae within this radius escape detection. As noted in §2, visual searches for novae are usually episodic, containing gaps in temporal coverage within which fast-rising and decaying novae are likely to be missed; coverage is particularly uneven in the southern hemisphere. Further, the southern Galactic plane contains regions of very strong interstellar absorption; $`A_V`$ values up to 4 are found within 2.5 kpc in some directions (Neckel & Klare 1980). There is therefore considerable uncertainty, not only about the global Galactic nova rate, but also about its spatial distribution; novae in a disk population may be more common than the spheroidal population modeled by Leising & Clayton (1985), but less observable due to obscuration (Hatano et al. 1997). There is evidence that the neon subclass may be more concentrated toward the plane, and thus more under-represented in global nova statistics (Della Valle et al. 1992).

Estimates of the success rate of a $`\gamma `$-ray line search must be very uncertain, for the above reasons; in fact, the best justification for such a search is precisely that it is the only way to obtain an unbiased census of novae on which to base global nova rates. Estimated rates for the Galaxy as a whole range from 11–260 yr-1, tending to cluster about values 35–50 (Shafter 1997), with about one-third occurring on ONeMg white dwarfs (Gehrz et al. 1998). To estimate TGRS’s rate of detections, we note that a full-scale search of almost 3 years of TGRS data will yield a much higher chance detection rate than that found in the present work ($`<1`$ detection in 70 d), so that the threshold for detection must be set higher than the value $`3\sigma `$ used here. Since the distribution of our null measurements is very close to the normal distribution expected by chance, it is easy to show that a threshold level $`4.6\sigma `$ yields a probability $`<1`$% of a single false detection by chance (Abramowitz & Stegun 1964). This threshold level corresponds to a distance $`2`$ kpc for a neon nova. A simple model of nova distribution following star counts (Bahcall 1986, integrating the disk over $`z`$) suggests that $`2`$% of Galactic novae occur within 2 kpc of the Sun in the southern hemisphere. Moderate estimates of the global nova rate $`50`$ yr-1 then imply that TGRS would detect one neon nova during $`3`$ yr; higher global rates, or a higher fraction of neon novae, would imply more detections.

Further papers in this series will utilize this search strategy, and will follow up the possibilities arising from detailed study of the already known novae as described in the previous section.

We acknowledge the work of the small but dedicated teams who monitor southern hemisphere novae: the Royal New Zealand Astronomical Society (communicated by F. Bateson) and the VSNET network of variable star observers (communicated by T. Kato). We are grateful to Dr. M. Hernanz for providing pre-publication results and to Drs. C. Woodward and M. Greenhouse for helpful discussions. Theresa Sheets (LHEA) and Sandhia Bansal (HSTX) provided assistance with the analysis software.


text:
# Quasispecies evolution on a fitness landscape with a fluctuating peak

## I Introduction

Ever since Eigen’s work on replicating molecules in 1971 , the quasispecies concept has proven to be a very fruitful way of modeling the fundamental behavior of evolution. A quasispecies is an equilibrium distribution of closely related gene sequences, localized around one or a few sequences of high fitness. The kinetics of these simple systems has been studied in great detail, and the formulation has allowed many of the techniques of statistical physics to be applied to replicator and evolutionary systems, see for instance .

The appearance in these models of an error-threshold (or error-catastrophe) as an upper bound on the mutation rate above which no effective selection can occur, has important implications for biological systems. In particular it places limits on the maintainable amounts of genetic information  which puts restrictions on theories for the origins of life.

Until now studies of quasispecies have mainly focused on static fitness landscapes. However, many organisms in nature live in a quickly changing environment. In this paper we will study how a population responds to changes in the fitness landscape. More precisely we will study the population dynamics on a fluctuating single peaked fitness landscape. Since the full theory turns out to be impossible to solve analytically, we introduce a simple approximation that makes the rate equations analytically tractable. The expression for the error threshold is then obtained from the expression in the static case by replacing the height of the fitness peak by the time average of the height of the fluctuating peak. We also study how the phase-shift between fitness oscillations and population dynamics depends on the frequency in the case of a small harmonic fluctuation.

## II Quasispecies in dynamic environments

A quasispecies consists of a population of self-replicating genomes, where each genome is represented by a sequence of bases $`s_k`$, $`\left(s_1s_2\mathrm{}s_\nu \right)`$. We assume that the bases are binary, $`s_k\{1,0\}`$ and that all sequences have equal length $`\nu `$. Every genome is then given by a binary string $`\left(011001\mathrm{}\right)`$, which also can be represented by an integer $`k=_js_j2^j`$ ($`0k<2^\nu `$).

To describe how mutations affect a population we define $`W_k^l`$ as the probability that replication of genome $`l`$ gives genome $`k`$ as offspring. We only consider point mutations, which conserve the genome length $`\nu `$.

We assume that the point mutation rate $`\mu =1q`$ (where $`q`$ is the copying accuracy per base) is constant in time and independent of the position in the genome. We can then write down an explicit expression for $`W_k^l`$ in terms of the copying fidelity:

$`W_k^l`$ $`=`$ $`\mu ^{h_{kl}}q^{\nu h_{kl}}=q^\nu \left({\displaystyle \frac{1q}{q}}\right)^{h_{kl}}`$ (1)
where $`h_{kl}`$ is the Hamming distance between genomes $`k`$ and $`l`$. The Hamming distance $`h_{kl}`$ is defined as the number of positions where the genomes $`k`$ and $`l`$ differ.

The equations that describe the dynamics of the population take a relatively simple form. Let $`x_k`$ denote the relative concentration and $`A_k(t)`$ the time-dependent fitness of genome $`k`$. We then obtain the following rate equations:

$`\dot{x}_k(t)`$ $`=`$ $`{\displaystyle \underset{l}{}}W_k^lA_l(t)x_l(t)e(t)x_k(t)`$ (2)
where $`e(t)=_lA_l(t)x_l(t)`$, and the dot denotes a time derivative. The second term ensures the total normalization of the population ($`_lx_l(t)=1`$) so that $`x_k(t)`$ describe relative concentrations.

In the classical theory introduced by Eigen and coworkers , the fitness landscape is static. The rate equations (2) can then be solved analytically by introducing a change of coordinates that makes them linear, and then solving the eigenvalue system for the matrix $`W_k^lA_l`$. The equilibrium distribution is given by the eigenvector corresponding to the largest eigenvalue.

If the fitness landscape is time-dependent, this method cannot be applied. A time-ordering problem occurs when we define exponentials of time-dependent matrices, since in general the matrix $`W_k^lA_l(t)`$ does not commute with itself at different points in time. Later in this paper we make a simple approximation that makes the rate equations one-dimensional; time-ordering is then no longer necessary.

Much of the work on quasispecies has focused on fitness landscapes with one gene sequence (the master sequence) with superior fitness, $`\sigma `$, compared to all other sequences. These are viewed as a background with fitness $`1`$. These landscapes are referred to as single peaked landscapes. The master sequence is denoted $`x_0`$. In this paper we focus on single peaked landscapes where the height of the fitness peak is time-dependent. The fitness landscape is then given by

$`A_k(t)`$ $`=`$ $`\{\begin{array}{cc}\sigma (t)\hfill & \text{if }k=0\hfill \\ 1\hfill & \text{otherwise}\hfill \end{array}`$ (5)
This class of time-dependent landscapes was studied by Wilke and co-workers . They investigated the behavior of a periodically fluctuating single peak landscape by numerically integrating the dynamics to find the limit cycle of the concentrations for a full period.

Fig. 1 shows how the concentration of the master sequence responds to a sudden, sharp jump in its fitness. When the fitness changes it takes some time for the population to reach the new equilibrium. It is this delay that causes a phase shift between a periodically changing fitness function and the response in the concentrations. The relaxation time of the population to the appropriate equilibrium distribution depends on both the fitness values of the landscape and the mutation rate. For extremely slow and smooth changes in the fitness the population will effectively reach equilibrium at every point in time. Thus the continued existence of a quasispecies will depend on the local dynamics of the landscape. When the landscape changes quickly the population will fail to follow the changes adequately and thus responds to the landscape dynamics in a way that is typical of a low pass filter. The following section examines the fluctuating single peak landscape in some detail. In particular, we introduce an approximation that lets us find an analytic form for the relaxation time of the population, and the phase lag it introduces in a periodic landscape.

## III Approximate quasispecies dynamics

We now introduce a simple approximation of the model presented above. In this approximation we can solve the rate equations and find an expression for the concentration of the master sequence $`x_0(t)`$. In the limit of long chain-length ($`\nu 1`$) we can neglect back-mutations from the background to the master sequence. This gives a simplified one-dimensional version of the rate equation of the following form:

$`\dot{x}_0(t)`$ $`=`$ $`Q\sigma (t)x_0(t)e(t)x_0(t)`$ (6)
where $`Q=q^\nu `$ is the copying fidelity of the whole genome and $`e(t)=(\sigma (t)1)x_0(t)+1`$.

Fig. 2 compares the concentration of the master sequence calculated by solving approximation 6 and by numerically integrating the full rate equation 2. The figure shows that the approximation is quite accurate.

Since this equation is one-dimensional there is no time-ordering problem and it can be solved analytically for non-periodic peak fluctuations. Equation 6 can be transformed to a linear form by introducing a new variable $`y(t)=\frac{Qx_0(t)}{(1Q)x_0(t)}`$. This gives

$`\dot{y}(t)`$ $`=`$ $`1(Q\sigma (t)1)y(t)`$ (7)
which can be solved. Substituting back gives the concentration of the master sequence

$`x_0(t)`$ $`=`$ $`x_0{\displaystyle \frac{e^{_0^t(Q\sigma (s)1)𝑑s}}{1+x_0_0^te^{_0^s(Q\sigma (u)1)𝑑u}(\sigma (s)1)𝑑s}}`$ (8)
Since we are only interested in the long time behavior of the system we can ignore transients carrying memory from initial values. Assuming $`e^{_0^t(Q\sigma (s)1)𝑑s}1`$ gives

$`x_0(t)`$ $`=`$ $`{\displaystyle \frac{Q}{1+(1Q)_0^te^{_s^t(Q\sigma (u)1)𝑑u}𝑑s}}`$ (9)
This is a generalization of the static expression for the asymptotic concentration:

$`x_0^s`$ $`=`$ $`{\displaystyle \frac{Q\sigma 1}{\sigma 1}}`$ (10)
On a static single peaked fitness landscape there is a phase transition in the concentration distribution when the copying fidelity decreases below a critical value . At high mutation rate the selective advantage of the master sequence due to its superior fitness is no longer strong enough for the gene sequences to be localized in sequence space. Instead they diffuse over the entire sequence space, and the distribution becomes approximately uniform. This is generally referred to as the error catastrophe or error threshold and is one of the main implications of the original quasispecies model. By making the same approximation as above, i.e. assuming no back-mutations onto the master sequence, the static landscape error threshold can be shown to occur when $`Q=\frac{1}{\sigma }`$. In other words, the transition occurs when the selective advantage of the master-sequence no longer is able to compensate for the loss of offspring due to mutations. This can also be seen from Eq. 10 which defines the stationary distribution of the master sequence in the static case.

One has to be careful when discussing the error threshold on a fluctuating peak. The fitness can, for example, slowly move from being strong enough to localize the population around the peak, to beibg so weak that the population delocalizes, and then back again. If we however consider an average over a time scale much longer than the fluctuation time of the fitness peak, a sensible definition of the error threshold can be made based on the average concentration of the master sequence. The time average of the concentrations can be found by rewriting equation 6 as differentials

$`{\displaystyle \frac{dx_0}{x_0}}`$ $`=`$ $`{\displaystyle (Q\sigma (t)1(\sigma (t)1)x_0(t))𝑑t}`$ (11)
The concentration of the master-sequence is positive. The left hand side of Eq. 11 is therefore positive and the last term in the integral, $`(\sigma (t)1)x_0(t)`$, on the left hand side is negative. This implies that for $`x_0(t)`$ to be positive as time goes to infinity, we must assume $`\left(Q\sigma (t)1\right)𝑑t>0`$. The fluctuating time dependent equivalent to the static error threshold is therefore given by

$`Q_{fc}`$ $`=`$ $`{\displaystyle \frac{1}{\sigma _t}}`$ (12)
This shows that the error threshold on a fluctuating fitness peak is determined by the time average of the fitness, if the fluctuations are fast compared to the response time of the population.

Eq. 9 indicates that the response time of the system is approximately given by $`(Q\sigma (t)1)^1`$, i.e. the relative growth of the mastersequence compared to the background. For the time average mentioned above to be an interesting parameter the fluctuations of the fitness peak must therefore be faster than this response time; only for this kind of environmental dynamics is it sensible to talk in terms of the average concentration of the master-sequence. Thus if the fluctuations occur on a time-scale faster than the response-time of the quasispecies, then the error-threshold is defined by Eq. 12. For extremely slow changes the system will effectively be in equilibrium around the current value of the fitness. For slightly faster changes the response of the population will lag somewhat behind the changes in selective environment. In these cases it is more interesting to study the minimal concentration of the master sequence, which occurs when the fitness peak has a minimum (as we shall see later the phase-shift decreases when the fluctuation frequency decreases).

When the full replicator equations for a rapidly fluctuating peak are numerically integrated, the time-averaged quasispecies distribution displays an error catastrophe at high error rates $`\mu =1q`$. In figure 3 the fitness peak fluctuates periodically with $`\sigma (t)=10+5\mathrm{sin}(t)`$. The average fitness is given by $`\sigma _t=10`$ and the genome length $`\nu =25`$ and thus Eq. 12 predicts the error-threshold to occur at $`\mu =0.088`$, which agrees with the value found by numerically integrating the equations of motion directly. The analysis in this section demonstrates that by making the error tail approximation and reducing the dynamics to one-dimensional form, an analytic form exists for the error-threshold on fast moving landscapes. This one-dimensional formulation removes the need to time-order the changes in selective advantage of the landscape. This allows the integrals for the time history of the master-sequence concentration to be solved explicitly in equation 8.

## IV Phase-shifts on periodic landscapes

To study how the master sequence responds to changes in the height of the fitness peak it is convenient to assume that the fluctuations are periodic. It then makes sense to speak of the amplitude of the oscillations in concentration and of the phase-shift between the concentration and the fitness. It is intuitively clear that when the fitness peak is oscillating slowly (compared to the response time $`(Q\sigma (t)1)^1`$) there will be a very small phase-shift; the population will have time to reach an equilibrium about every value of $`\sigma (t)`$. The amplitude of changes in the master-sequence concentration will, for the same reason, be as large as possible. This result, together with the time-averaging effect found in the previous section, indicates that the population responds to the driving of the environment with a low pass filter effect. In one-dimensional population genetic models this phenomenon has been noted for some time . Wilke et al.  demonstrated via simulations that the same filtering occurred to quasispecies evolution on a periodically fluctuating single peak. Noting that the maxima and minima in concentration occurs when $`\dot{x}_0=0`$, we can find a relation between the phase-shift (between the concentration and fitness fluctuations), and the amplitude of the fitness fluctuations. Let $`t_{xmax}`$ be the time when the concentration has a maximum. Similarly the fitness is at a maximum at time $`t_{\sigma max}`$. Thus the phase-shift between the two is $`\delta =t_{xmax}t_{\sigma max}`$. From equation 6 the condition for the maximum value of $`x_0`$ during a full cycle can be derived

$`\underset{t}{\mathrm{max}}(x_0(t))`$ $`=`$ $`{\displaystyle \frac{Q\sigma (t_{\sigma max}+\delta )1}{\sigma (t_{\sigma max}+t_\delta )1}}`$ (13)
In general there is no closed analytic expression for this phase-shift ($`\delta `$), or the response amplitude of the master-sequence concentration. When the fluctuations of the fitness peak is a small harmonic oscillation equation 13 becomes analytically tractable. For such fluctuations

$`\dot{y}(t)`$ $`=`$ $`1(Q\sigma (t)1)y(t)`$ (14)
$`\sigma (t)`$ $`=`$ $`\overline{\sigma }+ϵ\mathrm{sin}(\omega t)`$ (15)
From equation 9 it is reasonable to assume the solution to be of the form $`y(t)=\frac{1}{Q\overline{\sigma }1}+u(t)`$, where $`u(t)`$ is small compared to the average. Ignoring higher order terms equation 6 can be written in terms of the perturbation $`u(t)`$ as

$`\dot{u}(t)`$ $`=`$ $`(1Q\overline{\sigma })u(t){\displaystyle \frac{ϵQ\mathrm{sin}(\omega t)}{Q\overline{\sigma }1}}`$ (16)
This differential equation can be solved to obtain

$`u(t)`$ $`=`$ $`{\displaystyle \frac{ϵQ}{(Q\overline{\sigma }1)\sqrt{(Q\overline{\sigma }1)^2+\omega ^2}}}\mathrm{sin}(\omega t\delta )`$ (17)
$`\mathrm{tan}(\delta )`$ $`=`$ $`{\displaystyle \frac{\omega }{Q\overline{\sigma }1}}`$ (18)
In eq. 17 and 18 transients have been ignored since they decay exponentially as $`e^{(Q\overline{\sigma }1)t}`$. Thus the frequency of the oscillations is normalized by the (average) response rate of the population $`Q\overline{\sigma }1`$.

Substituting this back into the expression for $`x_0(t)`$ gives

$`x_0(t)`$ $`=`$ $`{\displaystyle \frac{\overline{x}}{1\frac{ϵ(1Q)\mathrm{sin}(\omega t\delta )}{(\overline{\sigma }1)\sqrt{(Q\overline{\sigma }1)^2+\omega ^2}}}}`$ (19)
where $`\overline{x}=\frac{Q\overline{\sigma }1}{\overline{\sigma }1}`$.

The characteristic behavior of a low pass filter is clearly shown in equation 18 and 19. As the frequency of the fluctuations increases, the amplitude of the concentration response decreases and the phase shift converges to $`\frac{\pi }{2}`$. Figure 4 shows how a population responds to harmonic oscillations of the fitness peak. The phase-shift makes the concentration of the mastersequence reach its maximum when the actual fitness has already decreased below maximum.

## V Conclusions

In this paper we have shown that the time dynamics of a quasispecies on a fluctuating peak can be studied under the standard no back-mutation approximation. The general time ordering problem stemming from a time dependent landscape disappears since the rate equation becomes one–dimensional. We show that the time dependent equivalent to the static error threshold is determined by the time average of the fluctuations of the fitness peak. An expression for the typical response time for a population is given in terms of copying fidelity and selection pressure. We also show that for small periodic fluctuations the time dynamics of the population has a phase shift and a low pass filter amplitude response. Analytic expressions for the phase shift and the amplitude are derived in the special case of small harmonically oscillating fluctuations.

When doing this work Nigel Snoad and Martin Nilsson were supported by SFI core funding grants. N.S. would also like to acknowledge the support of Mats Nordahl at Chalmers University of Technology while preparing this manuscript. We also Mats Nordahl for valuable comments and discussions.