text:
# Finite-size-scaling analyses of the chiral order in the Josephson-junction ladder with half a flux quantum per plaquette

## 1 Introduction

By means of a field-theoretical technique, Nersesyan, Gogolin and Eßler claimed Nersesyan98  that the ground state of the $`XXZ`$ spin chain with sufficiently strong next-nearest-neighbor interaction is in the chiral phase, where the spin-screw chirality is broken spontaneously. Their proposal is astonishing, because the $`XXX`$ counterpart, that has been studied very extensively so far Haldane82 ; Tonegawa87 ; Harada88 ; White96 ; Allen97 ; Harada99 , is known to be in the dimer phase, where the translational symmetry is broken spontaneously so long as the frustration is strong enough. New treatments devised specifically for $`XXZ`$ were reported Cabra00 ; Allen99  so as to support the aforementioned prediction. Meanwhile, as for other frustrated $`XXZ`$ chain, namely, the three-leg ladder with diagonal interchain interaction, Azaria et al. argued Azaria98 ; Azaria99  that a unique critical phase possessing high central charge ($`c=2`$) might be realized in the ground state. According to their argument, the chiral-symmetry breaking is significant to stabilize such high-$`c`$ criticality. These recent developments tell that the $`XXZ`$ anisotropy with in-plane frustration may give rise to unexpected exotic phases. The above-mentioned proposals are all based on field-theoretical descriptions. Therefore, in order to support these scenarios, numerical simulations should be carried out. Kaburagi et al. employed the exact diagonalization and the density-matrix renormalization-group methods for the $`XX`$ model with the next-nearest-neighbor coupling Kaburagi99 ; Hikihara00 . They concluded that the chirality is not broken for the $`S=1/2`$ chain, whereas it would be broken for the $`S=1`$ chain; namely, numerical-simulation result for $`S=1/2`$ appears to contradict the field-theoretical description. Nevertheless, numerical simulation of such frustrated chain itself is a matter of serious methodological concern White96 ; main obstacles are due to the exponentially small energy gap and incommensurability of the spin-correlation function, which often miss-matches the total system sizes.

In order to avoid such complications that would arise from incommensurability, we have investigated the two-leg Josephson-junction ladder,

$``$ $`=`$ $`{\displaystyle \frac{t_{}}{2}}{\displaystyle \underset{i=1}{\overset{L}{}}}(\mathrm{e}^{\mathrm{i}\mathrm{\Phi }}a_{1i}^{}a_{1,i+1}+\mathrm{e}^{\mathrm{i}\mathrm{\Phi }}a_{1,i+1}^{}a_{1i})`$ (1)
$`{\displaystyle \frac{t_{}}{2}}{\displaystyle \underset{i=1}{\overset{L}{}}}(a_{2i}^{}a_{2,i+1}+a_{2,i+1}^{}a_{2i})`$
$`{\displaystyle \frac{t_{}}{2}}{\displaystyle \underset{i=1}{\overset{L}{}}}(a_{1i}^{}a_{2i}+a_{2i}^{}a_{1i})`$
$`V_{}{\displaystyle \underset{l=1,2}{}}{\displaystyle \underset{i=1}{\overset{L}{}}}a_{li}^{}a_{li}a_{l,i+1}^{}a_{l,i+1}`$
$`V_{}{\displaystyle \underset{i=1}{\overset{L}{}}}a_{1i}^{}a_{1i}a_{2i}^{}a_{2i},`$
with gauge-twist angle,

$$\mathrm{\Phi }=\pi .$$
(2)
Here, $`a_{li}`$ and $`a_{li}^{}`$ are the hard-core boson annihilation and creation operators, respectively, at the $`i`$-th site along the $`l`$-th leg. The parameters $`t_{}`$, $`t_{}`$, $`V_{}`$ and $`V_{}`$ control the intrachain hopping amplitude, the interchain hopping amplitude, the nearest-neighbor intrachain attraction and the interchain attraction, respectively. The number of bosons is set to $`N=L`$ (half filled). Throughout this paper, we choose $`t_{}`$ as the unit of energy; namely,

$$t_{}=1.$$
(3)
We have imposed the periodic-boundary condition along the ladder; $`a_{l,L+1}=a_{l,1}`$. The angle $`\mathrm{\Phi }`$ denotes the gauge twist around each plaquette. Because the angle is set to $`\mathrm{\Phi }=\pi `$, the ladder is subjected to a uniform magnetic field of half a flux quantum per plaquette. That is, suppose that the bosons are in the superconducting state (this is a subtle issue in one dimension), the magnetic flux would possibly be quantized so as to form rigid vortex-lattice order along the ladder. This order is confirmed to develop in the preceding semi-classical analysis Ciria99 . That is, the chiral order of our model has two-unit-cell periodicity. For the purpose of studying the stability of chirality against the quantum fluctuation numerically, our model is far more advantageous than the frustrated spin chain, because the spin-screw chiral order of the latter model has long-wavelength incommensurate structure, and thereby the numerical data suffer from insystematic finite-size-scaling behavior. Tendency to the formation of such vortex structure may become more transparent, if we transform the boson Hamiltonian (1) in terms of spin: Through utilizing the mapping relations between the hard-core boson and the spin-$`1/2`$ operators, namely, $`a_{li}=S_{li}^{}`$, $`a_{li}^{}=S_{li}^+`$ and $`a_{li}^{}a_{li}1/2=S_{li}^z`$, the above boson model is mapped to the $`XXZ`$ spin ladder,

$``$ $`=`$ $`t_{}{\displaystyle \underset{i=1}{\overset{L}{}}}(S_{1i}^xS_{1,i+1}^x+S_{1i}^yS_{1,i+1}^y)`$ (4)
$`t_{}{\displaystyle \underset{i=1}{\overset{L}{}}}(S_{2i}^xS_{2,i+1}^x+S_{2i}^yS_{2,i+1}^y)`$
$`t_{}{\displaystyle \underset{i=1}{\overset{L}{}}}(S_{1i}^xS_{2i}^x+S_{1i}^yS_{2i}^y)`$
$`V_{}{\displaystyle \underset{l=1,2}{}}{\displaystyle \underset{i=1}{\overset{L}{}}}S_{li}^zS_{l,i+1}^zV_{}{\displaystyle \underset{i=1}{\overset{L}{}}}S_{1i}^zS_{2i}^z,`$
apart from a constant term. Note that the signs of the $`XY`$-component magnetic interactions are different from one leg and the other, and thus an in-plane frustration does exist. Therefore, vortex-antivortex alignment gets favored. The semi-classical version $`S\mathrm{}`$ of the model (4) has been studied extensively Granato92 ; Granato93 ; Denniston95 ; Granato96 ; Ciria99 : In the limit, spins are represented by rotators, and through path-integral mapping, the system reduces to a two-dimensional classical rotator model. According to numerical-simulation studies Granato92 ; Granato93 ; Granato96 , that rotator model appears to be in the KT critical phase. To our surprise, surrounded by such critical background of in-plane spin components, rigid long-range chiral order develops. And so, it is our motivation to investigate the stability of the chiral order against the quantum fluctuations of $`S=1/2`$. We show in this paper, that the chiral order fails to develop except in a limited condition, where an attractive intrachain coupling and a repulsive interchain coupling are both turned on. This new phase is a demonstration of critical phase accompanying chiral-symmetry breaking predicted field-theoretically Nersesyan98 ; Cabra00 ; Allen99 .

Let us mention some remarks concerning the present ladder model described by either boson representation eq. (1) or spin representation eq. (4). One might wonder that the hard-core condition, in other words, the $`S=1/2`$ condition, would be too restrictive, and the model misses microscopic physical ingredients such as the intra-and-inter-grain charge capacitances. Actually, setup of our model might be rather phenomenological. We stress, however, that the one-dimensional $`XXZ`$ chains belong to the Tomonaga-Luttinger-liquid universality class irrespective of $`S`$ Alcaraz92 . (Heuristically, it is known that one-dimensional quantum systems possessing the $`\mathrm{U}(1)`$ symmetry are flown to that universality in many cases.) As a matter of fact, according to the proposal Alcaraz92 , the transverse-correlation exponent $`\eta _{}`$ is governed by the compact formula,

$$\eta _{}=\frac{\pi \mathrm{cos}^1J^z}{2\pi S}.$$
(5)
(In our notation, $`J^z`$ is given by $`J^z=V_{}`$.) That is, through varying $`J^z`$, one can cover all possible range of $`\eta _{}`$ with $`S=1/2`$ fixed. Lastly, we mention the possibility of the so-called ‘spin liquid’ state that often arises in quantum spin systems; here, the term ‘spin liquid’ denotes the state with exponentially decaying short-range correlation function. we use the term spin liquid that The conventional two-leg ladder Dagotto92 ; Dagotto96  (with non-frustrated interchain coupling) is a prototypical system exhibiting spin-liquid state. We will show that our frustrated ladder does exhibit spin-liquid phase for a certain parameter condition beside the chiral phase.

This paper is organized as follows. In the next section, we explore the ladder model by means of the exact diagonalization method. The last section is devoted to summary and discussions.

## 2 Numerical results

In this section, we present numerical results. We carried out exact-diagonalization calculations for the system (4) with up to $`N=2L=32`$ spins. The data are analyzed in terms of the finite-size-scaling theory.

In Fig. 1, we have plotted the Binder parameter Binder81a ; Binder81b ,

$$U=1\frac{M^4}{3M^2^2},$$
(6)
for the chiral order,

$`M=M_{\mathrm{chiral}}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{L}{}}}(1)^i\left[𝐒_{1i}\times 𝐒_{2i}\right]_z`$ (7)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{L}{}}}(1)^i(S_{1i}^xS_{2i}^yS_{1i}^yS_{2i}^x)`$
$`=`$ $`{\displaystyle \frac{1}{2\mathrm{i}}}{\displaystyle \underset{i=1}{\overset{L}{}}}(1)^i(a_{1i}^{}a_{2i}+a_{1i}a_{2i}^{}),`$
for the system (1) with $`t_{}=0.5`$, $`V_{}=0.6`$ and $`V_{}`$ varied.

Here, $`\mathrm{}`$ denotes the ground-state average. Note that in the boson language, $`M_{\mathrm{chiral}}`$ measures the staggered boson current through the rungs. Therefore, it detects the formation of the vortex-lattice order. Finite-size scaling behavior of the Binder parameter contains the information whether the order $`M`$ develops or not Binder81a ; Binder81b : If the order is long (short) ranged, the Binder parameter grows (becomes suppressed) through enlarging the system sizes. At critical point, the Binder parameter remains scale-invariant. From Fig. 1, we found that the chiral order develops for the parameter range $`1.8V_{}0`$, while it is short ranged otherwise. As would become more apparent in the succeeding analyses, the strong-coupling regions of $`V_{}2.5`$ and $`V_{}1.8`$ belong to insulator phases of different characters, and are thus rather out of present concern. For $`V_{}2.5`$, for instance, the bosons are so cohesive that they may constitute an island in the sea of vacant sites; namely, the system becomes ‘phase separated.’ On the other hand, for $`V_{}1.8`$, owing to the strong repulsion, a Mott gap opens between the bonding and anti-bonding excitation branches.

Our result shows that the region of the chiral phase is limited. Hence, we found that quantum fluctuations are dominating so that unlike the semi-classical case, assistance of the many-body correlations of $`V_{}`$ and $`V_{}`$ is vital for stabilizing the chiral order. After scanning the parameter space, we found that the optimal condition lies around $`t_{}=0.5`$, $`V_{}=0.6`$ and $`V_{}=1`$. It is fairly reasonable that the optimal intrachain interaction is attractive, because attractive interaction enhances the tendency toward superconductivity ($`XY`$ order); see eq. (5), for instance. On the contrary, the fact that the repulsive interchain interaction is favorable sounds astonishing. It is expected that the repulsive interchain coupling might enhance the particle exchange between the chains. Hence, we see that such particle exchange is significant in order to confine (pin) the magnetic fluxes at each plaquette.

In order to confirm the above phase diagram, we have calculated the scaled domain-wall energy,

$$L\mathrm{\Delta }E(L)=L\left(E_\mathrm{g}(L)\frac{LE_\mathrm{g}(L+1)}{L+1}\right).$$
(8)
$`E_\mathrm{g}(L)`$ denotes the ground-state energy of the system with size $`L`$, which is supposed to be odd integer. For the system with odd number of plaquettes, one domain wall is created in the alternating alignment of vortex-antivortex structure. As would be apparent from the definition (8), $`\mathrm{\Delta }E`$ measures the extra domain-wall energy cost due to the kink. Hence, for ordered phase, the scaled kink energy should increase linearly with respect to $`L`$, while for disordered phase, $`L\mathrm{\Delta }E`$ vanishes exponentially; namely, $`\mathrm{e}^{L/\xi }`$ with correlation length $`\xi `$. At critical point, $`L\mathrm{\Delta }E`$ should be scale-invariant, because $`\mathrm{\Delta }E`$ is of scaling dimension $`1/L`$; note the relations $`\mathrm{\Delta }E1/\xi 1/L`$. In Fig. 2, we plotted $`L\mathrm{\Delta }E`$ for the same parameter range as that of Fig. 1. In fact, in the chiral phase $`1.8V_{}0`$, which is estimated from the above Binder-parameter analysis, we observe clear signature of the chiral-domain-wall energy cost. On the other hand, in the area $`0V_{}2.5`$, $`L\mathrm{\Delta }E`$ seems to be scale-invariant. In the former analysis of Fig. 1, we have concluded that in the region, the chirality is disordered. Hence, we conclude that in $`0V_{}2.5`$, the $`XY`$ order, rather than the chiral order, exhibits criticality. As are mentioned above, the strong-coupling regions of $`V_{}2.5`$ and $`V_{}1.8`$ are belonging to insulator phases with different characters. In terms of the pictures presented above, the behaviors of $`L\mathrm{\Delta }E`$ in these regions are readily understandable. The blowup of $`L\mathrm{\Delta }E`$ for $`V_{}2.5`$ is due to a kink (dislocation) formed in the island of particles; note that the particle occupation number is forced to be half-filled. On the other hand, the rapid closure of $`L\mathrm{\Delta }E`$ for $`V_{}1.8`$ is precisely due to the fact that the system is a Mott insulator.

Let us turn to the $`XY`$ order. In Fig. 3, we plotted the Binder parameter (6) for the in-plane spontaneous magnetization ($`XY`$ order), which corresponds to the superconductivity (gauge coherence) order parameter in the boson language,

$$M^2=M_{XY}^2=\underset{limj}{}(S_{li}^xS_{mj}^x+S_{li}^yS_{mj}^y),$$
(9)
for the same parameter range as that of former figures.

From the plot, we observe that the $`XY`$ order is kept critical for considerably wide range of parameter $`1.8V_{}2.5`$. We notice that this $`XY`$ critical phase does contain the chiral phase. This feature is quite contrastive to that of the semi-classical case $`S\mathrm{}`$, where the chirality is stronger (more stable) than the $`XY`$ order, and thus the chiral phase contains the $`XY`$ phase Granato93 . Moreover, it should be noted that in the critical region, the fixed-point value of the Binder parameter is kept hardly changed. In consequence, we found that the $`XY`$ correlation function is not influenced very much by $`V_{}`$, while the chiral sector is affected significantly by $`V_{}`$. These results show that those two sectors behave independently; in the commonly referred terminology, those two sectors are ‘separated.’ Such the situation where each mode is described by respective low-energy effective theory occurs commonly in one-dimensional physics. This point is discussed in Section 3.

Lastly, let us turn our attention to the possibility of spin-liquid state. The simplest way to realize the spin-liquid state is given just by setting $`t_{}`$ to a very large value. In that strong-interchain-coupling limit, the bonding and anti-bonding excitation branches are separated so that a band gap opens and thereby the ground state becomes disordered. Other than that rather trivial way, we found that the spin-liquid state is accessible from the set of parameters treated in the former figures 1-3 just through setting $`V_{}`$ repulsive; see Fig. 4, where we plotted the Binder parameter for the $`XY`$ order for $`t_{}=0.5`$, $`V_{}=1`$ and various $`V_{}`$.

We see that for attractive coupling $`0V_{}0.7`$, the $`XY`$ sector stays critical as is presented before. (The rapid suppression of $`U_{XY}`$ in $`V_{}0.7`$ is due to phase separation caused by strong attraction.) On the other hand, for repulsive coupling $`V_{}0`$, the $`XY`$ order becomes disordered eventually and driven to a spin-liquid phase. This behavior tells that the interchain coupling is a relevant perturbation in $`V_{}0`$. That feature coincides with that of the conventional non-frustrated ladder Schulz86 ; Strong92 ; Strong94 . We discuss this similarity in the next section.

## 3 Summary and discussions

We have diagonalized the Josephson-junction ladder subjected to the uniform magnetic field of half a flux quantum per plaquette (1). Unlike the semi-classical case, the chiral order suffers from strong disturbances due to quantum fluctuations. In order to stabilize the chiral order, we need to tune carefully the many-body interactions such as $`V_{}=0.6`$ and $`V_{}=1`$. It is surprising that the repulsive interchain interaction gives rise to the stabilization of the chiral order. This fact indicates that the particle exchange across the chains over the rungs, rather than the gauge coherence, is vital for pinning the vortices at each plaquette. On the contrary, we found that the $`XY`$ order (9) is insensitive to $`V_{}`$; the $`XY`$ order is kept critical for considerably wide range of parameters. This result indicates that the chiral and $`XY`$ sectors are separated.

In quantum spin systems, owing to the strong quantum fluctuations, the so-called spin-liquid state can appear Dagotto92 . The spin-liquid state, in fact, has been the main concern in the course of studies of the (non-frustrated) two-leg ladder Dagotto96 . Close to the chiral phase mentioned above, we found that a spin-liquid state emerges just through changing the sign of $`V_{}`$. Surprisingly, this behavior coincides with that of the conventional ladder; with the bosonization method Schulz86 ; Strong92 ; Strong94 , it was shown that the interchain coupling becomes relevant for $`V_{}0`$ (antiferromagnetic intrachain interaction). Hence, it is suggested that the $`XY`$ mode of our frustrated ladder and that of the conventional ladder behave similarly. This similarity might be reasonable, if we remember the separation of the chiral and $`XY`$ sectors mentioned above. Hence, it is suggested that the present critical phase, as well as that of conventional ladder Schulz86 ; Strong92 ; Strong94 , belongs to the universality class of the central charge $`c=1`$. Schulz86 ; Strong92 ; Strong94 . Direct identification of the universality class in terms of the finite-size conformal-field theory, for instance, may be exceedingly troublesome, because the excitation structure of the $`XY`$ sector is smeared by that of the chiral sector. This would remain for future study.

## Acknowledgments

The author is grateful to Professor I. Harada for helpful discussions.


format: text

text:
# Optical Spectrum of Main–, Inter– and Off-pulse Emission from Crab Pulsar

## 1 Introduction

The Crab pulsar (PSR J0534+2200 = PSR B0531$`+`$21) was the first pulsating radio source identified at optical wavelengths (Cocke, Disney, & Taylor 1969; Lynds, Maran, & Trumbo 1969). Most subsequent studies concentrated on the pulse shape in various photometric bands (Percival et al. 1993; Eikenberry et al. 1996; Perryman et al. 1999; Lundqvist et al. 1999). They generally suggested small color variations in the light curve. No spectroscopic data existed until Nasuti et al. (1996) obtained the first optical spectrum. They were able to measure the general flux distribution as well as hint to an unidentified intrinsic absorption feature at 5920 Å. Their spectrum is phase averaged, so no information on spectral variation with pulse phase can be extracted.

Here we present the first flux-calibrated optical spectrum, separately for the main-pulse, the inter-pulse and the off-pulse components. The spectra of both pulses are compared with phase-resolved photometric studies. Except for a very small ($`1`$% ) contribution from the pulsar (Golden et al. 2000), the off-pulse spectrum, also presented here, corresponds to the nebular emission along the line of sight of the pulsar.

## 2 Observations

Observations were obtained with the LFOSC spectrograph at the 2.12m telescope of the Observatorio Astrofísico Guillermo Haro in Cananea, Mexico. The useful wavelength range was $`\lambda \lambda 43709000`$ Å and the spectral resolution was $`15`$ Å. The slit of 2 arcsec width was kept in the N–S direction.

The light entered the spectroscope through a frequency and phase controlled chopper. The device, described by Čadež & Galičič (1996), was used before on the same telescope to study the optical light curve of the pulsar. A blade allows light to enter the spectrograph for 10% of each pulse period. This collection window is synchronized with the pulsar and kept on the desired pulse phase interval. We use the most recent of the Jodrell-Bank monthly Crab ephemerides, translated from the Solar-System Barycentre, including real-time corrections for the Earth orbital motion. Long photometric observations confirmed that the difference between the pulsar and the chopper phase remains within $`5\times 10^4`$ for at least 3 hours (Čadež et al. 2000). However the initial chopper–pulsar phase is set by maximizing the apparent magnitude of the pulsar, so that the absolute phase is only accurate to 0.03 (Galičič 1999).

Observations were obtained on Dec 13 and Dec 14 1998 and on Oct 18 and Oct 20 1999 under nearly photometric conditions and $`1\mathrm{"}`$ seeing. A total accumulated exposure on the main-pulse and the inter-pulse was 11700 sec and 9000 sec, while the total off-pulse observing time was 10200 sec. The off-pulse interval was chosen $`0.25`$ in phase after the main-pulse. Several direct images of the pulsar field were obtained in order to control the atmosphere transparency and phase of the pulsar. The wavelength calibration is generally good to 2 Å and somewhat deteriorating at the blue end. The flux calibration, obtained through several exposures of the star EG50 and Hiltner 600, is accurate to 5% in both color and absolute fluxes. All spectra were obtained between airmasses 1.01 and 1.3. The pulsar spectra were de-reddened assuming $`E_{BV}=A_V/3.1=0.51`$ (Percival et al. 1993) and

$$A_\lambda /E_{BV}=7.51\mathrm{log}(\lambda /\lambda _o)+1.15;\lambda _o=1\mu \mathrm{m}$$
(1)
which follows the average Galactic extinction curve (Savage & Mathis 1979) to better than 1% for $`0.48<\lambda /\lambda _o<1.0`$.

Spectra were reduced with standard IRAF routines used according to our particular purpose. First, all the spectra were aligned with respect to the tracing of the pulsar and of another star 11.4 arcsec South of the pulsar. In order to achieve an optimal flux calibration, a relatively wide aperture of 7 arcsec was used that included all pulsar-related flux and allowed for any wavelength variations of the focus of the spectrograph. All individual spectra were flux-calibrated but no nebular background was subtracted at this stage. The un-pulsed light was subtracted subsequently by subtracting the off-pulse flux-calibrated spectrum obtained just before or after. This technique makes an assumption-free subtraction of any nebular light, scattered light or telluric contributions. The chopper effectively enables us to take the pulsar off the sky and subtract the nebula emission at the exact pulsar’s position. This is important as the nebular emission is varying on arcsec scale (cf. Section 4).

## 3 Pulsar

The flux-calibrated de-reddened spectra of the main-pulse and inter-pulse are presented in the upper panel of Fig. 1. The surface flux of the underlying nebula is plotted in the lower panel. Both axes are plotted in logarithmic scale.

The colors of the main-pulse and inter-pulse continuum are identical, but the nebular continuum is redder (see Table 1).

The pulsar’s spectrum features no lines apart from telluric bands and weak nebular line residuals. Note that the $`[`$OI$`]`$ 5579 line is of telluric origin. The upper limits to equivalent widths of nebular lines are 5 and 10 Å for the main- and inter-pulse respectively. This translates to the very stringent limits on absolute fluxes of pulsed emission lines that are given in Table 2.

We did not find any trace of the 5920 Å absorption feature hinted by Nasuti et al. (1996). We conclude that they were right in attributing this feature to imperfect flux calibration of their spectrum.

Note that the nebular continuum is rather strong. The use of the 0.1-phase window with the chopper however suppressed it by a factor of 10. This emphasizes the advantage of our observing method.

## 4 Nebula

Davidson (1979) published the first study of selected condensations in the Crab nebula. Nasuti et al. (1996) discuss the nebular spectrum 2 arcsec East of the pulsar position. Our slit had a N–S orientation so a brief discussion of the spectrum along this direction as well as directly on the pulsar position is in order.

The shape of the (de-reddened) nebular continuum is constant along the slit ($`\alpha =0.4\pm 0.1`$) except for its strength. Spatial variations of surface flux continuum (de-reddened nebular continuum per square arc-sec) are shown in the lower panel of Fig. 2. The horizontal axis is the relative displacement in N-S direction with respect to the position of the pulsar and the vertical axis displays the continuum surface flux at 6000Å. There is a chance superposition of a star 11.4 arcsec South of the pulsar, so no flux is given around its position.

The upper panels of Fig. 2 give de-reddened surface fluxes of emission lines as a function of position relative to the pulsar. The lines are generally well-separated into approaching and receding components (Fig. 3). The exception is a blend of the red component of \[O III\] 4959 and the blue component of its 5007 counterpart which are plotted together. The H I 6563, \[N II\] 6548+6584 region was resolved assuming a standard 3 to 1 ratio of the \[N II\] line intensities.

All lines show well separated red and blue components, similar to those reported by Nasuti et al. (1996). Component velocities, given in Table 3, remain constant (to within 200 km/s) up to $`30`$ arcsec from the pulsar. All lines show similar velocity shifts, $`1800\mathrm{k}\mathrm{m}/\mathrm{s}`$, except for the very strong receding component of \[O III\] 4959 and 5007 lines, receding at double velocity.

The intensity of blue and red components is highly spatially variable and uncorrelated (Fig. 2). A simple interpretation is that of an expanding shell, proposed already by Davidson (1979). The receding components of forbidden nebular lines show a wide pronounced peak $`3`$ arcsec North of the pulsar. The slope of the condensation, however, extends well behind the pulsar itself. The approaching components are generally weaker and show a condensation $`4`$ arcsec South of the pulsar. The kinematics and strength of H$`\alpha `$ is similar to that of forbidden lines, but a pronounced receding component to the North of the pulsar is missing.

## 5 Discussion and conclusions

The main result of this work is the flux calibrated time resolved spectrum of the Crab pulsar. It shows no intrinsic spectral lines with the continuum following a power law with $`\alpha =+0.2\pm 0.1`$. The spectra of the main-pulse and the inter-pulse are identical: flux of the main-pulse multiplied by 0.32 equals that of the inter-pulse (Fig. 4), where the formal upper limit to the difference in the power law index for the main- and the inter-pulse is 0.01 ($`2\sigma `$).

The above value of $`\alpha `$ refers to de-reddened main- and inter-pulse spectra with nebular spectra subtracted at the exact area of pulsar’s PSF. This crucial subtraction was made possible by the dedicated chopper device that effectively enabled us to take the pulsar off the sky.

The underlying nebula is redder ($`\alpha =0.4\pm 0.1`$) and of varying intensity. We expect that this fact may explain why the power law index for the pulsar spectrum reported in the literature is slightly redder than deduced here. Nasuti et al. (1996) find $`\alpha =0.1\pm 0.01`$ in optical. Gull et al. (1998) give $`\alpha =0.1\pm 0.2`$ (if $`E_{BV}=0.51`$ is assumed) in the UV. Without the possibility of directly subtracting the off-pulse spectrum, one must make a model for the surrounding nebula. We simulated such modeling and took in an example the nebular spectrum 4 arcsec North and South from the pulsar as the background. This led to the derived value of $`\alpha =+0.1`$ in very good seeing (FWHM $`0.7`$ arcsec) or even 0.0 in more mediocre conditions. We also note that errors below 0.1 in $`\alpha `$ appear unrealistic, as they imply a color calibration accurate to better than 7% – a goal difficult to achieve on a patchy nebular background and given the uncertainties in the extinction law.

We detected no color variation in the pulsed light. Photometric investigation of Eikenberry et al. (1996) reports slight color variation on pulse slopes. Note, however, that our observations were obtained close to the centers of the main- and inter-pulse. So the two studies are somewhat difficult to compare: we present flux-calibrated spectra with slightly uncertain absolute phase, while Eikenberry et al. have an excellent phase definition but lack absolute color calibration. Romani et al. (1999) report on a relative decrease of the main pulse flux in the IR ($`1.4\mu `$$`<\lambda <1.8\mu `$m) with respect to the optical ($`0.38\mu `$$`<\lambda <0.83\mu `$m). We can not confirm a similar trend in our optical spectra, even if the braking of the power law beyond 800 nm could possibly be interpreted in this way. Since the nebular spectrum has a similar break, we would like another independent confirmation. A phase resolved spectrum of the Crab was presented by Perryman et al. (1999). In agreement with our results they detected no color dependence of the pulse shape in the optical. Note however that a limited resolution ($`100\mathrm{n}\mathrm{m}`$) and performance of the otherwise promising STJ junction detectors – together with modest seeing – prevented them from flux calibrating their spectrum and subtracting the nebular contributions.

By perfecting the stroboscopic technique we were able to obtain the spectrum of the pulsed emission of the Crab pulsar on a medium size telescope. We are now using the same technique to identify optical counterparts of known radio pulsars.

A.Č. and T.Z. acknowledge support from the Slovenian Ministry for Science and Technology. We would like to thank an anonymous referee whose comments improved the presentation of this paper.


text:
# On the Integrated Spectrum of the X-ray Binaries and the Origin of Soft X-ray Emission from the Bulge of M31

## 1 INTRODUCTION

M31, the Andromeda Galaxy, is the closest spiral galaxy and belongs to the same morphological type as the Milky Way. It provides us with a sample of X-ray sources at a uniform and relatively nearby range of $``$800 kpc (Stanek & Garnavich (1998)). Galactic absorption in the direction of M31 is rather low (Stark et al. (1992)), allowing the study of soft X-ray emission down to about 0.2 keV, which is not possible for the sources in the bulge or Galactic plane of our own Galaxy.

Observations with the Einstein and ROSAT satellites revealed multiple point X-ray sources within M31, and also significant unresolved X-ray emission in the bulge (Trinchieri & Fabbiano (1991); Primini, Forman & Jones (1993); Supper et al. (1997), hereafter Su97). The origin of the unresolved emission remains unknown. The possibilities include interstellar hot gas, or multiple faint point sources (stars or X-ray binaries) below the detection limit. Recently Irwin & Bregman (1999), hereafter IB99, suggested that whole X-ray emission of the bulge can be attributed to low-mass X-ray binaries (LMXBs), which are commonly assumed to be responsible for the hard (above several keV) X-ray emission of the X-ray faint early-type galaxies. Whether the same population of LMXBs might also be the origin for a significant very soft excess detected with Einstein and ROSAT in several X-ray faint early-type galaxies (e.g. Kim, Fabbiano & Trinchieri (1992); Pellegrini & Fabbiano (1994); Irwin, Sarazin & Bregman (2000)) remained, however, open and interesting question to explore.

In this letter we present results of our analysis of ROSAT PSPC data. We generated spatially resolved spectra of the bulge of M31 and for the first time extracted the spectrum of the unresolved X-ray emission. The results presented below demonstrate that the integrated spectrum of bright X-ray binaries is significantly harder than the spectrum of the unresolved emission in the bulge of M31.

## 2 OBSERVATIONS AND DATA REDUCTION

We analyzed the ROSAT PSPC observations discussed in detail by Su97. Descriptions of the instrument can be found in Aschenbach et al. (1981); Pfeffermann & Briel (1982), and Aschenbach (1988). We used all observations with a total exposure longer than 10,000 s that were taken before Oct 14, 1991111different response matrix must be used with data collected after this date, hence it is not possible to sum spectra together with previous observations. Most of our results were obtained from observation RP600068N00, when the central part of M31 was observed(Fig. 1). The angular resolution of ROSAT/PSPC does not allow a clean separation of point sources from the diffuse emission in the bulge of M31, because both are concentrated towards the center of the galaxy. However, we segregated these two components spatially as best as we could, collecting small regions around bright binaries into one spectrum, and photons from the areas without bright point-like sources into another. To identify point sources we used the ROSAT PSPC catalog published by Su97. To optimize the radius of photon collection for the point sources in the bulge, where diffuse emission is significant, we tried collection radii ranging from 30 to 60 arcsec for each source. As the radius increases we are collecting a great fraction of the source flux, but also the amount of diffuse contamination is increasing. We found that the choice of this radius does not change our conclusions described below. For X-ray binaries in the bulge, the typical radius was chosen to be 45 arcsec. Studying the integral spectrum of bright point sources outside of the bulge of M31, we collected all photons within 75 arcsec from each source.

We used only data within a central part of ROSAT’s field of view, that is inside the central ring of PSPC support grid with $``$ 20 arcmin radius. The background was measured in the same area of the detector, but outside of M31 bulge, in a region devoid of significant point sources and with much less unresolved emission than the bulge. The same background spectrum was used for all analyses of observation RP600068N00, both for bright point sources, and for unresolved emission, except when diffuse emission was subtracted from the point source spectrum as a background (as discussed below). We experimented with different background regions and found that our results do not depend on this choice. The spectral differences we discuss in this paper are much stronger than the expected background variations.

Various routines from the $`ftools`$ package (v.5.0) were used to generate spectra and other necessary files, which were then used for spectral fitting with $`xspec`$ (v.11.0) package. We ignored energy channels below 11 and above 219, and rebinned spectra so that each channel contained at least 25 counts.

## 3 SPECTRAL TEST RESULTS

Our choice of spectral model was based on recent spectral studies of M31 bulge reported by IB99 and Trinchieri et al. (1999), hereafter Tr99. Both studies revealed the presence of at least two spectral components in the integral spectrum of the bulge. The hard component was fitted with bremsstrahlung or power law model, while soft component - by a Raymond-Smith or MEKAL plasma emission model. We chose to use a power law for the hard component, mostly because the reported best-fit temperature of bremsstrahlung was well above the ROSAT upper energy limit. Unfortunately, ROSAT data cannot provide us with tight constraints for parameters of a two-component model. Table 1 presents several different combinations of parameters which all give reasonably good fits for the integrated spectrum of bulge emission, including all point sources and unresolved emission within inner 5 arcmin radius around the optical nucleus of M31 (see e.g. Crane, Dickel & Cowan (1992) for J2000 coordinates). The parameters of the fit are strongly correlated. There is the expected correlation between power law slope and absorption value, and also a correlation of the absorption with metallicity. The studies by IB99 and Tr99, which used ASCA and BeppoSAX data, also did not provide well-constrained fit parameters. The temperature parameter of soft component is 0.36–0.38 keV in IB99, but varies between 0.15 and 0.50 for the several fits presented by Tr99. Metallicity can be as low as 0.03 (IB99) or as high as 1 (fixed by Tr99). The power law slope from 1.2–1.3 (IB99) to 1.8–2.2 (Tr99), and absorption parameter $`N_H`$ from 4.4 $`\times `$ 1020 cm-2 (IB99) to 10–12 $`\times `$ 1020 cm-2 (Tr99). As one may see from Table 1, fraction of emission in the soft component depends strongly on fit parameters. Therefore, we need to freeze most of the parameters in order to track the relative contributions of the two components as a function of position.

For the most of our analysis we fixed the power law slope at 1.5, the absorption column at the Galactic value $`N_H`$ = 6.73 $`\times `$ 1020 cm-2 (Stark et al. (1992)), and the metallicity parameter at 0.03. This choice allowed us to get acceptable values of $`\chi ^2`$ value in all cases, and study relative significance of hard and soft components for various regions in the bulge of M31. It does not imply, though, that these values have a real physical sense. As illustrated by Table 2, one might, for example, choose to fix the metallicity at solar value and get different fit parameters. However, our main conclusion regarding the difference between the spectra of bright point sources and the unresolved emission is not model-dependent.

Using the model described above, we looked at the spectrum of the integral bulge emission (including all point sources) at various distances from the M31 center (less than 2 arcmin, between 2 and 5 arcmin, and between 5 and 8 arcmin). The contribution from the soft component decreased with distance from the center, in correlation with the decrease of the diffuse intensity in the image. The soft component fraction falls from 0.62 to 0.50 and then to 0.37, as we move from the inner 2 arcmin to the outermost annulus (Table 2).

For the next stage of our analysis, we separated bright point sources from unresolved emission. Overall approach was to draw small circular regions around each point source listed in Su97. Then we integrated these regions into one spectrum (‘XRBs’), and integrated the rest of the area into a second spectrum (‘diffuse’). Admittedly, we cannot separate point sources and unresolved emission near the core completely by this method, because the density of unresolved emission is maximal near the core of M31, where the concentration of point sources is also maximal, however, relative contribution of both kinds of sources is different in spectra of those two types. In Table 2 we present two spectra for each type, which differ by the area where the emission was collected. We tried to collect unresolved emission from the area where its density corresponds to average density of the same emission around bright sources. This allowed us to subtract ‘diffuse’ spectra from ‘XRBs’ spectra in an attempt to get ‘clean’ spectrum of point sources.

XRBs and diffuse emission near the core. For these spectra we extracted X-ray photons from the area near the center of M31, where the brightest point sources are concentrated and the density of diffuse emission is maximal. The separation of two types of progenitors is especially challenging in this area. 14 bright point sources from Su97 catalog, which are located within 6 x 8 arcmin ellipse around the center of M31 (see Inset b of Fig. 1), formed XRBs spectrum. Unresolved emission was collected over the ellipse 4 x 6 arcmin, but circles of 45 arcsec diameter around each catalogued source were excluded (Inset a of Fig. 1). For ‘XRBs near the core’ the contribution of soft component is large (0.48), but significantly lower than for the total emission from the same region (0.62). The soft component dominates the spectrum of unresolved emission collected in this area (0.8). If ‘diffuse near the core’ spectrum, normalized for area, is subtracted from ‘XRBs near the core’, then soft component is no longer significant in residual spectrum. In another analysis we fixed the metallicity parameter at 1. While, as expected, the absolute values of soft component fraction differ significantly for the different metallicity, soft component still contributes significantly more (0.47) to the spectrum of the unresolved emission than to the spectrum of ‘XRBs’ (0.22). The difference in the spectra illustrated by Fig. 2 does not depend on our choice of model parameters – it is easily seen in the distribution of counts over the spectral band.

XRBs and diffuse emission in the bulge. Here, we summed the spectra of point sources and collected an unresolved emission from a larger area, where they are more easily separable. The spectra of 14 bright sources that lie within 8 arcmin from the center, but outside of central 2 arcmin, were integrated into ‘XRBs in the bulge’(Inset c of Fig. 1). Unresolved emission was collected from inner 5 arcmin, less 45 arcsec radii around all point sources listed in Su97 (see Inset d of Fig. 1). Soft component is much stronger in the spectrum of unresolved emission (80% of total flux), than in the spectrum of ‘XRBs’ ($``$20% of total flux, not detected when spectrum of diffuse emission is subtracted).

We see that while all spectra can be satisfactorily fitted by the same model, the relative contribution of hard and soft component are significantly different for point sources and unresolved emission. Because we see the difference even with PSPC data which are seriously confused in this region, we expect that the difference would be even stronger if we obtained clean, well spatially resolved spectra. Furthermore, it motivates us to suggest that hard component in the bulge spectrum is generated by X-ray binaries of the bulge, while the soft component originates from a diffuse source of different nature. This suggestion can be checked with data coming from Chandra and XMM-Newton instruments.

Additionally we analyzed the integrated spectrum of other X-ray binaries in M31. In this case we added together spectra of all bright point objects located outside of central 8 arcminutes. We did not add spectra of the sources marked as foreground stars or background galaxies in Su97. Altogether 61 bright point sources were added together to the spectrum labeled as ‘XRBs in the disk’ in Table 2. The integrated spectrum can be fit by single power-law component with a slope around 1.4 and absorption slightly higher than the Galactic value (Fig. 2). No evidence for soft component was found. This spectrum is very close to the spectrum of ‘XRBs in the bulge’ with diffuse emission subtracted. We need to note that all spectra labeled as ‘XRBs’ in Table 2 include all bright compact sources in M31 and does not exclude supernova remnants or globular clusters. We suggest that a single power-law spectrum in ROSAT energy band is a more suitable template for the population of bright compact X-ray emitters, than the spectrum of M31 bulge, which includes a significant soft component.

## 4 DISCUSSION

Our spatially resolved spectral analysis revealed a qualitative difference between the unresolved emission and the bright point sources in the bulge of M31. We have performed several tests and that gave consistent results. We found that soft component is much more prominent in the spectrum of unresolved emission than in the integral spectrum of bright point sources (presumably LMXBs).

Previous studies were based on the spectrum of integral X-ray emission from the bulge, in comparison with spectra for individual sources. IB99 argued that LMXBs may be responsible for both hard and soft X-ray components in the spectrum of the M31 bulge. They suggested furthermore that unresolved emission in the bulge is composed of the same LMXBs and has the same energy spectrum. Our analysis shows conclusively that this is not the case. The spectrum of the unresolved emission is significantly different and hence it must be dominated by the sources of a different nature. These sources are responsible for most of soft excess observed from M31. We cannot rule out the possibility that LMXBs generate some of the soft component, but the contribution of this component to their spectra is lower than in the unresolved emission. The nature of the unresolved emission remains unclear, however, we believe that the observation of an extended radio source of a similar size at the same position in M31 (Hjellming & Smarr (1982)) supports an interstellar gas hypothesis. The reality of diffuse emission in M31 is, in a sense, not a surprise - in our own Galaxy, we have long known of diffuse emission from the Galactic plane down to 3 keV (Valinia & Marshall 1998, Kaneda et al. 1999, and references therein). Although our measurements of M31 did not extend to these energies, and the soft component that we detected would be hidden by absorption in our own Galaxy, the Galactic measurements give evidence for thermal emission from a perhaps hotter phase of the interstellar medium.

Our study suggests that the resolved bright sources in the bulge of M31 may not be responsible for any of the soft X-ray excess in the spectrum. In this case, their integral spectrum strongly resembles the bright X-ray sources outside of the M31 bulge. This looks surprising because it is commonly accepted that disk and bulge sources belong to different populations of X-ray binaries. It may appear, however, that LMXBs cannot be distinguished from high-mass X-ray binaries solely on the basis of their X-ray spectra in ROSAT band. After all, the compact objects in both types of binaries are of the same nature.

Our conclusion that the soft spectral component mainly emerges from the unresolved emission does not depend on a quantitative deconvolution of the bulge spectrum. This is important, because the exact fraction in the soft component depends strongly on the spectral model parameters, and cannot be well evaluated with existing data. Our poor state of knowledge is reflected by the difference between Tr99 ( 15% soft fraction) and IB99 (38% or more). The fraction of the total bulge emission from point sources is also unclear. Trinchieri and Fabbiano (1991) found that 75% of the bulge emission was resolved into 46 point sources, but Primini et al. (1003) estimated that 45 point sources accounted for only 58% of the emission from the inner 5’, with the rest unresolved. With these uncertainties, it is difficult to determine the fraction of soft emission that comes from LMXBs based on integrated spectra, motivating our spatially-resolved analysis. New data from XMM-Newton and Chandra promises to clarify the picture.

The interpretation of our results goes far beyond the case of M31, because this galaxy, due to its proximity, often serves as a testbed for other early-type galaxies (see e.g. Irwin & Sarazin, 1998). The X-ray spectra of many X-ray-faint galaxies had been found to be the sum of a hard component and a very soft component, similar to the spectrum of the M31 bulge. Our analysis suggests that this soft component comes mostly from a hot interstellar gas, or some unknown population of faint, soft X-ray sources, but not from LMXBs.

## 5 ACKNOWLEDGMENTS

The ROSAT public data were extracted from the HEASARC electronic archive operated by the Goddard Space Flight Center (NASA). The authors are thankful to an anonymous referee and to J.Irwin for useful comments. KB is glad to acknowledge helpful discussion of the draft by M.Revnivtsev and A.Vikhlinin.


text:
# The Harmonic and Sideband Structure of the Kilohertz Quasi-Periodic Oscillations in Sco X-1

## 1 Introduction

It is four years now since the kilohertz quasi-periodic oscillations (kHz QPOs) were discovered in the persistent flux of Scorpius X-1 \[1996a\] and 4U 1728–34 . In the meantime, similar kHz QPOs have been seen in some 20 other low-mass X-ray binaries (LMXBs; see van der Klis 2000 for a review). These QPOs often appear in pairs, with frequencies $`\nu _1`$ and $`\nu _2`$ ($`\nu _2>\nu _1`$) between $`400`$ and $`1300`$ Hz, which in a given source can shift by a few hundred Hz, apparently as a function of mass accretion rate.

Most of the models proposed so far assume that one of the kHz QPOs reflects the Keplerian orbital motion at some preferred radius in the accretion disc (e.g., Miller, Lamb & Psaltis 1998; Stella & Vietri 1999; Osherovich & Titarchuk 1999), but there are other explanations as well (Klein et al. 1996a,b; Jernigan, Klein & Arons 2000). In recent discussions (Lamb & Miller 1999; Stella 1999; Psaltis 1999) it was emphasized that empirical discrimination between two of the leading classes of models is possible in principle by studying the harmonic and sideband structure of the kHz QPOs. Here we concentrate only on these two model classes

In the ‘sonic-point’ model (SPM; Miller et al. 1998), the QPO at $`\nu _2`$ (the upper QPO) is produced at the radius where the radial flow velocity in the disc turns from subsonic to supersonic (the sonic radius), and the QPO at $`\nu _1`$ (the lower QPO) originates by a beat between the upper QPO and the spin frequency of the neutron star. In the ‘relativistic-precession’ model (RPM; Stella & Vietri 1999) the QPO at $`\nu _2`$ is also assumed to be Keplerian, but the QPO at $`\nu _1`$ is produced by the apsidal precession of a slightly non-circular inner accretion disc. The QPO frequencies in the RPM are calculated for test particles in purely geodesic relativistic motion, i.e., neglecting the hydrodynamical and radiative effects of the accretion flow. However, Psaltis & Norman  have recently proposed a dynamical model in which the QPOs are produced by oscillations in the accretion disk. In this model, which we will call ‘transition-radius’ model (TRM), there is a transition radius in the accretion disc that acts as a band-pass filter with resonances near the orbital and periastron-precession frequencies.

Besides the main peaks at $`\nu _1`$ and $`\nu _2`$, the SPM and the TRM predict other (weaker) harmonics and sidebands of these QPOs, at specific frequencies. For instance, the SPM predicts a relatively strong harmonic of the lower QPO at $`2\nu _1`$ (see Table 3 of Miller et al. 1998 for a list of other sideband peaks predicted by the SPM), whereas the TRM predicts a sideband at $`2\nu _2\nu _1`$ (cf. eq.  in Psaltis & Norman 2000). In principle, the detection of a QPO at $`2\nu _1`$ and a non-detection of a QPO at $`2\nu _2\nu _1`$ would tend to rule out the TRM, whereas the detection of a QPO at $`2\nu _2\nu _1`$ and a non-detection of a QPO at $`2\nu _1`$ would tend to rule out the SPM .

In this Letter, we use data from the Rossi X-ray Timing Explorer (RXTE) to search for the predicted harmonics and sidebands in the power spectrum of the kHz QPO source Sco X-1. Because of its very high flux, Sco X-1 has extremely significant kHz QPOs, and hence a very sensitive study of any harmonic structure is possible. We do not detect any of these secondary QPOs. The fact that we detect the kHz QPOs at $`\nu _1`$ and $`\nu _2`$, but none of these other peaks sets severe constraints on the models currently proposed to explain the kHz QPOs in LMXBs.

## 2 Observations

We used data from the Proportional Counter Array (PCA; Jahoda et al. 1996) on board RXTE (Bradt, Rothschild & Swank 1993) taken on 1996 February 14, 18, 19, 1996 May 24 – 28, 1997 March 15, 1997 April 18 – 24, 1997 August 22, 1998 January 2 – 8, 1998 February 27, 28, 1998 May 30, 31, 1998 June 1, 2, 1998 July 2 – 5, 1999 January 6, 8 – 11 and 13 – 16. Each observation (i.e., each part of the data with a unique RXTE ID number) consists of data blocks of $`60`$ s to $`3,700`$ s interrupted by passages of the satellite through the South Atlantic Anomaly and occultations of the source by the Earth. The total usable time was $`630`$ ks.

To avoid detector safety triggers, telemetry saturation, and to reduce the dead-time effects produced by the high count rate of Sco X-1, some observations were carried out with the source slightly off-axis, with some of the five proportional counter units of the PCA switched off, recording only photons detected by the upper anode chain of the PCA, recording only photons from a limited energy range, or using a combination of these constraints. In all cases, high-time resolution data were available with a time resolution of 0.25 ms or better. For our analysis below we combined all the available data, irrespective of whether they were collected using any of the above observational constraints. When single and double-event data were recorded in parallel (see van der Klis et al. 1996b), we combined them off-line to enhance the sensitivity.

## 3 Analysis and Results

We divided the high time resolution data into 16 s segments, and produced a power spectrum for each of these segments up to a Nyquist frequency of 2048 Hz (for $`50`$ % of the data we also produced power spectra up to a Nyquist frequency of 4096 Hz). For 69 % of the observations the power spectra were calculated using the full PCA energy band. For the rest of the power spectra we used data from selected energy bands (24 % of the power spectra were calculated between $`518`$ keV, 5 % between $`560`$ keV, and 2 % between $`218`$ keV) because we found that the kHz QPOs were more significant in those energy bands, or because those were the only energy bands available. Finally, we averaged together groups of 8 contiguous 16 s power spectra to produce average power spectra.

We searched these average power spectra for kHz QPOs, at frequencies $`\mathrm{¿}\mathrm{}250`$ Hz. We did this by first identifying those power spectra that showed a strong QPO, and visually estimating its frequency, $`\nu _0`$. When two QPOs were present in the power spectra, we always picked the one at higher frequency. It turned out that in the cases where only one kHz QPO was visible, it was always the upper kHz QPO (see below). We then fitted the power spectra, in the range $`\nu _0100`$ Hz to $`\nu _0+100`$ Hz, using a function consisting of a constant, a power law and one Lorentzian. We discarded the power spectra for which the QPOs were less than $`3\sigma `$ significant. We note that in this manner we might have discarded data with weak QPOs that could have been detected averaging together more data. There were 1,384 average power spectra with significant kHz QPOs, equivalent to 177,152 s of data.

For those observations where we detected only one QPO in the average power spectra we applied the ‘shift-and-add’ technique \[1998a\] to try and detect the second QPO: Based on the frequency of the detected QPO, we aligned and averaged all the power spectra of a single observation. In those cases, this procedure revealed a second kHz QPO, which in all cases was at a lower frequency than the QPO peak that we used to align the power spectra, showing that all the initial frequency measurements in the 1,384 average power spectra corresponded to the upper kHz QPO.

Next, we measured the frequency separation between the kHz QPOs as a function of the upper QPO frequency, $`\nu _2`$. We proceeded as follows: We aligned the 1,384 power spectra using the upper kHz QPO as a reference; we grouped the data in 49 sets, each of them containing $`40`$ to $`210`$ power spectra, such that $`\nu _2`$ did not vary by more than $`510`$ Hz within each set, and we combined these aligned spectra to produce an average power spectrum for each set. We fitted these 49 power spectra in the range $`4001300`$ Hz using a function consisting of a constant, a power law and two Lorentzians. The fits were good, with reduced $`\chi ^21.1`$ for 279 degrees of freedom, and the significance of both peaks was always $`>3\sigma `$.

Fig. 1 shows that $`\mathrm{\Delta }\nu =\nu _2\nu _1`$, the frequency difference between the two QPOs, decreases from $`310`$ Hz to $`240`$ Hz, as $`\nu _2`$ increases from $`840`$ Hz to $`1100`$ Hz. This figure can be compared to fig. 3a of van der Klis et al. , who first reported the decrease in $`\mathrm{\Delta }\nu `$ with $`\nu _2`$ in Sco X-1. Because we included more data, and because we used the shift-and-add technique to measure $`\mathrm{\Delta }\nu `$, the errors are smaller in the present figure than in the previous one, and some structure, particularly between 900 and 1020 Hz, becomes apparent (perhaps this structure is related to the ‘bump’ seen in the plot of $`\mathrm{\Delta }\nu `$ vs. $`\nu _1`$ of 4U 1608–52 at $`\nu _1700`$ Hz; see fig. 3 in Méndez et al. 1998b).

From Fig. 1 we can read off $`\nu _1`$ as a function of $`\nu _2`$; because we already know $`\nu _2`$ for each average power spectrum, we can calculate the expected frequencies of hypothetical signals at multiples of $`\nu _1`$ and $`\nu _2`$, or at frequencies that are combinations of these two frequencies, for each power spectrum. As we described in §1, some of these frequencies are important in the context of the models proposed to explain the kHz QPOs.

We calculated the following quantities: $`2\nu _1`$, $`2\nu _2`$, $`\nu _2\nu _1`$, $`\nu _2+\nu _1`$, $`2\nu _1\nu _2`$, $`2\nu _2\nu _1`$, $`2\nu _1+\nu _2`$, $`2\nu _2+\nu _1`$ and $`2(\nu _2\nu _1)`$. We then aligned, in turn, the power spectra on each of the above frequencies, and averaged them together to try and detect a signal present at any of these frequencies. We fitted these shifted and averaged power spectra in a segment of $`400`$ Hz that included the expected frequency, using a constant plus a Lorentzian with the centroid fixed in each case at the expected frequency, and with a fixed FWHM of 200 Hz, approximately equal to the sum of the FWHM of each of the QPOs \[1996b\]. This value is about two times larger than the largest FWHM measured in the QPOs of this and other sources, and yields conservative upper limits to the unobserved harmonics and sidebands; tighter upper limits could be obtained by fixing the FWHM to a smaller value. No detections result from choosing a smaller FWHM in the current data set. In a few cases we added a power law to the fit to take account of a sloping continuum. The shape of the power spectrum at high frequencies is dominated by dead-time effects  which, in the case of Sco X-1 at a count rate that exceeds $`25,000`$ c s-1 PCU-1, are large and not yet sufficiently well understood to predict the shape of the high-frequency part of the power spectrum accurately. We note that a better knowledge of this effect could yield tighter upper limits than those that we report in this paper.

We did not detect a significant QPO at any of the above frequencies. As an example, in Fig. 2 we show two power spectra that were shifted to $`2\nu _1`$ and to $`2\nu _2\nu _1`$, the frequencies of the main high-frequency peaks (besides the kHz QPOs) predicted by the models of Miller et al.  and Psaltis & Norman , respectively.

Upper limits are shown in Table 1. To calculate these upper limits we varied the amplitude of the Lorentzian until the $`\chi ^2`$ of the fit increased by 2.71 with respect to the best-fitting value (95 % confidence level for a single parameter); the upper limits we quote in Table 1 represent the change in the amplitude for this change in the $`\chi ^2`$. Because of the observational constraints described in §2, and because we combined power spectra computed from data in different energy bands, we cannot give upper limits in units of fractional rms amplitude ; instead, all the upper limits in Table 1 are in units of the observed fractional rms amplitude of the QPO at $`\nu _2`$. Any peaks at the predicted frequencies are between 11 and 156 times weaker than the peak at $`\nu _2`$, which translates into rms amplitude ratios between 0.08 and 0.30.

## 4 Discussion

We have measured the frequencies of the two simultaneous kHz QPOs, $`\nu _1`$ and $`\nu _2`$, in Sco X-1 in a large data set. We have used these measurements to calculate the expected frequencies of hypothetical signals at harmonics and sidebands of these QPOs, as predicted by two classes of kHz QPO models. We used a sensitive technique to search for these harmonics and sidebands in the power spectra of Sco X-1, but none was detected. Their power is between 10 and 100 times less than the power of the upper kHz QPO (95 % confidence).

Oscillations produced close to the surface of the neutron star can be attenuated in a scattering corona. The attenuation is generally larger for oscillations produced by a pencil beam sweeping the surroundings (beaming oscillations), than for oscillations due to actual changes of the luminosity with time (luminosity oscillations; Brainerd & Lamb 1987; Kylafis & Klimis 1987; Kylafis & Phinney 1989; Miller et al. 1988). The attenuation factor (the ratio of the amplitude $`A_{\mathrm{}}`$ of the oscillations at infinity to their original amplitude $`A_0`$) as a function of frequency $`\nu `$, for a scattering corona of a radius $`R`$ and optical depth $`\tau `$ is given by:

$$\frac{A_{\mathrm{},\mathrm{lum}}}{A_{0,\mathrm{lum}}}=2^{3/2}xe^x+e^\tau $$
(1)
for a luminosity oscillation, and by

$$\frac{A_{\mathrm{},\mathrm{beam}}}{A_{0,\mathrm{beam}}}=\frac{2^{5/2}xe^x}{1+\tau }+e^\tau ,$$
(2)
for a beaming oscillation, with $`x=(3\pi \nu R\tau /c)^{1/2}`$ .

We can use these relations to calculate the largest unattenuated amplitudes allowed by the data in Table 1. Because at each frequency the upper limit on the ratio of observed amplitudes is $`A_{\mathrm{},\mathrm{y}}(\nu )/A_{\mathrm{},\mathrm{beam}}(\nu _2)`$ (Table 1), the upper limit on the ratio of unattenuated amplitudes $`A_{0,\mathrm{y}}(\nu )/A_{0,\mathrm{beam}}(\nu _2)`$, can be obtained by minimizing

$$\alpha _\mathrm{y}(\nu )=\frac{A_{\mathrm{},\mathrm{y}}(\nu )/A_{0,\mathrm{y}}(\nu )}{A_{\mathrm{},\mathrm{beam}}(\nu _2)/A_{0,\mathrm{beam}}(\nu _2)}$$
(3)
with respect to $`\tau `$ and $`R`$, with the additional constraint that $`A_{\mathrm{},\mathrm{beam}}(\nu _2)/A_{0,\mathrm{beam}}(\nu _2)>0.01`$ (because the QPO at $`\nu _2`$ is detected at an amplitude $`\mathrm{¿}\mathrm{}1`$ %; van der Klis et al. 1997). In eq. , $`y`$ stands either for luminosity or beaming (see eq.  and ), and we have assumed that the QPO at $`\nu _2`$ is a beaming oscillation (Miller et al. 1998).

The upper limits to the unattenuated amplitudes of the sidebands and harmonics of the kHz QPOs, both for luminosity and beaming oscillations, are shown in Table 2. In units of the amplitude of the QPO at $`\nu _2`$, the amplitudes of these harmonics and sidebands range from 0.11 to 0.54 for a luminosity oscillation, and from 0.21 to 2.16 for a beaming oscillation (assuming beaming for the QPO at $`\nu _2`$). From the minimization of $`\alpha _\mathrm{y}(\nu )`$ at the frequencies listed in Table 2, we get that optical depths $`\tau 110`$ and radii $`R1350`$ km for the putative scattering cloud in Sco X-1 would be most effective at suppressing the sidebands and harmonics. These values are comparable to those obtained by by Vaughan et al. , $`\tau \mathrm{¡}\mathrm{}5`$ and $`R\mathrm{¡}\mathrm{}80`$ km, for 4U 1608–52 from measurements of the time lags of the kHz QPOs (see also Kaaret et al. 1999). We want to mention that the upper limits given in Table 2 have been calculated assuming properties of the scattering corona that maximize the attenuation at each frequency independently. In reality, if the corona exists it would have fixed values of $`\tau `$ and $`R`$, so that if it is effective at suppressing the oscillations at some frequencies, it will be less effective at others.

In conclusion, there is a remarkable lack of harmonic and sideband structure in the kHz QPOs of Sco X-1. None of the frequencies predicted by the two main classes of models are detected, with upper limits implying they are one to two orders of magnitude weaker in power than the main peaks. This sets severe constraints on either the purity of the kHz QPO generation mechanism in Sco X-1, or the parameters of the scattering corona that surrounds the neutron star.

## Acknowledgments

We would like to thank various participants of the 1999 workshop on X-ray Probes of Relativistic Astrophysics at the Aspen Center for Physics for pleasant and extremely fruitful discussions. We are specially grateful to Coleman Miller, Fred Lamb, Luigi Stella, and Dimitrios Psaltis. We thank Coleman Miller, Fred Lamb and Tomaso Belloni for comments that helped us to improve the original manuscript. MM acknowledges comments from Mikhail Revnivtsev concerning the dead-time effects of the PCA. This work was supported by the Netherlands Research School for Astronomy (NOVA), the Netherlands Organization for Scientific Research (NWO) under contract number 614-51-002 and the NWO Spinoza grant 08-0 to E. P. J. van den Heuvel. MM is a fellow of the Consejo Nacional de Investigaciones Científicas y Técnicas de la República Argentina. This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center.


text:
# Using top-c for Commodity Parallel Computing in Cosmic Ray Physics Simulations

## 1 Introduction

Ultra high energy cosmic rays are observed indirectly through detection of the extensive air showers that are produced when they travel through the atmosphere. To adequately interpret the measured observables and to be able to infer the properties of the incident primary particle, a full Monte Carlo treatment of the extensive air shower is neeeded. The CPU time required rises with the primary energy. For example, for primary energies around $`10^{20}`$ eV a shower contains about $`10^{11}`$ secondary particles. The amount of computing time required to follow all the particles seems to be prohibitive. Traditionally, sampling techniques are used to reduce the number of particles tracked .

In this article we describe an ongoing program to use commodity parallel computing for fast Monte Carlo simulations . The aim is to go beyond the simple event-level parallelism which is commonly used today and actually run individual events faster than would be possible on a single workstation or PC.

## 2  geant4

For a variety of reasons, in no small part driven by the wish to work with software which is likely to see use in the future, we decided to try to parallelize geant4 , the C++ rewrite of the older (fortran77) geant3. geant4 is an object-oriented simulation package that provides general-purpose tools for defining and simulating detector geometry, material properties, particle transport and interactions, visualization, and all relevant physics processes. Its versatility allows it to be employed in applications beyond its traditional usage in High Energy Physics experiments, from the medical and biological sciences to Cosmic Ray Physics .

## 3 top-c

top-c (Task Oriented Parallel C)  was initially designed with two goals in mind:

1. to provide a framework for easily developing parallel applications;
2. to build in the ability to tolerate the high latency typically found on Beowulf clusters. 111The term “Beowulf cluster” refers to a cluster of systems running Linux and connected by ethernet.

The package is freely available . The same application source code has been run under shared and distributed memory (SMP, IBM SP-2, NoW, Beowulf cluster). A sequential top-c library is also provided to ease debugging. The largest test to date was a computer construction of Janko’s group over three months using approximately 100 nodes of an IBM SP-2 parallel computer at Cornell University .

The top-c programmer’s model  is a master-slave architecture based on three key concepts:

1. tasks in the context of a master/slave architecture;
2. global shared data with lazy updates; and
3. actions to be taken after each task.

Task descriptions (task inputs) are generated on the master, and assigned to a slave. The slave executes the task and returns the result to the master. The master may update shared data on all processes. Such global updates take place on each slave after the slave completes its current task. The programmer’s model for top-c is graphically described below.

## 4  Parallelization of geant4 Using top-c

The task-oriented approach of top-c is ideally suited to parallelizing legacy applications. The tactic in parallelizing geant4 was to perturb the existing software as little as possible and to modify just the section of the code which handles particle tracking and interaction (a frequent operation) to allow it to run on multiple CPU’s. The largest difficulty was in marshalling and unmarshalling the C++ geant4 track objects that had to be passed to the slave processes. Marshalling is the process by which one produces a representation of an object in a contiguous buffer suitable for transfer over a network, and unmarshalling is the inverse process.

We developed a 6-step software methodology to incrementally parallelize geant4, allowing us to isolate individual issues. The six steps were:

1. the use of .icc (include) files to isolate the code from the original geant4 code;
2. collecting the code of the inner loop in a separate routine, DoTask(), whose input was a primary particle track, and whose output was the primary and its secondary particles;
3. marshalling and unmarshalling the C++ objects for particle tracks
4. integrating the marshalled versions of the particle tracks with the calls to DoTask();
5. adding calls to top-c routines such as TOPC\_init(), TOPC\_submit\_task\_input and then testing as the marshalled particle tracks were sent across the network;
6. and finally adding CheckTaskResult(), which inspected the task output, and added the secondary tracks to the stack, for later processing by other slave processes.

Prior to the fifth step, all debugging was in a sequential setting. The maturity of the top-c library then allowed us to create fully functioning parallel code in less than a day.

## 5 Discussion

geant4 (approximately 100,000 lines of C++ code) was successfully parallelized using top-c. In the future we plan to perform timing tests on a long run using many processors. Initial results for the example described indicate that a single task in our application requires approximately 1 ms of CPU time. Hence, it will be essential to submit approximately 100 particles for a single slave process to compute in order to overcome network overhead. Optimization of the parallel implementaion is underway.

top-c seems to be well-suited to the problem of parallelizing geant4, and would likely be well-suited to other high energy physics and cosmic ray applications as well. Its flexibility and simplicity makes it possible to envision enormous speedups for geant4 within a single event, something not often considered in high energy experiments, but offering advantages over the usual event-by-event parallelism, especially during interactive data analysis and code or hardware design.

Of particular interest is the parallelization of existing cosmic ray simulation programs such as aires  and corsika . Although written in fortran, such programs are in fact often converted to C for compilation using f2c, and can certainly be linked with other C programs, so we anticipate no major obstacles. We are always interested in collaboration with other groups who may have needs for the speedups that our methodology offers.


text:
# Nonuniqueness of the Potentials of Spin-Density-Functional Theory

## Abstract

It is shown that, contrary to widely held beliefs, the potentials of spin-density-functional theory (SDFT) are not unique functionals of the spin densities. Explicit examples of distinct sets of potentials with the same ground-state densities are constructed, and general arguments that uniqueness should not occur in SDFT and other generalized density-functional theories are given. As a consequence, various types of applications of SDFT require significant corrections or modifications.

The Hohenberg-Kohn (HK) theorem , which guarantees that the ground-state single-particle density on its own is sufficient to uniquely determine all observables of a many-body system, is one of the most remarkable theorems of quantum mechanics. It is also at the heart of density-functional theory (DFT), one of the most popular many-body methods . In the present paper a fundamental problem in the extension of this theorem to spin-polarized situations is analysed, and consequences for practical applications of SDFT are pointed out.

In the case of the original formulation of DFT, in which the basic variable is the particle-density $`n(𝐫)`$, the HK theorem can be cast in the form of two logically independent one-to-one maps . Quantum mechanics guarantees that, for given particle-particle interaction $`\widehat{U}`$ and particle number $`N`$, the potential in which the particles move determines the ground-state many-body wave function (via solution of Schrödinger’s equation), which in turn determines the ground-state density (by simple integration). The essence of the original Hohenberg-Kohn theorem for nondegenerate ground states is that both of these maps are invertible: the ground-state density $`n(𝐫)`$ uniquely determines the ground-state wave function $`\mathrm{\Psi }(𝐫_1,\mathrm{}𝐫_N)`$, which in turn determines, up to an additive constant, the potential ,

$$v(𝐫)\stackrel{(1)}{}\mathrm{\Psi }(𝐫_1,\mathrm{}𝐫_N)\stackrel{(2)}{}n(𝐫).$$
(1)
Map two implies that every ground-state observable is a unique functional of the density, whereas the combined map allows the much stronger statement that the potential (and, given the interaction, thus every observable) is a unique density functional. These abstract maps have found extremely important practical applications in the form of the Kohn-Sham (KS) formulation of DFT , which is used for almost all band-structure calculations in solid-state physics, and a rapidly increasing number of electronic-structure calculations in quantum chemistry. Many of these applications, however, do not employ the original formulation of DFT, but spin-density-functional theory , in which the fundamental variables are the spin-resolved particle densities $`n_\sigma (𝐫)`$. In SDFT the map from spin densities to wave functions is easily established, but that from wave functions to potentials could, in spite of considerable effort  not be proven and remains an, albeit popular, conjecture. In the early days of SDFT von Barth and Hedin  already pointed out that the uniqueness of the spin-dependent potentials is not guaranteed, and explicitly constructed two potentials which, when used in a one-body Hamiltonian, have common eigenstates. The question whether a similar construction is possible in the many-body case, however, remained open, and even in the one-body case it was objected that these common eigenstates were not necessarily common ground states .

In this paper the general problem of uniqueness of the SDFT potentials is solved by (i) giving a general characterization (not limited to SDFT) of possible sources of nonuniqueness in DFT, (ii) explicitly constructing different SDFT potentials with a common ground state, and (iii) pointing out the physical reason behind nonuniqueness.

To address point (i), let $`v^{}=v+\mathrm{\Delta }v`$ and $`𝐁^{}=𝐁+\mathrm{\Delta }𝐁`$ be electrostatic and magnetic fields that are supposed to give the same (many-body) ground state $`\mathrm{\Psi }`$ as $`v`$ and $`𝐁`$. $`\mathrm{\Delta }v`$ and $`\mathrm{\Delta }𝐁`$ then necessarily satisfy

$$d^3r\left[\widehat{n}(𝐫)\mathrm{\Delta }v(𝐫)\widehat{𝐦}(𝐫)\mathrm{\Delta }𝐁(𝐫)\right]\mathrm{\Psi }=\mathrm{\Delta }E\mathrm{\Psi }.$$
(2)
A necessary condition for nonuniqueness is thus that, given $`\mathrm{\Psi }`$, one can find a linear combination of the density operators ($`\widehat{n}`$ and $`\widehat{m}`$ in the case of SDFT) of which $`\mathrm{\Psi }`$ is an eigenfunction. This is automatically the case if a linear combination exists that is a constant of motion, since constants of motion commute with the Hamiltonian and thus have the same eigenfunctions (assumed nondegenerate for convenience). If the Hamiltonian expressed in terms of $`v`$ and $`𝐁`$ has a gap between its ground state and its first excited state, one can always make sure that $`\mathrm{\Psi }`$ remains ground state of the new Hamiltonian, containing $`v^{}`$ and $`𝐁^{}`$, by making the coefficients of the linear combination, $`\mathrm{\Delta }v`$ and $`\mathrm{\Delta }𝐁`$, sufficiently small.

From any extensive constant of motion (i.e., one which is a linear combination of the density operators) one can thus systematically construct families of potentials with the same ground state. As an example consider a system with an energy gap, for which the total magnetization $`\widehat{M}_z=d^3r\widehat{m}_z(𝐫)`$ is a constant of motion. The choice $`\mathrm{\Delta }v=0`$, $`\mathrm{\Delta }𝐁=\overline{B}𝐮_z=const`$, where $`𝐮_z`$ is a unit vector in the z-direction, then is clearly a suitable linear combination, as long as $`\overline{B}`$ is not large enough to induce level crossings . The presence of a gap is crucially important, since it guarantees that for sufficiently small values of $`\overline{B}`$ such level crossings do not take place, so that $`\mathrm{\Psi }`$ remains the ground state.

It is important to emphasize at this point that the exchange-correlation (xc) functionals are still guaranteed to be unique density functionals by the second part of the HK theorem. Since the Hartree potential, too, is manifestly a unique density functional, the nonuniqueness can only arise from the external potential contribution to the effective potentials. Practical consequences of this are pointed out towards the end of this paper.

Examples for nonuniqueness obtained from constants of motion are below referred to as systematic nonuniqueness. Apart from these, Eq. (2) may, under suitable circumstances, also permit solutions not systematically associated with conserved quantities. To construct an example for such accidental nonuniqueness first recall how SDFT describes a spin-polarized system. After choosing the xc functional and specifying the corresponding spin-dependent Kohn-Sham potential $`v_{s,\sigma }(𝐫)`$, one solves selfconsistently the spin-dependent Kohn-Sham equation,

$$\left[\frac{\mathrm{}^2}{2m}^2+v_{s,\sigma }(𝐫)\right]\varphi _{n\sigma }(𝐫)=ϵ_{n\sigma }\varphi _{n\sigma }(𝐫),$$
(3)
and calculates the spin densities from

$$n_\sigma (𝐫)=\underset{n}{\overset{N_\sigma }{}}\varphi _{n\sigma }^{}(𝐫)\varphi _{n\sigma }(𝐫),$$
(4)
where $`N_\sigma `$, the number of particles with each spin, is determined by distributing the $`N`$ particles over the $`N`$ lowest eigenfunctions $`\varphi _{n\sigma }`$, such that the total energy is minimized. If the number of filled spin up levels is equal to that of filled spin down levels, the system is unpolarized ($`N_{}=N_{}`$), while it aquires a finite spin polarization if they are different. The maximum possible spin polarization is realized if all particles have the same spin (say spin up), so that $`n_{}=n`$ and $`n_{}=0`$.

In the case of a fully polarized system (e.g., saturated ferromagnetism), the N lowest eigenvalues of the spin up Kohn-Sham equation are all lower than the lowest eigenvalue of the spin down Kohn-Sham equation, so that no eigenfunction of the latter is occupied. The full density is thus obtained from the selfconsistent solution of

$$\left[\frac{\mathrm{}^2}{2m}^2+v_{s,}(𝐫)\right]\varphi _n(𝐫)=ϵ_n\varphi _n(𝐫)$$
(5)
with

$$n(𝐫)=n_{}(𝐫)=\underset{n}{\overset{N}{}}\varphi _n^{}(𝐫)\varphi _n(𝐫).$$
(6)
In this case the spin up Kohn-Sham potential of SDFT already determines the full particle density $`n(𝐫)`$. The spin down potential, on the other hand, is a completely arbitrary function of $`𝐫`$, as long as its lowest eigenvalue $`ϵ_0`$ is higher than the $`N`$’th eigenvalue of the corresponding spin up equation, $`ϵ_N`$ (which is precisely the condition for the system to remain fully polarized) . The set of spin densities $`\{n_{}=n,n_{}=0\}`$ thus does not uniquely determine the set of potentials $`\{v_s,v_s\}`$, because there trivially is an infinite number of possible spin down potentials. Since spin densities and wave functions are still in one-to-one correspondence, this implies that the potentials of SDFT are not determined completely by either the KS Slater determinant $`\mathrm{\Phi }`$ or the full many-body ground state $`\mathrm{\Psi }`$.

In the more general case of a partially polarized system the above simple argument for nonuniqueness breaks down, but it is not hard to show that partial polarization does not restore uniqueness. Consider a given representation of the ground state KS Slater determinant in terms of single-particle orbitals $`\varphi _n`$. In ordinary DFT, based on the density only, each of these orbitals determines the corresponding KS potential (up to the irrelevant additive constant $`ϵ_n`$ ) according to

$$v_s(𝐫)=ϵ_n\frac{\widehat{t}\varphi _n(𝐫)}{\varphi _n(𝐫)},$$
(7)
where $`\widehat{t}`$ is the single-particle kinetic energy operator. Eq. (7) is simply the Kohn-Sham equation solved for $`v_s(𝐫)`$. If one now tries to repeat this reasoning in the spin-dependent case, one must solve the Kohn-Sham equation of SDFT for the external potentials $`v_{s,}(𝐫)`$ and $`v_{\sigma ,}(𝐫)`$ in terms of the single-particle orbitals $`\varphi _n(𝐫)`$ and $`\varphi _n(𝐫)`$,

$$v_{s,\sigma }(𝐫)=ϵ_{n\sigma }\frac{\widehat{t}\varphi _{n\sigma }(𝐫)}{\varphi _{n\sigma }(𝐫)}.$$
(8)
Here we restricted us to collinear spin configurations, so that the magnetic field can be expressed in terms of the spin-dependent potentials as $`B_z(𝐫)=[v_{s,}(𝐫)v_{s,}(𝐫)]/\mu `$, where $`\mu =q\mathrm{}/(2mc)`$ is the Bohr magneton.

Again, any of the single-particle orbitals determines the corresponding potential up to a constant, which remains unspecified. But, unlike in the previous case one can choose only one constant freely by adjusting the zero of energy. That is, given $`\varphi _n(𝐫)`$, say, one can determine $`v_{s,}(𝐫)`$ up to an irrelevant constant, but then the zero of energy has been fixed, and $`\varphi _n(𝐫)`$ determines $`v_{s,}(𝐫)`$ only up to a (now physically relevant) constant, rigidly shifting spin up levels with respect to spin down levels. Consequently, the KS Slater determinant does not completely determine the spin-dependent potentials. By appealing to the second part of the HK theorem one readily concludes that the same holds for the many-body wave function $`\mathrm{\Psi }`$.

The only caveat to the above construction is that for the original KS determinant to remain the ground state after adding the free constant to $`v_{s,}`$, this shift must not induce the presence of new (previously unoccupied) single-particle orbitals in the Slater determinant. As above, this is guaranteed if the KS N-particle ground state is separated from the first excited state by a gap, and the shift of $`v_{s,}`$ is sufficiently small compared to this gap.

Adding a constant term, $`\overline{v}`$, to $`v_{s,}(𝐫)`$ amounts to adding $`\overline{v}`$ to the electrostatic potential $`v_s(𝐫)`$ and $`\overline{v}/\mu `$ to the magnetic field $`B_z(𝐫)`$. The nonuniqueness associated with the fact that there is only one zero of energy, but two constants to be determined in Eq. (8), has thus brought us back to a particular example of systematic nonuniqueness, in which the linear combination in Eq. (2) is given by $`\mathrm{\Delta }v=\overline{v}`$ and $`\mathrm{\Delta }𝐁=(\overline{v}/\mu )𝐮_z`$. In general we expect that systematic nonuniqueness is the only generic source of nonuniqueness, while accidental nonuniqueness requires certain ‘pathological’ features, such as complete spin polarization.

Reflecting on the previous examples, what all cases of nonuniqueness, accidental or systematic, have in common is that the wave functions do not change if certain changes are made to the potentials. This implies a degree of ’rigidity’ or ’incompressibility’ of the former with respect to the latter. The relation to incompressibility is, in fact, more than a mere analogy: an essential feature of all nontrivial examples of nonuniqueness (i.e., all except the additive constant one can choose in the electrostatic potential) is a gap between ground and first excited state, as a consequence of which sufficiently small changes to the potentials cannot bring down an excited state to have lower energy than the ground state. The existence of a gap between ground and first excited state, on the other hand, is tantamount to incompressibility in the usual (quantum liquid) sense of the word. Such incompressibility thus constitutes a necessary (but not sufficient) condition for nontrivial nonuniqueness.

In the remainder of this paper we look at several applications of SDFT in order to determine what changes become necessary in view of the above results. First of all, in the standard Kohn-Sham formulation of SDFT the nonuniqueness of $`v_{s,\sigma }`$ is not critical, since it arises from the external potentials (which are known from the beginning), whereas the xc potentials remain unique. Fields in which significant changes are required if the first map is no longer available, include (i) the calculation of excited states, (ii) the determination of exact KS potentials from numerically exact densities, (iii) perturbative strategies to construct approximate density functionals, (iv) the optimized-effective potential method of SDFT, and (v) several recent generalizations of DFT.

(i) Let us consider excited states first. A standard argument  to justify the calculation of excited states in DFT (albeit not by direct identification of the Kohn-Sham eigenvalues with true excitation energies) is that these, too, are functionals of the ground-state density; afterall, according to the combined map $`v[\mathrm{\Psi }[n]]`$ this density determines the full potential, which, via Schrödinger’s equation, determines all eigenstates, ground state as well as excited ones. Clearly, this argument is only correct in the original formulation of DFT, in which the combined map is available. The above discussion shows that it cannot be used to justify attempts to employ SDFT for the calculation of excited states .

(ii) Recently several methods for constructing accurate KS potentials, using as input highly precise numerical densities (obtained from Monte Carlo or configuration interaction calculations), have been proposed . Uniqueness of the resulting KS potentials is taken for granted in these prescriptions. Again, this holds only in the original formulation of DFT, but not in the (much more widely used) SDFT. The generalization of these methods to determine the KS potentials of SDFT must thus be reconsidered. Even the xc potentials, which are unique, are not completely free of this problem: in an exact calculation the nonuniqueness of the KS potentials must of course cancel in the final result, but this is not guaranteed to be the case in an actual (necessarily approximate) implementation, e.g., using finite basis sets. Special care must thus be taken in such calculations to ensure that the xc potential is not contaminated by nonuniqueness.

(iii) A useful method for obtaining LDA-type approximations for the xc functional is based on perturbative diagrammatic calculations of the xc energy of the electron gas . The results of such calculations are necessarily obtained in terms of the potentials acting on the unperturbed system. As long as the first map of the HK theorem holds, they are thus automatically also functionals of the densities. Obviously this ceases to be true if that map is no longer available. It follows from the above that the SDFT counterpart of this procedure works only because the electron gas does not have a gap between its ground and first excited state. Attempts to apply this procedure to systems with a gap (such as the superconducting electron gas considered in Ref. ) thus face unanticipated and fundamental problems in their generalization to the spin-resolved case.

(iv) In the optimized-effective potential (OEP) method one minimizes an energy functional that is an explicit functional of the single-particle orbitals and only an implicit functional of the densities. In order to calculate the functional derivative with respect to the densities, one then employs the functional chain rule, to write

$`v_{xc}[n](𝐫)={\displaystyle \frac{\delta E_{xc}[\{\varphi _n\}]}{\delta n(𝐫)}}=`$ (9)
$`{\displaystyle 𝑑r^{}𝑑r^{\prime \prime }\underset{n}{}\frac{\delta E_{xc}[\{\varphi _n\}]}{\delta \varphi _n(𝐫^{})}\frac{\delta \varphi _n(𝐫^{})}{\delta v(𝐫^{\prime \prime })}\frac{\delta v(𝐫^{\prime \prime })}{\delta n(𝐫)}}+c.c.`$ (10)
The first functional derivative on the right-hand side can be calculated explicitly, the second is easily found from perturbation theory, and the last is usually identified with an inverse response function. In SDFT the potentials are not unique, so that this inverse does not exist for all densities, and hence the solution to the OEP integral equation following from Eq. (10) is not guaranteed to be the correct $`v_{xc}`$. The problem here arises from using the KS potential as an intermediate step in performing the functional derivative of the orbital functional, the xc potential itself is of course still unique. This intermediate step can be avoided by directly using $`\delta \varphi _n/\delta n`$ on the right-hand side, but the resulting equation will be more complicated to solve than the standard OEP equation.

(v) Finally, consider generalized density-functional theories. The fact that even in the oldest and most widely used of all generalizations of DFT, namely SDFT, uniqueness of the potentials is not guaranteed, casts serious doubts on uniqueness claims in other, more recent, generalizations of DFT. It can be suspected that in, e.g., relativistic DFT, time-dependent DFT, polarization-dependent DFT, DFT for superconductors, etc, similar problems to those discussed here in the context of SDFT are lurking. Indeed, in current-DFT  this is known to be the case, and examples of both systematic and accidental nonuniqueness can be found .

In summary, it has been shown, both by explicit examples and by general considerations, that, unlike in ordinary DFT, in SDFT and other generalized DFTs the effective and the external potentials are not uniquely determined by the spin densities alone. Nonuniqueness can arise accidentally, via special features of the ground state, and systematically, via extensive constants of motion. In both cases it depends on the presence of a gap in the spectrum of the system, and on the resulting incompressibility. As a consequence many previous applications of SDFT and other generalized DFTs in condensed matter physics and quantum chemistry must be critically reexamined.

Acknowledgments KC thanks L. N. Oliveira, M. Lüders, and E. K. U. Gross for valuable discussions, the Physics Department at the University of Columbia, Missouri, for generous hospitality, and the FAPESP for financial support. GV acknowledges support from NSF Grant DMR 9706788.