text:
# The network topology of a potential energy landscape: A static scale-free network

## Abstract

Here we analyze the topology of the network formed by the minima and transition states on the potential energy landscape of small clusters. We find that this network has both a small-world and scale-free character. In contrast to other scale-free networks, where the topology results from the dynamics of the network growth, the potential energy landscape is a static entity. Therefore, a fundamentally different organizing principle underlies this behaviour: The potential energy landscape is highly heterogeneous with the low-energy minima having large basins of attraction and acting as the highly-connected hubs in the network.

Energy landscapes have been at the forefront of many of the recent theoretical developments in our understanding of biomolecules,Brooks01  clustersWalesMW98 ; WalesDMMW00  and the glass transition.Debenedetti01  For example, this research has provided important new insights into how proteins foldBryngelson95  and the origin of the unusual properties of supercooled liquids, such as the distinction between “strong” and “fragile” liquids.Sastry01 ; SaikaVoivod01  There has also been a surge of interest in modelling complex systems as networks,Strogatz01  inspired by Watts and Strogatz’s discovery that many networks behave as “small worlds”.Watts98  Intriguingly, a diverse range of such networks, e.g. the world-wide web,Albert99  the internet,Faloutsos99  scientific collaborationNewman01a  and citationRedner98  networks, and biochemical networks,Jeong00 ; Jeong01  also have a “scale-free” topology, where the distribution of the number of connections to each node, the degree, follows a power law. This topology results from the dynamics of the network growth in these systems.Barabasi99  Here we draw these two strands of research together by applying the techniques of network analysis to probe the global structure of potential energy landscapes of clusters.

The potential energy landscape is a multi-dimensional surface representing the dependence of the potential energy on the positions of all the atoms of the system. For a system with many atoms the landscape will have a complex topography with higher-dimensional analogues of mountain ranges, valleys and passes. Although the potential energy landscape determines the system’s structure, thermodynamics and dynamics, the nature of this relationship is complex. A particularly successful means of elucidating this relationship is the inherent structure approach of Stillinger and Weber,StillW84a  in which the landscape is divided into basins of attraction surrounding each minimum (See Fig. 1). This partition provides a natural way to describe the dynamics of the system, because, except at high temperature, the system spends most of the time vibrating in the well surrounding a minimum and only occasionally hops to a different well by passing over a transition state. The interbasin dynamics can then be represented as a walk on a network whose nodes correspond to the minima and where edges link those minimum which are directly connected by a transition state. Figure 1 provides an illustration of such an inherent structure network (ISN) for a two-dimensional energy surface.

The ISN should provide the starting point for an energy landscape view of the global dynamics of a system. Indeed, as is increasingly being done, it is relatively easy to calculate the dynamics from this network using a master equation approach.BerryK95 ; Cieplak98a ; WalesDMMW00  However, fundamental questions about the topology of the ISN have received little attention. By contrast the global topography of energy landscapes has been the focus of much research.Brooks01 ; WalesMW98 ; Bryngelson95  To give one example, this emphasis has revealed that when a landscape is like a ‘funnel’Bryngelson95  the system is guided towards the global minimum, be it the native state of a protein,Bryngelson95  an ordered nanoparticleBall96  or a bulk crystal. However, the topological aspects of this idea remain open despite their importance: if the average number of steps to reach the global minimum from an arbitrary starting minimum scales unfavourably with size, the location of this structure would become significantly hindered at large size.

To characterize the topology of the ISN we study small clusters for which we are able to locate nearly all the minima and transition states on the potential energy landscape.WalesDMMW00 ; Tsai93a  The atoms of the cluster interact with a Lennard-Jones potential, which provides a reasonable description for rare gas clusters. The numbers of minima and transition states are expected to increase roughly as $`N_{\mathrm{min}}e^{\alpha N}`$ and $`N_{\mathrm{ts}}Ne^{\alpha N}`$, respectively,Still99 ; Doye02a  where $`N`$ is the number of atoms in the cluster. Therefore, the largest network that we are able to consider is for a 14-atom cluster for which we have located 4196 minima and $`\mathrm{87\hspace{0.17em}219}`$ transition states

Small-world networks have characteristics typical of both random graphs and lattices. The average separation between nodes scales logarithmically with network size, while the network is highly clustered, i.e. any two neighbours of a node are also likely to be connected. From Fig. 2 it is clear that the ISNs for the clusters show both these features and so are small worlds. The clustering is unsurprising given that the connections between basins on a potential energy landscape are based on adjacency in configuration space,cluster  but to interpret Fig. 2(a) properly we must take into account the increase in both the dimension of configuration space and the average degree, $`k`$, as the size of the network increases.

For example, for a hypercubic lattice with a constant number of lattice points, $`L`$, along each edge and dimension $`3N`$, the number of lattice points, $`N_{\mathrm{latt}}`$, scales exponentially with $`N`$ and the average number of steps between any two lattice points is $`3N(L+1)/3=(L+1)\mathrm{log}N_{\mathrm{latt}}/3\mathrm{log}L`$. By contrast, if the network behaves as a random graph, the average separation should scale as $`\mathrm{log}N_{\mathrm{min}}/\mathrm{log}k\mathrm{log}N_{\mathrm{latt}}/\mathrm{log}(\mathrm{log}N_{\mathrm{latt}})`$ because $`kN_{\mathrm{ts}}/N_{\mathrm{min}}N`$. The sublogarithmic scaling suggested by Fig. 2(a) points to the latter scenario. This result is somewhat surprising. In Watts and Strogatz’s small-world model the random-graph character results from the introduction of random links into the lattice, which can potentially connect up distant nodes, but there is no obvious equivalent of the random links on the potential energy landscapes.

If we now examine the distributions for the numbers of connections for each node we find that as the size of the cluster increases a clear power-law tail develops, which has a universal form independent of the cluster size (Fig. 3). The exponent of the power law, $`\gamma `$=2.78, is similar to other scale-free networks.Albert01  The cause of the random-graph like scaling of the average separation is thus the scale-free topology of the ISN. The network is extremely heterogeneous with a few hubs that have a very large number of connections, but with the majority of nodes only connected to a relatively small number of other minima.

This is a particularly surprising result because all other scale-free networks are dynamic in origin. They grow and change over time, be it on an almost continuous basis as in the WWW or on evolutionary time scales in the case of biochemical networks. Even the recently introduced deterministic scale-free networks are based on an iterative growth procedure.Barabasi01 ; Dorogovtsev02  Furthermore, models of network growthBarabasi99 ; Albert01  and studies on the time evolution of real networksNewman01d ; Jeong01b  strongly suggest that the heterogeneity at the heart of the scale-free topology develops as a result of new nodes preferentially linking to those nodes which have many connections, be they much-cited papers or popular web-sites. However, the network associated with a potential energy landscape is static. It is simply determined by the form of the interatomic interactions and the number of atoms in the cluster.

The source of the heterogeneity in the ISNs is apparent from Fig. 4, where we see that the degree of a node increases somewhat faster than exponentially as the energy of the minimum decreases.mindist  The low-energy minima act as the hubs in the network. Thus, for the 14-atom cluster, 76% of the nodes in the network are connected to the global minimum. To measure the extent of the catchment basin around a minimum, we can calculate the distance in configuration space to all the transition states connected to a minimum. For LJ14 we find that this distance is 2.7 times larger for the global minimum than for the surrounding minima. When the multi-dimensionality of configuration space is taken into account this result suggests that the hyperarea of this catchment basin is many orders of magnitude larger than the average. Similarly, it has been previously found that on average the basin area falls off approximately exponentially with the energy of a minimum.Doye98e

These results show that the global topology and topography of the potential energy landscape are intimately connected because the deep minima have very large catchment basins which are connected to lots of smaller basins around their edges. By contrast, if a potential energy landscape were flat and all basins of attraction had the same area the scale-free topology would be lost. For example, an investigation of the network topology of the configuration space of a non-interacting lattice polymerScala01b  found the connectivity distribution to be a Gaussian.Amaral00

This contrasting example naturally leads one to ask how general is the topology that we have found for the Lennard-Jones clusters. Although such a question can only be definitively answered by similar analyses for a variety of systems, there is nothing “special” about the Lennard-Jones potential and so there is no reason why similar behaviour should not be seen for other materials, as long as there are no constraints present that would prevent the formation of the high degrees associated with the hubs. A polymer provides an example of the latter, because the connectivity of the chain limits the number of transition states that can surround a minimum. For example, a similar analysis for Lennard-Jones polymersCalvo02a  did not find a power-law tail to the degree distribution (although it was still longer than exponential).JDunpub  There is no equivalent of many of the transition states for the equivalent Lennard-Jones cluster because they involve the breaking of the polymer chain.

The scale-free topology of the ISN is potentially good news for global optimization, the task of locating the global minimum. Even though the number of minima increases exponentially with the size of the system,Still99  the average number of steps in the shortest path to the global minimum grows sublinearly with system size. Of course, finding this path is not necessarily easy. Our calculations of the shortest paths required information on the global structure of the potential energy landscape, whereas a global optimization algorithm usually takes a step based on only local information.

Some path finding strategies to efficiently navigate scale-free networks have been suggested that make use of the fact that most of the shortest paths pass through the highly-connected hubs.Adamic01 ; Kim02  In our case the clear link between the topology and topography of the potential energy landscape provides an additional advantageous strategy. A downhill step to a lower-energy minimum is likely to take one to a minimum that is more connected and closer to the global minimum. How well-obeyed the latter correlation is, depends upon the global topography of the potential energy landscape and is a good indicator of the difficulty of global optimization. Thus, when the landscape is like a single funnel, global optimization algorithms can achieve near to the ideal scaling. For example, the basin-hopping algorithm can locate the global minimum of the 55-atom Lennard-Jones cluster after on average less than 150 minimizations when started from a random configuration, even though there are an estimated $`10^{21}`$ minima.Doye98e  The increasing number of links as the energy decreases evident in Fig. 4 further adds to the efficacy of a funnel in guiding the system towards its bottom, and provides clear evidence of the convergence of pathways into the hub at the funnel bottom that is often postulated. By contrast, when an energy landscape has multiple funnels and there is a tendency to enter a funnel that leads the system to a low-energy minimum that is far from the global minimum, global optimization can be very difficult.

The topology of the ISN will of course significantly affect the dynamics. This connection can be probed for very small systems where the network can be completely characterized and the inherent structure dynamics obtained by a master equation approach. However, this approach is not practical for the system sizes that are of most interest. Therefore, models of protein folding and the glass transition usually have to assume a simplified topology for the interstate dynamics,Shak89 ; Kohen00  or relate the dynamics to static quantities through phenomenological equations, such as the Adam-Gibbs equation which relates the relaxation time in supercooled liquids to the configurational entropy.AdamG65  To fully unlock the potential insights from the inherent structure view of the dynamics, a means of statistically modelling the network topology from a partial characterization of the potential energy landscape is therefore needed. Our results could significantly advance this goal.

The author is grateful to Emmanuel College, Cambridge and the Royal Society for financial support.


format: text

text:
# Coulomb Energy of Nuclei

## I Introduction

Relatively weak Coulomb interaction substantially affects the properties of nuclei owing to its long-range character. It is reliably established that the electrostatic Coulomb energy of the distributed charge is the dominant term in the Coulomb energy of nuclei. This term, also known as the Hartree energy, is proportional to $`Z^2e^2/R`$, where $`Z`$ is the number of intranuclear protons, $`R`$ is the nuclear radius, and $`e`$ is the proton charge; at large values of $`Z`$, it induces the breakup of nuclei. At the same time, there are a few less significant contributions to the Coulomb energy. These include the Fock exchange contribution and a number of terms associated with the correlated motion of nucleons. It will be shown below that, in constructing the density functional determining the Coulomb energy of a nucleus, it is sufficient to calculate all these terms to the first order in $`e^2`$. The presence of these terms is clearly illustrated by the well-known Nolen–Schiffer anomaly , which was discovered more than 30 years ago. This anomaly was deduced from a comparison of the calculated mass difference between two mirror nuclei with its experimental value. The comparison revealed that, on average, the results of the calculations performed without taking consistently into account the correlation Coulomb energy (CCE) fall short of relevant experimental values by 10% . This discrepancy can be removed by introducing charge-dependent forces  whose strength is assumed to be preset by this discrepancy. However, this is possible if we are sure that nuclear dynamics is taken properly into account and that all the Coulomb contributions are consistently included. A new mechanism that enhances the contribution of Coulomb interaction to the energy of the nuclear ground state and which is caused by the presence of a surface (a general property of equilibrium finite Fermi systems) was found in . This contribution, which is proportional to the nuclear surface $`(Z^{2/3})`$, made it possible to explain a dominant part of the above anomaly.

The main objective of this study is to calculate consistently the density functional determining the Coulomb energy. As the result of this calculation, which will be performed to the first order in $`e^2`$ and which will take into account the mechanism enhancing the CCE contribution, we will be able to clarify the effect of this mechanism on the single-particle spectrum of nuclei, on the effective nucleon mass, and on the position of the proton drip line. In what follows, we do not highlight the difference between this functional and the Coulomb energy unless this leads to confusion.

The ensuing exposition is organized as follows. A general formulation of our approach to calculating the Coulomb energy is given in Section 2. Section 3 is devoted to calculating the systematic CCE-induced shifts of the single-particle spectrum and the proton drip line. The effect of Coulomb interaction on the effective proton mass is considered in Section 4. Section 5 summarizes the basic results.

## II  Coulomb correlation energy

In order to formulate the aforementioned mechanism and to consider its effect, we will make use of the Hartree–Fock method employing effective forces, which was substantiated within density-functional theory . In this case, the energy $`E`$ of the nuclear ground state has the form

$$E=F_0[\rho _p,\rho _n]+F_c[\rho _p,\rho _n].$$
(1)
where $`F_0`$ is that part of the functional which is due to nuclear forces. Since we assume that the charge symmetry of nuclear forces takes place, $`F_0`$ depends identically on the single-particle densities $`\rho _p`$ and $`\rho _n`$ of protons and neutrons (these densities are determined by minimizing $`E`$). The term $`F_c`$ is due to other interactions, which are weaker than nuclear forces. The smallness of these interactions can be characterized by the ratio a of the potential energy of these interactions between two nucleons to the Fermi energy $`\epsilon _F40`$ MeV. For the Coulomb interaction, this ratio is $`\alpha =e^2/r_0\epsilon _F0.03`$, where $`r_0`$ is the mean spacing between nucleons in equilibrium nuclear matter. Therefore, it is sufficient to calculate the functional $`F_c`$ in the first order in the relevant coupling constant. In the case of Coulomb interaction, which is the main subject of the present study, the functional $`F_c`$ will be calculated to the first order in $`e_2`$. This functional could be approximated by some simple expression, as was done, for example, for the functional $`F_0`$. However, experience gained in calculating the mass differences between mirror nuclei (these mass differences are directly determined by the functional $`F_c`$) proves the inefficiency of such attempts. Moreover, this oversimplified approximation does not reveal the physics behind the functional $`F_c`$. Therefore, it is preferable to calculate the functional $`F_c`$. In doing this, it is assumed that the functional $`F_0`$ is known. It can be taken in the form specified by the Skyrme interaction  or in the form proposed in . The latter is advantageous in that it is characterized by a high precision, is simple in use, and involves the bare nucleon mass. By construction, the functional $`F_0`$ must contain the unnrenormalized mass $`M`$ of the proton (neutron) rather than its effective mass $`M^{}`$ because the effective mass, as well as the nuclear spectrum, cannot be directly included in it . The single-particle distributions of the proton and neutron densities ($`\rho _p`$ and $`\rho _n`$, respectively) are specified as

$$\rho _p(𝐫)=\underset{l}{}n_p^l|\varphi _p^l(𝐫)|^2;\rho _n(𝐫)=\underset{l}{}n_n^l|\varphi _n^l(𝐫)|^2,$$
(2)
where $`n_p^l`$ and $`n_n^l`$ are the occupation numbers for single-particle proton and neutron levels, respectively, while $`\varphi _p^l`$ and $`\varphi _n^l`$ are, respectively, the proton and neutron single-particle wave functions determined from the Hartree–Fock equation . The term $`F_c[\rho _p]`$ is usually taken in the first order in the Coulomb interaction in the Hartree–Fock approximation :

$$F_c=\frac{e^2}{2}\left[\frac{\rho _p(𝐫_1)\rho _p(𝐫_2)}{|𝐫_1𝐫_2|}𝑑𝐫_1𝑑𝐫_2\left[\frac{\chi _p^0(𝐫_1,𝐫_2,i\omega )+2\pi \rho _p(𝐫_1)\delta (𝐫_1𝐫_2)\delta (\omega )}{|𝐫_1𝐫_2|}\right]\frac{d𝐫_1d𝐫_2d\omega }{2\pi }\right].$$
(3)
Here, $`\chi _p^0(𝐫_1,𝐫_2,\omega )`$ is the linear response function for noninteracting protons moving in the self-consistent single-particle field $`V_p`$. The first and the second term on the right-hand side of Eq. (3) are, respectively, the Hartree and the Fock density functional ($`F_c^H`$ and $`F_c^F`$). The exchange term is usually approximated by the expression ,

$$F_c^F=\frac{3e^2}{4}\left(\frac{3}{\pi }\right)^{1/3}\rho _p^{4/3}(𝐫)𝑑𝐫.$$
(4)
However, the equality in (3) is not exact even in the first order in the Coulomb interaction ($`e^2`$ order) because it does not include the functional $`F_c^{corr}[\rho _p(𝐫)]`$ taking into account the contribution to the Coulomb energy from the correlated motion of protons under the effect of the effective (residual) nuclear interaction $`R_{\alpha \beta }(𝐫_1,𝐫_2)`$. In the first order in $`e^2`$, the CCE functional $`F_c^{corr}`$ has the form

$$F_c^{corr}[\rho _p]=\frac{e^2}{2}\left[\frac{\chi _{pp}(𝐫_1,𝐫_2,i\omega )\chi _p^0(𝐫_1,𝐫_2,i\omega )}{|𝐫_1𝐫_2|}\right]\frac{d𝐫_1d𝐫_2d\omega }{2\pi },$$
(5)
where $`\chi _{pp}(𝐫_1,𝐫_2,\omega )`$ is the exact function representing the linear response of intranuclear protons to an external electric field. It should be noted that the function $`\chi _{pp}`$ is completely determined by the functional $`F_0`$  and that it involves no smallness associated with the Coulomb interaction. The total response function contains three independent components — $`\chi _{pp}`$, $`\chi _{nn}`$ and $`\chi _{np}`$ — and satisfies the equation

$$\chi _{lm}=\chi _l^0\delta _{lm}+\underset{k}{}\chi _l^0R_{lk}\chi _{km}.$$
(6)
where integration signs are omitted for the sake of simplicity. The response function $`\chi _p^0`$ for noninteracting nucleons (protons) moving in the self-consistent field $`V_p`$ has the form

$$\chi _p^0(𝐫_1,𝐫_2,\omega )=\underset{i,k}{}n_p^i(1n_p^k)\varphi _p^i(𝐫_1)\varphi _p^i(𝐫_2)\varphi _p^k(𝐫_2)\varphi _p^k(𝐫_1)\left[\frac{1}{\omega \omega _p^{ik}+i\eta }\frac{1}{\omega +\omega _p^{ik}i\eta }\right].$$
(7)
The single-particle eigenvalues $`\epsilon _p^i`$, the frequencies $`\omega _p^{ik}=\epsilon _p^k\epsilon _p^i`$, the single-particle wave functions $`\varphi _p^i(𝐫)`$, and the corresponding densities $`\rho _p`$ and $`\rho _n`$ are determined by solving the single-particle Hartree–Fock–like equations, or the Kohn–Sham equations ,

$$\left(\frac{^2}{2M}+V_p(𝐫)\right)\varphi _p^i(𝐫)=\epsilon _p^i\varphi _p^i(𝐫),$$
(8)
which are derived by means of a conventional procedure that amounts to varying Eq. (1) with respect to the single-particle functions $`\varphi _{n,p}^l`$ . The self-consistent potential $`V_p`$ and the effective interaction $`R_{\alpha \beta }`$ are given by

$$V_p(𝐫)=\frac{\delta (F_0[\rho _p,\rho _n]T_0[\rho _p,\rho _n])}{\delta \rho _p(𝐫)};R_{\alpha \beta }(𝐫_1,𝐫_2)=\frac{\delta ^2(F_0[\rho _p,\rho _n]T_0[\rho _p,\rho _n])}{\delta \rho _\alpha (𝐫_1)\delta \rho _\beta (𝐫_1)},$$
(9)
where $`T_0[\rho _p,\rho _n]`$ is the kinetic-energy functional for noninteracting nucleons. Thus, we conclude that, to the first order in $`e^2`$, the eventual form of the functional $`F_c`$ is

$$F_c=F_c^H+F_c^F+F_c^{corr}.$$
(10)
Equations (1), (3), (5), and (10) provide a basis for consistently calculating the Coulomb energy and for taking into account the Coulomb interaction effect on the properties of nuclei. The poles of the response function c $`\chi _{pp}`$ determine the collective spectrum of nuclei, while the residues at these poles govern the relevant transition probabilities. As can be seen from Eq. (5), a dominant contribution to $`F_c^{corr}`$ comes from collective isoscalar surface vibrations whose excitation energies are much lower than those of corresponding isovector modes. Thus, the functional $`F_c^{corr}`$ is controlled by the isoscalar components of the effective interaction.

It is instructive to consider the CCE for symmetric infinite nuclear matter , in which case the relevant equations are considerably simplified. The function c $`\chi _{pp}`$ assumes the form

$$\chi _{pp}(q,\omega )=\frac{1}{2}\left[\frac{\chi ^0(q,\omega )}{1R_+(q,\omega ,\rho )\chi ^0(q,\omega )}+\frac{\chi ^0(q,\omega )}{1R_{}(q,\omega ,\rho )\chi ^0(q,\omega )}\right],$$
(11)
and determines the functional $`F_c^{corr}`$:

$$F_c^{corr}[\rho _p]=2\pi e^2\left[\frac{\chi _{pp}(q,i\omega )\chi _p^0(qi\omega )}{q^2}\right]\frac{d𝐪d\omega }{(2\pi )^4}.$$
(12)
For symmetric matter, we have $`\chi _p^0=\chi _n^0=\chi ^0/2`$; the effective interaction $`R(\rho )`$ is similar to the local interaction $`f(\rho )`$ between nucleons that was introduced by Migdal :

$$\frac{p_FM}{\pi ^2}R_+(\rho )=f(\rho );\frac{p_FM}{\pi ^2}R_{}(\rho )=f^{}(\rho ).$$
(13)
Considering that the isoscalar amplitude $`f(\rho 0)`$ is approximately equal to $`2.5`$ and that $`p_FM/\pi ^2=\chi ^0(0,i0)`$, one can see from Eq. (11) that, at relatively low densities corresponding to the surface region of nuclei, the denominator $`1R_+\chi ^0`$ vanishes at small values of $`q`$ and $`\omega `$. Therefore, the function $`\chi _{pp}`$ develops poles on the imaginary axis. These poles, which evince instability of low-density nuclear matter , lead to the divergence of the integral in Eq. (12), and this prevents a calculation of the CCE at these nuclear-matter densities. We can conclude that the main contribution to this integral comes from the isoscalar response function. From the above, it also follows that the conventional procedure based on calculating various quantities in the local density approximation (LDA) in infinite uniform matter fails in this case. In this connection, it is worth noting that, although an attempt was made in  to construct the functional $`F_0`$ within the LDA, this is impossible for the same reason. It is convenient to calculate the CCE in semi-infinite nuclear matter, where there is a region of low-density nuclear matter near the surface. In the generalized LDA, we then obtain ,

$$F_c^{corr}[\rho _p]=\rho _p(𝐫)e_c(\rho (𝐫))𝑑𝐫,$$
(14)
where $`e_c(\rho )`$ is the CCE per proton and $`\rho _p`$ is the single-particle proton density. The calculation of the energy $`e_c(\rho _p)`$ for semi-infinite nuclear matter revealed that, at the surface, this energy has a pronounced positive peak, which corresponds to smoothing the divergence in (12) for uniform nuclear matter . It is convenient to approximate the energy $`e_c`$ by the simple expression

$$e_c(\rho _p(r))=D\left[\frac{\rho _p(r)}{\rho _0}\left(1\frac{\rho _p(r)}{\rho _0}\right)\right]^{4/3},$$
(15)
where $`2\rho _0=0.16fm^3`$ is the equilibrium nuclear density and $`D6`$ MeV. In contrast to the terms $`F_c^H`$ and $`F_c^F`$, the energy $`e_c(\rho _p)`$, which has a pronounced peak at the nuclear surface, makes a noticeable contribution to the surface tension $`\sigma _c`$,

$$\sigma _c=_{\mathrm{}}^{\mathrm{}}\rho _p(z)[e_c(\rho _p(z))e_c(\rho _p(\mathrm{}))]𝑑z.$$
(16)
Hence, the Weizsöcker mass formula must be supplemented with the term $`\mathrm{\Delta }E`$ whose contribution to the total binding energy of a nucleus can be represented as

$$\mathrm{\Delta }E=a_1^{}Z+a_2^{}Z^{2/3},$$
(17)
where $`a_1^{}0.1`$ MeV and $`a_2^{}1.0`$ MeV. Thus, we can see from Eq. (17) that the surface tension has an isovector component. It will be shown below that this is so for the effective mass as well. It follows from Eq. (17) that the surface tension of the proton nuclear Fermi liquid is effectively greater than that of neutron matter. This means that, when protons are added to a nucleus, their binding energy becomes somewhat smaller because of an increase in the surface of the proton liquid component. It will also be proven below that this decrease in the proton binding energy is sufficient for explaining the mass difference between mirror nuclei. Further, the additional surface tension $`\sigma _c`$ and the addition of extra neutrons are expected to induce opposite modifications in the root-mean-square radius $`<r_p>`$ of the intranuclear-proton distribution: the former would moderate an increase associated with the latter. The preliminary calculations reveal that this could explain the anomalously small increase in the radius $`<r_p>`$ in going over from the 40Ca to the 48Ca nucleus. The contribution $`\sigma _c`$ must also be taken into account in considering fission barriers for heavy nuclei.

## III Single-particle spectrum and the proton drip line

Let us calculate the CCE-induced shift of the single-particle proton excitation spectrum $`\epsilon _p^l`$. For this, we will use the well-known Landau equation

$$\frac{\delta E}{\delta n_p^l}=\epsilon _p^l.$$
(18)
It follows from Eqs. (5) and (18) that this shift $`\mathrm{\Delta }\epsilon _p^l`$ can be represented as

$$\mathrm{\Delta }\epsilon _p^l=\frac{\delta F_c^{corr}}{\delta n_p^l}=\frac{e^2}{2}\frac{\delta }{\delta n_p^l}\left[\frac{\chi _{pp}(𝐫_1,𝐫_2,i\omega )\chi _p^0(𝐫_1,𝐫_2,i\omega )}{|𝐫_1𝐫_2|}\right]\frac{d𝐫_1d𝐫_2d\omega }{2\pi }.$$
(19)
The variational derivative $`\delta \chi _p^0/\delta n_p^l`$ has the simple functional form

$$\frac{\delta \chi _p^0(𝐫_1,𝐫_2,\omega )}{\delta n_p^l}=\left[G^p(𝐫_1,𝐫_2,\omega +\epsilon _p^l)+G^p(𝐫_1,𝐫_2,\omega +\epsilon _p^l)\right]\varphi _p^l(𝐫_1)\varphi _p^l(𝐫_2),$$
(20)
where $`G^p(𝐫_1,𝐫_2,\omega )`$ is the single-particle Green’s function for the system of $`Z`$ noninteracting protons moving in a self-consistent single-particle nuclear potential. The functional derivative $`\delta \chi _{pp}/\delta n_p^l`$ is determined by the matrix equation

$$\frac{\delta \chi _{lm}}{\delta n_p^l}=\frac{\delta \chi _l^0}{\delta n_p^l}\delta _{lm}+\underset{k}{}\left[\frac{\delta \chi _l^0}{\delta n_p^l}R_{lk}\chi _{km}+\chi _l^0\frac{\delta R_{lk}}{\delta n_p^l}\chi _{km}+\chi _l^0R_{lk}\frac{\delta \chi _{km}}{\delta n_p^l}\right].$$
(21)
which can be derived by directly varying Eq. (6). Integration with respect to spatial coordinates is implied in Eq. (21) in just the same way as in Eq. (6). As can be seen from Eq. (9), the effective interaction $`R_{\alpha \beta }`$ is determined by the form of the functional $`F_0`$. In order to simplify the calculations, we took, however, the interaction $`R_{\alpha \beta }`$ in the separable representation ,

$$R_{lk}(𝐫_1,𝐫_2)=\lambda \frac{dV_l(r_1)}{dr}\frac{dV_k(r_2)}{dr}\delta (\mathrm{\Omega }_1\mathrm{\Omega }_2),$$
(22)
where $`V_l(r)`$ is the self-consistent single-particle proton (neutron) potential. The value of the parameter $`\lambda `$ is chosen in such a way that the dipole response has a pole at $`\omega =0`$. The above separable form of the effective (residual) interaction is extensively used and, as was shown in , provides a good description of collective nuclear excitations. The value calculated here for the shift $`\mathrm{\Delta }\epsilon _p^l`$ of single-particle proton levels located near the Fermi level is $`(0.20.4)`$ MeV both for medium-mass and heavy nuclei, in overall agreement with the value of the Nolen–Schiffer anomaly. It is worthwhile to verify these results by using the simple LDA expressions (14) and (15) for $`F^{corr}[\rho _p]`$. The shift $`\mathrm{\Delta }\epsilon _p^l`$ then assumes the form

$$\mathrm{\Delta }\epsilon _p^l=\frac{\delta F_c^{corr}[\rho _p]}{\delta \rho _p}|\varphi _p^l(𝐫)|^2𝑑𝐫.$$
(23)
It is convenient to approximate the single-particle density $`\rho _p(r)`$ by the Fermi distribution

$$\rho _p(r,R)=\frac{\rho _0}{1+\mathrm{exp}((rR)/a)}.$$
(24)
where $`R`$ is the nuclear radius and the diffuseness parameter $`a`$ is taken to be 0.6 fm. In this case, the functional $`F_c^{corr}[\rho _p]`$ given by Eq. (14) reduces to a function $`F_c^{corr}(R)`$ of the nuclear radius. An addition of a proton increases the radius $`R`$ by the quantity $`\mathrm{\Delta }R`$ that is determined from the normalization condition

$$\left[\rho _p(r,R+\mathrm{\Delta }R)\rho _p(r,R)\right]𝑑𝐫=1.$$

The CCE-induced shift of a proton level occurring near the Fermi surface is given by

$$\mathrm{\Delta }\epsilon _p=F_c^{corr}(R+\mathrm{\Delta }R)F_c^{corr}(R).$$

After simple transformations of the relevant integrals of the Fermi functions, we find that, for medium-mass and heavy nuclei, the sought shift of a proton level near the Fermi level is

$$\mathrm{\Delta }\epsilon _p\frac{Da}{2R}(0.30.4)\mathrm{M}eV.$$
(25)
Since corrections of order $`D(a/R)^3`$ were discarded, the result in (25) slightly overestimates the shift.

The calculated mass differences between mirror nuclei are quoted in the table. In this calculation, the functional $`F_0`$ was constructed on the basis of the $`SIII`$ interaction , while the functionals $`F_c`$ and $`F_c^{corr}`$ were, respectively, taken in the form (10) and specified by Eqs. (14) and (15).

Mass differences between mirror nuclei (in MeV)
Nuclei 15O – 15N 17F – 17O 39Ca – 39K 41Sc – 41Ca 48Ni – 48Ca    Theory 3.48 3.56 7.23 7.24 66.70  Experiment  3.54 3.54 7.30 7.28 67.06

The table demonstrates that the mass differences between mirror nuclei are closely reproduced. The remaining disagreement can be used to determine the coupling constants for forces violating the charge symmetry of nuclear interaction. Thus, traditional nuclear physics can be advantageous in calculating the constants concerning elementary-particle physics. Equation (25) shows that the shift $`\mathrm{\Delta }\epsilon _p`$ of a single-particle level approximately compensates for the level shift caused by exchange Coulomb interaction. It should be recalled that the exchange Coulomb interaction reduces the energy of a single-particle level, as can be seen from Eq. (4). Therefore, the total shift of a single-particle level due to Coulomb interaction can be derived by taking into account solely the direct Coulomb interaction—that is, by retaining only the term $`F_c^H`$ in Eq. (10). A similar procedure was postulated in . It should be noted that the recipe assuming mutual cancellation of the exchange Coulomb interaction and the CCE does not lead to a noticeable value of the coefficient $`a_2^{}`$ in Eq. (17) (recall this coefficient determines the isovector component of the surface tension) because the contribution from $`F_c^H`$ and $`F_c^F`$ \[see Eq. (4)\] to the surface tension is small. On the contrary, our calculations indicate that the coefficient $`a_2^{}`$ is not small. It also follows from Eq. (25) that, upon taking the CCE into account, the last filled single-particle proton level is shifted upward by 0.3 to 0.4 MeV. As a result, the maximal charge $`Z`$ of a nucleus occurring near the proton drip line decreases by 2 or 3 units, as follows from . Therefore, we conclude that, in calculating the proton drip line, it is important to take the CCE into account.

## IV Effective mass

Let us consider the change $`\mathrm{\Delta }M`$ that Coulomb interaction induces in the effective proton mass $`M^{}`$. In the case of uniform nuclear matter, the single-particle spectrum depends on the momentum $`p`$, while the effective mass, as follows from Eq. (18), is determined by the expression

$$\frac{1}{M^{}}=\frac{1}{p_F}\frac{d\epsilon _p(p)}{dp}|_{p=p_F},$$
(26)
where $`p_F`$ is the Fermi momentum. In order to calculate $`\mathrm{\Delta }M`$, we use Eq. (26), replacing $`\epsilon _p(p)`$ by the shift of the single-particle spectrum due to Coulomb interaction. We have

$$\frac{\mathrm{\Delta }M}{M^{}(M^{}+\mathrm{\Delta }M)}=\frac{e^2}{p_F}\frac{d}{dp}\frac{\delta }{\delta n_p}\left[\frac{\chi _{pp}(q,i\omega )+2\pi \rho _p\delta (\omega )}{q^2}\right]\frac{dqd\omega }{(2\pi )^3}.$$
(27)
Evaluating variations in Eq. (27) and considering that, in the case of uniform matter, Eqs. (6) and (21) reduce to algebraic equations, we obtain ,

$$\frac{\mathrm{\Delta }M}{M^{}(M^{}+\mathrm{\Delta }M)}=\frac{e^2d}{p_Fdp}\frac{\delta \chi _0(q,i\omega )}{\delta n_p}\frac{1}{(1R_+(q,i\omega )\chi _0(q,i\omega ))^2}\frac{d𝐪d\omega }{q^2(2\pi )^3}.$$
(28)
where $`M^{}`$ is the effective proton mass in the absence of Coulomb interaction and $`R_+=(R_{pp}+R_{pn})/2`$ is the effective interaction. The derivative $`d/dp`$ was taken at the point $`p=p_F`$. Let us consider the effective-mass variation $`\mathrm{\Delta }M`$ for the case where the system under consideration is characterized by the compressibility parameter $`\kappa `$, $`1/\kappa =\rho ^2d^2E/d\rho ^2`$, tending to infinity. This state is similar to states existing in the surface region of nuclei. This consideration will enable us to draw qualitative conclusions on the behavior of the effective mass in a nucleus. It can be shown that

$$\frac{d}{dp}\frac{\delta }{\delta n_p}\chi _0(q,\omega )|_{pp_F}=\frac{4\pi }{p_F^2}\delta (p_F|𝐩+𝐪|)\delta (\omega )𝐩(𝐩+𝐪).$$
(29)
Substituting (29) into (28) and retaining only the leading term, we obtain

$$\frac{1}{M^{}+\mathrm{\Delta }M}=\frac{1}{M^{}}+\frac{e^2}{2\pi p_F}_1^1\frac{xdx}{(1x)[1R(q(x),0)\chi _0(q(x),0)]^2}.$$
(30)
where $`q(x)=p_F\sqrt{2(1x)}`$. At the point where the compressibility diverges, the denominator $`(1R_+\chi _0)^2`$ of the integrand in Eq. (30) vanishes at $`x=1(q=0)`$. The relevant integral, which is positive, diverges, which results in the vanishing of the effective mass—that is, $`M^{}+\mathrm{\Delta }M0`$. This result demonstrates that, in equilibrium nuclear systems, the CCE noticeably affects the effective mass, the decrease $`\mathrm{\Delta }M`$ in the effective proton mass being dependent on the degree to which the density of the system,$`\rho `$, deviates from its equilibrium value $`\rho _0`$. In a nucleus, the divergence of the integral in Eq. (30) — it could be due to the decrease in the density near the surface of the nucleus — gives way to a finite value. As a result, the effective proton mass decreases in relation to $`M^{}`$ and becomes smaller than the effective neutron mass. It can be said that the effective nucleon mass acquires an isovector component . This possibility was discussed in . A theoretical validation of this effect has been given here.

## V  Conclusion

The basic results of the present study can be summarized as follows. A consistent scheme for constructing the density functional determining the Coulomb energy of a nucleus has been developed to the first order in $`e^2`$. Using this functional, we have calculated the single-particle spectra of nuclei and the systematic shift of these spectra that is induced by the CCE. It has been shown that the Nolen–Schiffer anomaly is removed to a considerable extent by this systematic shift of the single-particle spectrum. Owing to the same mechanism, the proton drip line undergoes a shift of 2 or 3 units toward smaller values of the charge $`Z`$ of a nucleus occurring near this line. It has been shown that the CCE must be taken into account in calculating the effective nucleon mass. This opens the possibility of estimating the coupling constants for forces violating charge symmetry. The contribution of these forces must be treated in the same manner as this has been done for the Coulomb interaction. It can be expected that the proposed procedure will make it possible to construct a density functional applicable to describing various properties of nuclei occurring both in the valley of stability and beyond it.

## VI Ackowledgments

I am grateful to A. Bulgac for stimulating discussions on the problems analyzed in this study. This work was supported in part by INTAS-OPEN-97-603.


text:
# Non-Radiative Accretion and Thermodynamics

## 1 An unsolved problem: the rate of non-radiative quasi-spherical accretion

We want to calculate the rate of quasi-spherical non-radiative accretion onto a black hole of Schwarzschild radius $`R_S=2GM/c^2`$. We assume that $`R_SR_A`$, where the accretion radius $`R_A=2GM/c_s^2`$. Here $`c_s`$ is the speed of sound at infinity, where there is gas at rest, with mass density $`\rho `$.

A standard answer is the Bondi (1952) rate:

$$\dot{M}_{\mathrm{Bondi}}R_A^2c_s\rho .$$
(1)
This indeed is the correct rate for an ideal spherically symmetrical stationary flow of compressible fluid with adiabatic index less than 5/3. Under spherical symmetry, for an ideal gas, one can get Bondi in a direct numerical simulation.

But: (i) the spherical symmetry is almost certainly broken by turbulence, (ii) turbulence makes the flow non-stationary, (iii) magnetic fields are probably generated, (iv) dissipative processes appear. Then Bondi is not immediately applicable, but it can still be a good estimate. Indeed, the rate of mass accretion should be about $`\dot{M}R_S^2c\rho _S`$, where $`\rho _S`$ is the density near the hole, where everything goes in at about the speed of light. Now we want to estimate $`\rho _S`$ in terms of $`\rho `$. This can be done by requiring that the mass flux be constant, $`R^2v(R)\rho (R)=`$const. Assuming that $`v(R)`$ is close to the free fall, $`R^{1/2}`$, we get $`\rho (R)R^{3/2}`$, and therefore the Bondi rate.

Let us do it again: the rate of mass accretion should be about $`\dot{M}R_S^2c\rho _S`$. To estimate $`\rho _S`$, assume that thermal energy generated by accretion is transported outward by a Shakura-Sunyaev-like turbulence (the flow is non-radiative) at a constant rate $`R^2c_s^3(R)\rho (R)=`$const. This gives a 1/2-law density profile $`\rho (R)R^{1/2}`$, and a much smaller mass accretion rate

$$\dot{M}_{1/2\mathrm{law}}R_AR_Sc_s\rho .$$
(2)
This small accretion rate, or equivalently a shallow density profile, is the main feature of CDAF, which is a particular worked out example of the 1/2-law flow (Gruzinov 2001, Quataert & Gruzinov 1999, Balbus & Hawley 2002 and references therein). I take CDAF to mean constant energy flux and not all of math details in all the papers (including ours, Quataert & Gruzinov 1999) e.g., marginal stability, …. Very reasonable to admit that math details are different in the magnetic case, but one can still have a constant energy flux as a defining feature of the flow.

Which answer is correct? The three dimensional time-dependent MHD (magnetohydrodynamics) simulation is the right way to solve the problem (Balbus & Hawley 2002). Or is it?

To see what are the typical plasma parameters, consider accretion onto the black hole at the center of the Galaxy. A recent Chandra detection of quiescent Sgr A (Baganoff et al 2001) can be understood as thermal bremsstrahlung coming from about the vicinity of accretion radius $`R_A10^{17}cm`$, with gas number density at accretion radius $`n10^3\mathrm{cm}^3`$. (Quataert & Gruzinov 2000, Quataert 2002). Then, assuming Bondi, we can estimate the gas temperature and density at smaller radii. Say at $`R=100R_S=10^{14}\mathrm{cm}`$, we have a plasma of temperature $`T(R)1`$MeV, and $`n(R)(R/R_A)^{3/2}n3\times 10^7\mathrm{cm}^3`$. This gives a mean free path (from Coulomb collisions) $`\lambda 3\times 10^{16}\mathrm{cm}R`$. 111 $`\lambda 10^{12}\mathrm{cm}(T/1\mathrm{e}\mathrm{V})^2/(n/1\mathrm{c}\mathrm{m}^3)`$ The mean free path is even larger if the 1/2-law density profile is used. Thus we have a collisionless plasma. MHD is not applicable (because particles make large excursions along magnetic field lines).

A direct numerical simulation of a collisionless plasma seems impossible today – even the MHD simulations have inadequate resolution to see if the Bondi, $`nR^{3/2}`$, or the 1/2-law, $`nR^{1/2}`$, flow forms (Hawley et al 2001). In this situation, it would be of great importance, if thermodynamic arguments could be used to estimate the rate of accretion (or equivalently the density profile). Indeed, Balbus and Hawley (2002) have argued that thermodynamics rules out CDAFs.

This cannot be correct because CDAFs were obtained in a thermodynamically consistent simulation. (A simulation is thermodynamically consistent if entropy is growing, that is “all viscosities are positive”). Thus the 1/2 law is thermodynamically allowed. But this is not the point of the present note.

Here we will show that all accretion rates smaller than Bondi but larger than the 1/2-law are thermodynamically admissible. Therefore, until a collisionless plasma simulation is done, one won’t be able to say what is the true accretion rate in Sgr A and in many other Dim Black Holes.

## 2 Thermodynamically admissible collisionless non-radiative quasi-spherical accretion flows

Bondi: We start with the standard lore. In a collisionless plasma, electromagnetic field will appear. Thermodynamics does not contradict the assumption that a weak small-scale magnetic field is the only field generated. This field will make an effective mean free path small. Plasma becomes a fluid. Thermodynamics does not contradict the assumption that the resulting flow will be spherical, stationary, with negligible dissipation. Under these assumptions, we get Bondi, equation (1).

1/2 law: Thermodynamics does not contradict the assumption that no electromagnetic fields are generated. Then the plasma is just a collection of non-interacting particles in a gravitational potential. At large times $`t`$, particles at distances $`Rc_st`$ will have a chance of flying into the hole. To hit the hole, the particle’s velocity should lie within the loss cone of opening angle $`\theta \sqrt{R_AR_S}/R`$. The total mass accreted into the hole at times $`t`$ is $`\rho R^3\theta ^2`$, and the corresponding mass accretion rate is the 1/2-law rate given by equation (2).

Clearly, by adjusting our (thermodynamically admissible) assumptions about the nature of the flow, we can get any accretion rate between Bondi and 1/2-law.

## 3 Conclusions

* Thermodynamics alone cannot solve the problem of non-radiative quasi-spherical accretion flows. Very different flow regimes, Bondi, 1/2-law, and intermediate (Gruzinov 1999) are thermodynamically admissible. This is true for fluids as well as for collisionless plasmas.
* Three dimensional MHD numerical simulations is the most reliable approach today. If adequate resolution is achieved (to reach the self-similar flow regime which almost certainly exists), the resulting rate of mass accretion will have to be taken seriously.
* However, in many astronomically interesting cases the plasma is collisionless, and even 3D MHD is potentially misleading. Large non-MHD terms (say Shakura-Sunayev-like transport coefficients included into the basic model) might be present. The only way to get a true rate of mass accretion in Sgr A and in other Dim Black Holes is direct simulation of collisionless plasma, which today looks like sheer impossibility.

I thank Eliot Quataert for useful discussions; several sentences in this note are his. I thank my NYU colleagues Roman Scoccimarro and Matias Zaldarriaga for bringing the preprint of Balbus and Hawley to my attention and demanding explanations.


text:
# Strongly decelerated expansion of SN 1979C

## 1 Introduction

Supernova SN 1979C in M100 was discovered on 1979 April 19 by Gus E. Johnson (Mattei et al. mat79 (1979)) at around magnitude 12 in the visible, and it is thought to have exploded on April 4, 1979 (Weiler et al. wei86 (1986)). Optical spectra first showed a featureless continuum, which later evolved to exhibit strong H$`\alpha `$ emission (e.g., Schlegel sch96a (1996)). SN 1979C was extraordinarily luminous, $`M_B^{\mathrm{max}}20`$ (e.g., Young & Branch you89 (1989)), making it among the most luminous type II supernova ever observed. SN 1979C was also found to be a radio supernova (RSN)(Weiler and Sramek, 1980). Indeed, at a distance of $`16.1\pm 1.3`$ Mpc estimated for M100 (Ferrarese et al. fer96 (1996)) it is one of the intrinsically strongest RSN (Weiler et al. wei96 (1996)). From H$`\alpha `$ widths an expansion speed of 9 200$`\pm `$500 km s-1 was estimated (Panagia et al. pan80 (1980) and G. Vettolani, priv. comm.) for epoch around 45 d.

SN 1979C has been classified as a SN of type II-L. Fesen et al. (1999), hereafter F99, have studied the late-time optical properties of SN 1979C and other three type II-L supernovae. Late-time optical emission lines from type SNe II-L appear surprisingly steady over long time intervals. These authors find similar expansion velocities in the range 5100–6000 km s-1 in the strongest lines present (H$`\alpha `$, \[Oi\] 6300, 6364 Å, \[Oii\] 7319, 7330 Å, and \[Oiii\] 4959, 5007 Å) for all type II-L supernovae. From the asymmetries of the lines these authors find evidence of dust formation. Strong radio emission from a SN appears to be strongly correlated with late time bright optical emission, with SN 1979C as a prime example. In SN 1979C, all detected lines but H$`\alpha `$ show emission peaks at $`\mathrm{5\hspace{0.17em}000}`$ and $`\mathrm{1\hspace{0.17em}000}`$ km s-1 suggestive of clumpy emission from material coming mostly from ejecta, and physically separated from the material emitting H$`\alpha `$. The H$`\alpha `$ emission tracks well the peak 6 cm emission (Weiler et al. wei91 (1991)) indicating that both emissions are related to the material from the swept-up shell.

The evolution of the radio emission from SN 1979C has been interpreted within the mini-shell model (Chevalier che82 (1982)) as due to synchrotron emission from the outer part of the shock swept-up shell of circumstellar material attenuated by thermal absorption from the even more distant ionized CSM (Weiler et al. wei86 (1986), Montes et al. mon00 (2000)). Synchrotron self absorption which has been shown to be relevant in SN 1993J (Fransson & Björnsson fra98 (1998), Pérez-Torres et al. per01 (2001)) might also play a role in SN 1979C, although Chevalier (che98 (1998)) is not of this opinion. The modulation present in the radio light curves of SN 1979C has led Weiler et al. (wei92 (1992)) and Montes et al. (mon00 (2000)) to model the progenitor as possibly being in a detached eccentric binary system with a less massive companion. ¿From a study of the emission environment of SN 1979C, Van Dyk et al. (vdy99 (1999)) have estimated the mass of the progenitor to be $`1718(\pm 3)`$ M. Hydro-dynamical simulations by Schwarz & Pringle (sch96 (1996)) confirm that a binary system similar to that proposed by Weiler et al. (wei92 (1992)) with a spiral shaped stellar wind around the progenitor is feasible for SN 1979C. Such modulation of the density of the gas surrounding the RSN and now being shocked, might have dramatic effects on the shape of the radio structure. A determination of such radio structure would be a powerful way to further understand the pre-supernova phase of the progenitor of SN 1979C.

The radio emission of this RSN, monitored by Montes et al. (mon00 (2000)), has started to show an interesting trend which can be seen in Fig. 1: after years of steady decline the radio emission has been constant, or perhaps increasing, since about 1990 ($`10`$ yr) in a way somewhat reminiscent of the radio emission increase for the case of SN 1987A (Ball et al. ball95 (1995)). It appears that after a standard RSN radio emission decline phase (Weiler et al. wei96 (1996)) the expanding shock may have encountered denser ionized gas which has boosted the radio emission. These results may be in accord with those of F (99) which show that high electron densities (105-7 cm-3) are necessary to explain the absence of several lines, such as \[Oii\] 3726, 3729 Å in their optical spectra.

Soft X-rays have been detected from SN 1979C by Immler et al. (imm98 (1998)), Ray et al. (ray01 (2001)), and Kaaret (kaa01 (2001)). The soft X-ray luminosity appears constant or slightly declining from epoch 16 yr through 20 yr. This result appears compatible with the model predictions of circumstellar interactions (Chevalier & Fransson che94 (1994)).

VLBI observations at 3.6, 6, 13 and 18 cm by Bartel et al. (1985), hereafter B85, determined a size of SN 1979C and its growth for epochs ranging from Dec 1982 to May 1984 assuming a synchrotron emitting optically thin spherical model. At those early epochs the source could not be resolved by any VLBI array. At an average expansion speed of, for example, $``$8 000 km s-1, one should expect a shell radius growth rate of $``$ 0.3$`\mu `$as day-1 which over 20 years should yield a shell radius of $``$2.2 mas which could be resolved by $`\lambda `$6 cm VLBI if detected on the long baselines. However, due to sensitivity limitations of the available arrays, the chances of detecting it at $`\lambda `$ 6 cm at the present time are null. On the other hand, at $`\lambda `$18 cm the emission strength of SN 1979C and the sensitivity of the arrays make the detection possible. We have carried out state-of-the-art $`\lambda `$18 cm observations in an attempt to learn as much as possible about the source structure.

## 2 Observations

SN 1979C was observed on 4 June 1999 from 19:30 UT to 01:30 UT at 1.6 GHz with a sensitive transatlantic VLBI array. The antennas used were (diameter, location): Effelsberg (100m, Germany), Robledo (70m, Spain), phased-VLA (130m-equivalent, NM, US), and Goldstone (70m, CA, US). Data were recorded using a standard VLBA mode with a data rate of 256Mb/s, covering a bandwidth of 64 MHz. The data were correlated at the VLBA Correlator Center of the National Radio Astronomy Observatory (NRAO)111NRAO is operated under cooperative agreement with the National Science Foundation.in Socorro, NM, US.

The observations were made with a duty cycle of about 26 min consisting of 3 min observations of the calibrator source PKS B1157+156, followed by 5 min observations of the calibrator source TXS 1214+161, in turn followed by 14 min observations of the source SN 1979C. The source 3C 274 was observed twice as a calibrator, once at the beginning and once at the end of the observations.

The correlated data were analyzed using the Astronomical Image Processing System (aips). A basic integration time of 2 sec was used. The phase slopes across each of the 8 MHz bands and over the 64 MHz synthesized band were aligned using a phase calibration determined from observations of the calibrator source 3C 274. The data of the calibrators PKS B1157+156 and TXS 1214+161 were then fringe-fitted in a standard manner. In the next step, phases, delays, and delay-rates corresponding to the sources PKS B1157+156 and TXS 1214+161 were interpolated to the source SN 1979C, to provide an initial model for the fringing of the data of this radio source. A solution was obtained by fringing the data using narrow search windows set around the initial model. Fringe detections of SN 1979C were obtained for all baselines at almost all times. However, some SN 1979C data had to be removed from the analysis after visual inspection of the data of the two main calibrators. Indeed, for some time ranges in the observations there were calibration problems or data with lower than acceptable quality. For instance, high winds at the phased-VLA produced an auto-stowing of some antennas and problems in the phasing during the second half of the observations. Therefore, data from the VLA from 23:30 UT to 01:15 UT were edited.

Adjacent to the VLBI observations reported in this paper, some cross-scans on PKS B1157+156 were carried out at 1.6 GHz to fix the pointing parameters of the antenna in Effelsberg. As a byproduct, those observations allowed the total flux density of PKS B1157+156 to be measured as 0.79$`\pm `$0.02 Jy (A. Kraus, priv.comm.). This value was adopted to set our flux density scale. Then, by readjusting the gains of the antennas, the VLBI map of PKS B1157+156 was forced to have the same flux density as the total flux density (that is, PKS B1157+156 was assumed unresolved by our VLBI array.) With this calibration procedure for the antennas of the array a total flux density of 5.4$`\pm `$0.2 mJy was obtained for SN 1979C.

Further model fitting analysis was carried out with modelfit in difmap (Shepherd et al. she94 (1994)), and with qfit in the Caltech VLBI Analysis Package (Pearson & Readhead pea84 (1984)).

## 3 The size of SN 1979C

We assumed spherical symmetry in the modeling. The data remaining after the editing were of the highest quality, but from a look at the data it was not fully clear that all points in the visibility vs. uv-radius plot are on the main lobe of the interferometric response. We had to prove to ourselves that such was indeed the case. Good measurements of the closure phases (all being around zero and none around 180) were a great help. We conducted several tests and model fits to convince ourselves that indeed the array at $`\lambda `$18 cm was not large enough to be able to resolve the structure of SN 1979C. Unfortunately, there are then many models which can reproduce the data equally well (see Marscher mar85 (1985) and B (85)).

Thus, we decided to determine the size of SN 1979C for the following putative source shapes: an optically thick source, an optically thin shell of width 30% of the external radius, and a ring that is an infinitely bright shell of zero width. The last one is obviously unphysical and its determination is only valuable as a limit. Even the first of the three putative source shapes is not expected given that the spectrum of the source is optically thin with $`S\nu ^{0.63\pm 0.03}`$ (Montes et al. mon00 (2000)). The diameters determined by model fitting from our data for such models are $`4.57\pm 0.25`$, $`3.60\pm 0.17`$, and $`3.10\pm 0.14`$ mas, respectively. The uncertainties were approximately estimated from the statistical errors given by the programs qfit and modelfit, taking into account discrepant results yielded by them.

A fit for the preferred model of an optically thin shell of diameter $`3.60\pm 0.17`$ mas with width 30% of the external radius is shown, along with the measured visibilities, as the solid line in Fig. 2.

## 4 Discussion

Marcaide et al. (1995a ; 1995b ; mar97 (1997)) determined for SN 1993J a shell structure with a shell width 30% of the outer radius, somewhat larger than that predicted by Chevalier & Fransson (che94 (1994)). Therefore, we have adopted a shell of width 30% of the external radius as a model for SN 1979C in this paper and all our discussion will be based on such model. For the present discussion all previous results by B (85) will be transformed (normalized) to this model.

At a distance of 16.1 Mpc to M100 an angular size of 1 mas corresponds to a linear size of $`2.4110^{17}`$ cm. Hence, the source outer radius (1.80 mas) corresponds to $`4.3310^{17}\mathrm{cm}`$. Since 20.12 yr have elapsed since the explosion, the corresponding average supernova expansion velocity is $``$6 800 km s-1. Since early expansion speeds were estimated at epoch 45 d from the H$`\alpha `$ lines as $``$9 200 km s-1 it follows that the expansion has been strongly decelerated by the action of the circumstellar matter. How strongly? It is difficult to make a detailed analysis to answer such a question, since intermediate size estimates are not available. Nonetheless, we now discuss our VLBI results, in conjunction with previously published VLBI and optical spectroscopy data, to show that SN 1979C has indeed suffered a heavy deceleration in its expansion, starting $`6\pm 2`$ yr after its explosion.

B85 determined for SN 1979C a diameter of $`2.16_{0.50}^{+0.44}`$ mas 5.17 yr after explosion, from 1.67 GHz VLBI observations using an optically thin uniform sphere model. This value translates into a diameter of $`1.56_{0.36}^{+0.31}`$ mas for our discussion shell model. This size corresponds to an average speed for the first 5.17 yr of $``$11 500 km s-1. On the other hand, the average speed of the supernova over the first six weeks after explosion as determined from the widths of H$`\alpha `$ lines was $`\mathrm{9\hspace{0.17em}200}\pm 500`$ km s-1 (Panagia et al. pan80 (1980), G. Vettolani, priv. comm.), which is about 20% lower than the average velocity estimated for the first 5 years from the B (85) estimates normalized to our discussion model. These two results are nicely consistent. Indeed, in the mini-shell model of Chevalier che82 (1982), the radio emission comes from the forward shock, while the optical, UV, and X-ray emissions come mainly from the reverse shock, which moves slower than the forward shock. Hence the velocity ratio determined from the VLBI and optical observations is roughly what is expected for the radii ratio of the forward and reverse shocks, in a model of 20-30% shell width (Chevalier & Fransson che94 (1994), Marcaide et al. mar97 (1997)) if SN 1979C followed a free, self-similar expansion for the first 5 years after its explosion. A free expansion for the first 4 years, followed by a decelerated expansion would also be compatible with the mentioned results. These results should be understood as evidence of free expansion of SN 1979C for the first 4-5 years and of the presence of a shell-like structure, and not as evidence of a given shell width, since different model normalizations of the results of B85 would yield equally consistent results and a similar conclusion.

An estimate of $`m`$ ($`Rt^m`$) using the size estimate of B (85) and our size estimate results in $`m=0.62_{0.17}^{+0.22}`$. This results shows that SN 1979C has indeed dramatically decreased its growth rate after a likely initial phase of free expansion ($`m1.0`$; see Fig. 3). The optical data, although as scarce as the VLBI data, can be used as a check of the reliability of our estimate of $`m`$. In a self-similar expansion scenario $`Rt^m`$ and, consequently, $`vt^{m1}`$. If we now assume that the supernova did not decelerate for the first 5.17 years, we can take the value of 9 200 km s-1 as the expansion velocity just prior to the phase of strongly decelerated expansion. Combining this value with that of 6 200$`\pm `$300 km s-1 at epoch 14.09 yr (F99 and R. Fesen, priv. comm.), gives $`m=0.61_{0.11}^{+0.10}`$, in agreement with our estimate of $`m=0.62_{0.17}^{+0.22}`$.

Hence, $`m=0.62`$ appears to be an adequate value to characterize a phase of strong deceleration in the expansion of SN 1979C, assuming it started about 5.17 years after its explosion. If the supernova started to decelerate later than 5.17yr, then $`m`$ should be smaller, and correspondingly stronger the deceleration of SN 1979C. It seems somewhat arbitrary to assume that the deceleration set in at the same time as the observations of B85. Indeed, there is nothing apparent in the integrated radio light curves which could point to a change at that particular epoch. It is intriguing, though, that the integrated radio light curves of SN 1979C (Fig. 1) seem to show an increase in flux at around t=2200 days, or 6 yr after its explosion. One could speculate that the flux rise at this epoch is somehow related to a change in the deceleration parameter, likely due to the shock front of the supernova entering a high-density region of the circumstellar medium (CSM).

Our estimate of the angular diameter for SN 1979C can also be used to obtain an upper limit for the epoch at which the phase of strong deceleration started. We find that the latest epoch of free expansion which is still compatible with both the optical observations of F99 (epoch 14.09yr) and our VLBI determination at 20.12 yr, and with the fact that a lower limit is set by SN 1979C entering a Sedov phase ($`m=0.4`$), corresponds to $``$8 yr.

The integrated radio emission curves have been well modeled for the first 10 years by Weiler et al. (wei92 (1992)) and Montes et al. (mon00 (2000)) who assumed self-similar expansion and a modulation in the density of the CSM due to the binary interaction, to explain the modulation in the flux density curve. Montes et al. (mon00 (2000)) found that starting around epoch 10 yr the radio emission started to grow after a long decline predicted by the radio emission model. Can our VLBI observations shed any light on what happened at epoch 10 yr? The VLBI results by B (85) and us require a strong deceleration of the expansion. Such deceleration is compatible with the optical results by F99 for a self-similar expansion. Our results rule out any initial free expansion longer than 8 years. Thus, the onset of decelerated expansion took place well before the epoch 10 yr where the radio emission decline turned onto a slight growth. The sooner the onset of the decelerated expansion the milder the deceleration needed to be compatible with the observations. The epoch of that onset is restricted to the range 4-8 years and the corresponding deceleration parameter $`m`$ to the range 0.67-0.44. (Strictly speaking an onset at 3 years with a corresponding $`m`$=0.72 is not excluded, but it appears unlikely since it yields a radio to optical velocity ratio of 0.68 at epoch 5.17 yr. Such a ratio is small, but not yet inconsistent with the models used. In any case, the deceleration required for an early onset is comparable to the strongest decelerations reported in other cases by Marcaide et al. (mar97 (1997)) for SN 1993J in M81 and Mc Donald et al. (mcd01 (2001)) for 43.31+592 in M82). An onset of decelerated expansion at epoch 6 yr (corresponding to a deceleration parameter $`m`$=0.57) appears most likely and it may be associated with a small halt in the otherwise decreasing flux density trend, a halt which may be partly related to encountering a denser CSM.

If the expansion of SN 1979C has significantly decelerated, as our VLBI estimate shows, then the mass of the CSM swept up by the shock front, $`M_{\mathrm{sw}}`$, must be comparable to or larger than the mass of ejected hydrogen-rich envelope, $`M_{\mathrm{env}}`$. Our estimated angular diameter translates into a forward shock radius of $`4.33\times 10^{17}`$ cm. For a presupernova wind of velocity $`v_\mathrm{w}`$=10 km s-1, this radius corresponds to $``$13 700 yr prior to the explosion of SN 1979C. Weiler et al. (1991) estimated a mass-loss rate of $`1.2\times 10^4M_{\mathrm{}}`$ yr-1 for the progenitor of SN 1979C. If we now assume that the circumstellar medium up to this distance is a standard one ($`\rho _{\mathrm{cs}}r^2`$), we then obtain $`M_{\mathrm{sw}}1.6M_{\mathrm{}}`$. Momentum conservation therefore implies that $`M_{\mathrm{env}}0.59M_{\mathrm{sw}}`$, which results in a mass for the ejected hydrogen-rich envelope of $`0.94M_{\mathrm{}}`$.

In the case of the Type II supernova SN1993J, whose progenitor was a $`15M_{\mathrm{}}`$ star in a binary system, mass exchange with its companion meant that the mass of the hydrogen-rich envelope ejected at the explosion was $`0.20.4M_{\mathrm{}}`$ (Woosley et al. 1994, Houck & Fransson 1996), much less than the $`3.3M_{\mathrm{}}`$ one would expect from a single-star model (Höflich et al. 1993). Weiler et al. (1986) have suggested that SN1979C was born in a binary system consisting of a 15-18 $`M_{\mathrm{}}`$ red supergiant and a 10 $`M_{\mathrm{}}`$ B1 main-sequence star. If we assume that SN 1979C lost, like SN 1993J, most of its hydrogen-rich envelope prior to explosion, then the $`0.9M_{\mathrm{}}`$ decelerated ejecta suggested from our observations represents at least all of the remaining hydrogen-rich envelope.

¿From VLBI observations of SN 1979C more than 20 years after its explosion, we conclude that an initial free expansion of SN 1979C was followed by a strong deceleration which set on $``$6 yr after the supernova explosion. We estimate the shocked swept-up mass to be $`M_{\mathrm{sw}}1.6M_{\mathrm{}}`$, assuming a standard density profile for the CSM. Such a swept-up mass for SN 1979C suggests a mass of the hydrogen-rich envelope ejected at explosion as low as $`M_{\mathrm{env}}0.94M_{\mathrm{}}`$ or lower. If, as suggested by Weiler et al. (wei86 (1986)) and Schwarz & Pringle (sch96 (1996)) SN 1979C originated in a binary star (alike SN 1993J) this value of $`M_{\mathrm{env}}`$ suggests that the companion of the progenitor star stripped off most of the hydrogen-rich envelope mass of the presupernova star prior to its explosion. New observations at $`\lambda `$18 cm around 2005 should be able to test the above conclusion and would also be essential to image the true structure of SN 1979C. On the other hand, a continued rise in the level of flux density at $`\lambda `$6 cm and improved VLBI instrumentation and techniques like phase-referencing may allow an earlier determination of the intrinsic structure of SN 1979C with high resolution and a test of the binary star scenario which the present results tend to favor.

###### Acknowledgements.

We are grateful to A. Kraus for the help provided during the observations and the flux density measurements on PKS B1157+156. This work has been partially supported by the Spanish DGICYT Grants No. PB96-0782 and PB97-1164. NRAO is operated under license by Associated Universities, Inc., under cooperative agreement with NSF. We acknowledge support from the European Commission’s TMR-LSF programme, Contract No. ERBFMGEST950012. JMM is grateful to the Max-Planck Institut für Radioastronomie for the support received for this research during his visit. KWW wishes to thank the Office of Naval Research for the 6.1 funding supporting this work.


text:
# The X-ray Spectra of VY Scl Stars Are Not Blackbodies

## 1 Introduction

Cataclysmic variables (CVs) are a diverse class of semidetached binaries including novae, dwarf novae, and novalike variables, composed typically of a low-mass main-sequence secondary and a white dwarf. With the exception of novae in outburst, the engine for all CVs is the release of gravitational potential energy as material accretes onto the white dwarf. In nonmagnetic systems accretion is mediated by a disk, and simple theory predicts that half of the gravitational potential energy of the accreting material is liberated in the disk and half is liberated in the boundary layer between the disk and the surface of the white dwarf, with luminosities $`L_{\mathrm{disk}}L_{\mathrm{bl}}GM_{\mathrm{wd}}\dot{M}/2R_{\mathrm{wd}}=4\times 10^{34}(\dot{M}/10^8\mathrm{M}_{}\mathrm{yr}^1)(M_{\mathrm{wd}}/\mathrm{M}_{})(R_{\mathrm{wd}}/10^9\mathrm{cm})^1\mathrm{erg}\mathrm{s}^1`$, where $`\dot{M}`$ is the mass-accretion rate and $`M_{\mathrm{wd}}`$ and $`R_{\mathrm{wd}}`$ are respectively the mass and radius of the white dwarf. When $`\dot{M}`$ is low (e.g., dwarf novae in quiescence), the boundary layer is optically thin and quite hot (of order the virial temperature $`T_{\mathrm{vir}}=GM_{\mathrm{wd}}m_\mathrm{H}/6kR_{\mathrm{wd}}20`$ keV); when $`\dot{M}`$ is high (e.g., novalike variables and dwarf novae in outburst), the boundary layer is optically thick and quite cool (of order the blackbody temperature $`T_{\mathrm{bb}}=[GM_{\mathrm{wd}}\dot{M}/8\pi \sigma R_{\mathrm{wd}}^3]^{1/4}10`$ eV). In the context of this theory, high-$`\dot{M}`$ CVs are X-ray sources only because the upper “atmosphere” of the boundary layer remains optically thin.

The observational picture of the X-ray emission of CVs was first sketched in the 1970s and has become clearer through the years as instruments with higher effective area, broader bandpass, and better spectral resolution have flown on X-ray satellites. Summaries of HEAO 1, Einstein, EXOSAT, and ROSAT investigations are provided by Córdova, Jensen, & Nugent (1981); Córdova & Mason (1984); Eracleous, Halpern, & Patterson (1991); Mukai & Shiokawa (1993); Richman (1996); van Teeseling, Beuermann, & Verbunt (1996); and Verbunt et al. (1997). That the X-ray emission region of nonmagnetic CVs is compact and centered on the white dwarf is established directly by the X-ray light curves of eclipsing systems. That the X-ray spectra of CVs are due to hot thermal plasma was established by ASCA, which had the combination of large effective area, broad bandpass, and high spectral resolution needed to resolve the emission lines of K-shell Mg, Si, S, Ar, and Fe at high energies, and to establish the presence of emission lines of K-shell O and Ne and L-shell Fe at low energies (Nousek et al., 1994; Mukai, 2000; Baskill, Wheatley, & Osborne, 2001). As evidenced by EUVE light curves and spectra of the dwarf novae SS Cyg, U Gem, VY Hyi, and OY Car in outburst, the blackbody component of high-$`\dot{M}`$ nonmagnetic CVs is typically not observed in the canonical X-ray bandpass (SS Cyg in outburst is the only clear exception, see Beuermann & Thomas 1993; Ponman et al. 1995; Mauche, Raymond, & Mattei 1995) because its temperature is too low and its luminosity is sometimes anomalously weak (Mauche, 2002, and references therein). While it is not yet possible to be entirely certain, the working hypothesis is that the X-ray spectra of all nonmagnetic CVs are due to hot thermal plasma in the boundary layer between the accretion disk and the surface of the white dwarf.

## 2 VY Scl Stars

Given this understanding, it is unsettling that the X-ray spectra of VY Scl stars have been repeatedly described as absorbed blackbodies. VY Scl stars are a class of novalike variables that occasionally dim by 3–5 mag in the optical; they are sometimes referred to as “anti-dwarf novae.” It is thought that the downward transitions of VY Scl stars are the result of temporary reductions in the mass-transfer rate from the secondary, possibly due to the passage of star spots over the $`L_1`$ point. VY Scl stars are more than oddities in the CV zoo. Livio & Pringle (1994) have proposed that the same mechanism that causes the intensity dips in VY Scl stars (with $`P_{\mathrm{orb}}=3`$–4 hr) ultimately results in the $`P_{\mathrm{orb}}=2`$–3 hr gap in the orbital period distribution of CVs. Leach et al. (1999) have proposed that the lack of dwarf nova outbursts in VY Scl stars in quiescence is the result of irradiation of the inner disk by the hot ($`T\genfrac{}{}{0pt}{}{\text{ }>\text{ }}{\text{ }\text{ }}40`$ kK) white dwarfs known to exist in these systems. The white dwarfs in novalike variables should be similarly hot because of their relatively high average mass-accretion rates. Otherwise, novalike variables, VY Scl stars in their high states, and dwarf novae in outburst should all be quite similar. In particular, we expect that their X-ray spectra are produced by hot thermal plasma in the boundary layer between the accretion disk and the surface of the white dwarf.

The first indication that the X-ray spectra of VY Scl stars might differ from other high-$`\dot{M}`$ CVs was provided by Schlegel & Singh (1995), who found that the ROSAT PSPC spectra of the VY Scl stars MV Lyr and KR Aur are best described by absorbed blackbodies. Next, van Teeseling, Beuermann, & Verbunt (1996) published a figure showing that the ROSAT PSPC colors of the VY Scl stars TT Ari, KR Aur, BZ Cam, and MV Lyr lay in a region of parameter space inhabited by blackbodies (and highly absorbed thermal plasmas). Finally, Greiner (1998) concluded from a comprehensive study of ROSAT pointed and all-sky survey PSPC observations that the X-ray spectra of VY Scl stars in their high optical states are best described by absorbed blackbodies. Greiner fit the ROSAT PSPC spectrum of the brightest VY Scl star, TT Ari, with absorbed power law, thermal bremsstrahlung, Raymond-Smith thermal plasma, and blackbody models and obtained $`\chi _\nu ^2=2.23`$, 1.99, 1.87, and 1.67, respectively. Although noting that even the absorbed blackbody model fit was poor, Greiner advocated the blackbody interpretation of the ROSAT PSPC spectra of VY Scl stars as a class because it resulted in similar parameters for the brightest members of the class: blackbody temperatures $`kT_{\mathrm{bb}}0.25`$–0.5 keV and emitting sizes of 50–120 m (hence fractional emitting areas $`f=\mathrm{Area}/4\pi R_{\mathrm{wd}}^2\genfrac{}{}{0pt}{}{\text{ }<\text{ }}{\text{ }\text{ }}1\times 10^{11}`$).

There are a number of problems with the conclusion that the X-ray spectra of VY Scl stars are blackbodies. First, it runs counter to our theoretical expectation that the X-ray emission of nonmagnetic CVs is due to hot thermal plasma. Second, the implied emitting areas are far smaller than any structure known in CVs. The area of the boundary layer $`A_{\mathrm{bl}}2\pi R_{\mathrm{wd}}\mathrm{\hspace{0.17em}2}H_{\mathrm{bl}}`$, where $`H_{\mathrm{bl}}`$ is the vertical height of the boundary layer. As argued by Patterson & Raymond (1985), a lower limit to $`H_{\mathrm{bl}}`$ is the density scale height of the inner disk $`H_{\mathrm{disk}}=3\times 10^6(\dot{M}/10^8\mathrm{M}_{}\mathrm{yr}^1)^{0.18}(M_{\mathrm{wd}}/\mathrm{M}_{})^{1.2}`$ cm (Pringle & Savonije, 1979). Hence, the boundary layer fractional emitting area $`f\genfrac{}{}{0pt}{}{\text{ }>\text{ }}{\text{ }\text{ }}H_{\mathrm{disk}}/R_{\mathrm{wd}}=3\times 10^3(\dot{M}/10^8\mathrm{M}_{}\mathrm{yr}^1)^{0.18}(M_{\mathrm{wd}}/\mathrm{M}_{})^{1.2}(R_{\mathrm{wd}}/10^9\mathrm{cm})^1`$. It is theoretically possible to produce even smaller accretion spots by magnetically channeling the flow of material down to the white dwarf surface, but even in CVs with the highest magnetic field strengths (AM Her stars or polars, with $`B10`$–100 MG), the fractional emitting area $`f10^3`$ (Mauche, 1999). If the spot sizes were this small in VY Scl stars, it would be hard to avoid flux modulations at the white dwarf spin period. Third, the blackbody temperatures are too high for a white dwarf accretor. With $`kT_{\mathrm{bb}}0.25`$–0.5 keV, the implied local energy flux $`\sigma T_{\mathrm{bb}}^4=4`$$`64\times 10^{21}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$, orders of magnitude greater than the Eddington rate: $`L_{\mathrm{edd}}/4\pi R_{\mathrm{wd}}^2=1\times 10^{19}(M_{\mathrm{wd}}/\mathrm{M}_{})(R_{\mathrm{wd}}/10^9\mathrm{cm})^2\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$.

## 3 ASCA Observations

Given the apparent discrepancy between ROSAT observations and theory, we searched the ASCA archive for data to test the hypothesis that the X-ray spectra of VY Scl stars in their high optical states are blackbodies. We found two VY Scl stars in the archive: TT Ari and KR Aur, observed by ASCA on 1994 January 20/21 and 1996 March 6, respectively. By examining the AAVSO records, we determined that on the dates of the ASCA observations TT Ari was in a high optical state with $`V10.5`$, while KR Aur was in an intermediate brightness state with $`V16`$ (whereas in its high and low states $`V14`$ and $`V19`$, respectively). Therefore, the only high-state VY Scl star with data in the ASCA archive is TT Ari, but this is the brightest VY Scl star, so provides the most stringent test of the nature of the X-ray spectra of these high-$`\dot{M}`$ nonmagnetic CVs. For completeness, we also analyzed the ASCA spectrum of KR Aur.

ASCA carries four co-aligned X-ray telescopes, two with gas scintillation imaging proportional counter (Gas Imaging Spectrometer; GIS) detectors and two with CCD (Solid-state Imaging Spectrometer; SIS) detectors (Tanaka, Inoue, & Holt, 1994). During both observations, all four instruments were operating normally, with the GIS detectors in standard PH mode. We screened the GIS data by excluding intervals during Earth occultations (up to $`5^{}`$ above the horizon), passages through the South Atlantic Anomaly (SAA), and times of high background as indicated by the in-situ monitors. After screening, we obtained 18.5 ks and 21.6 ks of good GIS data for TT Ari and KR Aur, respectively. During the TT Ari observation, the SIS detectors were operated in 2-CCD clocking mode, switching between Faint and Bright data mode depending on the available telemetry rate. We applied a similar set of screening criteria (only differing in the details of the monitor counts used) to the SIS data, and further excluded data taken within $`20^{}`$ of the bright Earth limb and 64 s of day/night and SAA transitions. After the screening, we obtained 17.5 ks of good SIS data for TT Ari. During the KR Aur observation, the SIS detectors were used in 1-CCD Faint mode throughout, and after a similar screening (but with $`15^{}`$ limit for the bright Earth limb and 32 s for day/night and SAA transitions), we obtained 21.1 ks of good SIS data.

Next, we extracted source and background spectra for each detector for each observation. For TT Ari, source counts were extracted from circular regions centered on the source ($`6^{}`$ radius for the GIS, $`2.5^{}`$$`3.5^{}`$ radius for the SIS to stay within a single CCD chip), while background counts were extracted from adjoining source-free regions. Similar extraction regions could not be used for KR Aur because of the presence of faint sources near the VY Scl star. We therefore reduced the GIS source extraction region to $`4^{}`$ radius, and used a local background region that avoids both KR Aur and the contaminating sources. For the SIS, we used a similar source extraction region as for TT Ari, but extracted background spectra from blank-sky observations, since a source-free local background region of sufficient size could not be defined. For the spectral responses of the instruments, we have taken GIS redistribution matrix v4.0, and generated SIS redistribution matrix and auxiliary response files using the standard (HEAsoft v5.1) software. For each observation, we co-added spectra and responses for the two GIS instruments into one set, and the two SIS instruments into another set. Finally, the SIS and GIS data sets were fit simultaneously after ignoring bad channels, grouping the channels to produce a minimum of 20 counts per channel, and restricting the energy range to 0.5–10 keV. Each data set was fit with three spectral models: a single-temperature blackbody, a single-temperature solar-abundance MEKAL thermal plasma, and a two-temperature solar-abundance MEKAL thermal plasma, all absorbed by a column density of neutral material. To account for known residual discrepancies between the SIS and GIS response functions, we allowed the absorbing column density to differ between the SIS and GIS.

The best-fit parameters (and 90% confidence intervals for the thermal plasma models) resulting from these fits are shown in Table 1. The single-temperature blackbody model produces very poor fits ($`\chi _\nu ^2=5.77`$ and 2.24) to the 0.5–10 keV ASCA spectra of TT Ari and KR Aur because a blackbody produces too few counts both at high and low energies. Even worse fits are realized with the blackbody parameters $`kT_{\mathrm{bb}}0.3`$–0.5 keV and $`N_\mathrm{H}1`$$`2\times 10^{20}\mathrm{cm}^2`$ favored by fits to the ROSAT PSPC spectra of these stars. Far better fits ($`\chi _\nu ^2=1.31`$ and 0.75) are obtained with the same number of degrees of freedom with the single-temperature thermal plasma model. In particular, the Fe K$`\alpha `$ line is clearly detected near 6.7 keV in the GIS data of both sources, which is decisive evidence for the presence of a hot ($`kT5`$–10 keV) plasma. With two additional degrees of freedom, the two-temperature thermal plasma model further improves the fits—by a significant amount in TT Ari ($`\mathrm{\Delta }\chi ^2=116`$, $`\chi _\nu ^2=1.07`$), and by a modest amount in KR Aur ($`\mathrm{\Delta }\chi ^2=22`$, $`\chi _\nu ^2=0.66`$)—by better accounting for systematic residuals at low energies. The best-fit two-temperature thermal plasma model is shown in Figure 1 superposed on the SIS and GIS spectra of TT Ari.

## 4 Summary

Based on an analysis of ASCA SIS and GIS spectra, we find that the 0.5–10 keV X-ray spectra of TT Ari and KR Aur are poorly fit by an absorbed blackbody model but are well fit by an absorbed thermal plasma model. The TT Ari spectra are adequately described by a two-temperature solar-abundance thermal plasma, with one component at $`kT7`$ keV and another at $`kT0.7`$ keV, with a relative emission measure of approximately 20:1 in favor of the high-temperature component. The lower-quality KR Aur spectra are adequately described by a single-temperature solar-abundance thermal plasma with $`kT6`$ keV, and allow a second temperature component only if the column density is significantly increased and the emission measures of the two components are comparable. These results should be understood to be simply parameterizations of the X-ray spectra of these two stars, not a definitive determination of the nature of their X-ray spectra. Higher-quality data (e.g., Chandra High-Energy Transmission Grating and XMM Reflection Grating Spectrometer and European Photon Imaging Camera spectra) are required to determine—by resolving the L-shell emission lines of Fe and the K-shell emission lines of H- and He-like ions of abundant elements from C to Fe—the true emission measure distribution, abundances, density, and Doppler broadening of the plasma in these stars. It is clear, however, that the X-ray spectra of neither of these stars is that of a blackbody, contrary to the conclusions of Schlegel & Singh (1995) and Greiner (1998). The different conclusions about the nature of the X-ray spectrum of KR Aur may be due to differences in the accretion rate, since this star was in a high optical state during the ROSAT observation, but in an intermediate optical state during the ASCA observation. TT Ari, on the other hand, was in a high optical state during both observations, so directly contradicts the claim that the X-ray spectra of VY Scl stars in their high optical states are blackbodies. It remains possible that TT Ari, the brightest VY Scl star, is anomalous, and that the X-ray spectra of other high-state VY Scl stars are blackbodies, but based on theoretical expectations and the ASCA, Chandra, and XMM X-ray spectra of other high-$`\dot{M}`$ nonmagnetic CVs, we believe that they are not. Instead, we believe that the X-ray spectra of VY Scl stars in their low and high optical states are due to hot thermal plasma in the boundary layer between the accretion disk and the surface of the white dwarf, and appeal to the acquisition of Chandra and XMM grating spectra to test this prediction.

In this research, we have used data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA’s Goddard Space Flight Center. We have also used, and acknowledge with thanks, data from the AAVSO International Database, based on observations submitted to the AAVSO by variable star observers worldwide. C. W. M.’s contribution to this work was performed under the auspices of the U.S. Department of Energy by University of California Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48.