text:
# PX Andromedae: Superhumps and variable eclipse depthbased on observations obtained at Rozhen National Astronomical Observatory, Bulgaria and at Hoher List, Germany

## 1 Introduction

PX And is probably one of the most complicated SW Sex stars (Thorstensen et al. th91 (1991); Hellier & Robinson hel94 (1994); Still et al. still (1995)). Thorstensen et al. (th91 (1991)) reported shallow eclipses with highly variable eclipse depth and repeating with a period of $``$0$`\stackrel{d}{.}`$1463533. The authors assumed a steady-state accretion disk effective temperature distribution $`T_{\mathrm{eff}}r^{3/4}`$ and $`q0.46`$, and adjusted the inclination and disk radius to match the width and depth of the mean eclipse. They obtained $`i73.8^{}`$ and $`r_\mathrm{d}0.6R_{L_1}`$. No other system parameters estimations have been published and all authors used the above values. Apart from the distinctive characteristics of SW Sex stars (reviewed recently by Hellier helrev (2000)), PX And shows some other interesting peculiarities. Patterson (patt99 (1999)) reported PX And to show simultaneously “negative” and “positive” superhumps, and signals with typical periods of 4–5 days. This implies that PX And most probably possesses an eccentric and tilted accretion disk. In the view of the expected stream overflow (Hellier & Robinson hel94 (1994)) all this would result in very complex accretion structures.

In this paper we report the results of a photometric study of PX And.

## 2 Observations and data reduction

Photometric CCD observations of PX And were obtained with the 2.0-m telescope of Rozhen Observatory. A Photometrics 10242 CCD camera and a Johnson $`V`$ filter were used. The atmospheric conditions were stable with seeing $``$2$`^{^{\prime \prime }}`$ and a 2$`\times `$2 binned camera was used. This resulted in $`3`$ pixels per FWHM and $``$13 s read-out dead-time. In total 8 runs were obtained in 2000 and 2001. The exposure time used was between 20 and 40 s. In addition, two unfiltered runs were obtained with the 1-m telescope at Hoher List Observatory. Some details of the observations are given in Table 1. After bias and flat-field corrections the photometry was done with the standard DAOPHOT aperture photometry procedures (Stetson ste (1987)). The magnitude of PX And was measured relatively to the star And1-5 ($`V`$=15.908) and And1-8 ($`V`$=16.613) served as a check (Henden & Honeycutt comp (1995)). The runs are shown in Fig. 1.

## 3 Results

### 3.1 The ephemeris

The observations cover 14 usable eclipses. The corresponding eclipse timings were determined by fitting Gaussian to the lower half of the eclipses and are given in Table 1. Combining them with the eclipse timings published by Li et al. (li (1990)), Hellier & Robinson (hel94 (1994)), Andronov et al. (andr (1994)) and Shakhvskoy et al. (shak (1995)) we determined the following orbital ephemeris:

$$T_{\mathrm{min}}[HJD]=49238.83662(14)+0\stackrel{d}{.}146352739(11)E.$$
(1)
### 3.2 Periodogram analysis

#### 3.2.1 Low-frequency periodicities

The mean out-of-eclipse magnitude of PX And shows large variations. The five October 2000 runs obtained within 6 nights suggest that the mean magnitude is modulated with a period of $``$5 days and a full amplitude of $``$0.5 mag. After the eclipses have been masked, the October 2000 data were searched for periodic brightness modulations by computing the Fourier transform. The complex spectral window of the data results in a power spectrum which is completely dominated by the aliases of the strong 5-days wave (Fig. 2$`a`$). The computed power spectrum was deconvolved with the spectral window by using the CLEAN algorithm (Roberts et al. rob (1987)). The CLEANed spectrum is shown in Fig. 2$`b`$. Apart from the strong low-frequency peak there are a number of peaks remaining after the cleaning. PX And has already been reported to show “negative superhumps” with a period of 0$`\stackrel{d}{.}`$1415 and cycles with periods of $`45`$ days (Patterson patt99 (1999); Center for Backyard Astrophysics (CBA) web site111 http://cba.phys.columbia.edu/results/pxand). The origin of these signals is not fully understood, but they are believed to be caused by a retrograde precession of an accretion disk which is tilted with respect to the orbital plane. In this model the relation between the orbital period $`P_{\mathrm{orb}}`$, the “negative superhumps” period $`P_{\mathrm{SH}}^{}`$ and the precession period of the disk $`P_{\mathrm{prec}}`$ is given by

$$\frac{1}{P_{\mathrm{prec}}}=\frac{1}{P_{\mathrm{SH}}^{}}\frac{1}{P_{\mathrm{orb}}}.$$
(2)
If the peak at $`P0\stackrel{d}{.}142`$ is associated with the “negative superhumps”, then the expected precession period of the disk is $``$4$`\stackrel{d}{.}`$8. Thus, the low-frequency wave seen in the data is most likely a manifestation of the disk precession. We note that CLEAN is an iterative algorithm, which is initiated with the assumption that the strongest peak in the power spectrum corresponds to a real periodicity in the data. Thus, to apply CLEAN algorithm to our PX And data we have to assume that the $``$5-days modulation is real. Although the baseline of our observation is short, they strongly indicate the presence of a $``$5-days modulation. In addition, CBA observations also reveal such a modulation. Thus, our assumption seems reasonable.

To search for additional weaker signals, the $``$0$`\stackrel{d}{.}`$142 and $``$4$`\stackrel{d}{.}`$8 waves were subtracted from the data by performing a non-linear least-squares multi-sinusoidal fit with unknown amplitudes and phases. With observations over a baseline of only 6 days, $`P_{\mathrm{prec}}`$ could not be accurately determined; because of this $`P_{\mathrm{prec}}`$ was calculated from Eq. (2). The CLEANed spectrum of the residuals (not shown) shows that all strong peaks are removed. The only exception is the peak at $`P0\stackrel{d}{.}207`$. The appearance of this peak is quite unusual because it could not be related to any clock in a system with $`P_{\mathrm{orb}}`$=0$`\stackrel{d}{.}`$146352739. This peak should be regarded with caution since the corresponding period is roughly equal to the mean duration of the runs. Peaks corresponding to periods of the order of the duration of the runs or harmonics of a period of 1 day may appear in the power spectrum because of the color difference between the target and the comparison star. The full amplitude of the 0$`\stackrel{d}{.}`$207 signal is however $``$0.15 mag; too large to be caused by the second-order extinction. Moreover, the second-order extinction in $`V`$ band is negligible. To check if this peak is not an artifact of the cleaning we applied several ways of pre-whitening. We subtracted fits with the $``$4$`\stackrel{d}{.}`$8 and $``$0$`\stackrel{d}{.}`$142 periods only, a fit with both periods but allowing them to vary in a narrow interval around peak center, etc. The results show that the peak at $`P`$0$`\stackrel{d}{.}`$207 is always present. We further estimated what is the probability for this peak to appear in the power spectrum by chance alone. This estimate is complicated by the fact that the noise in PX And light curves is red rather than white (Sec. 3.2.2). Because of this we applied the following procedure. We have generated 1000 artificial red noise time series sampled exactly as the observations, applied CLEAN to each of them and then estimated the probability to see peaks as strong as the peak at $``$0$`\stackrel{d}{.}`$207 by chance. The artificial time series were generated by the method suggested by Timmer & König (tim (1995)). This method allows one to generate random realizations of a given red noise process. The only input information needed is the power spectrum of the process. We used the mean red noise power spectrum shape derived in Sec. 3.2.2, which ensures that the artificial time series have on average the same power spectrum as the PX And light curves. The analysis of the CLEANed spectra of the artificial time series shows that the peak at 0$`\stackrel{d}{.}`$207 is statistically significant by more than 99%. Further, we examined whether the peak at 0$`\stackrel{d}{.}`$207 could be introduced in the CLEANed spectrum by the presence of the other two periodic signals, in conjunction with the red noise and the gaps in the data. To do this we added to the artificial red noise time series two sinusoids with periods of 4$`\stackrel{d}{.}`$8 and 0$`\stackrel{d}{.}`$142, and amplitudes and phases as determined from the fit of the original data. None of the CLEANed spectra of these time series show strong peaks around 0$`\stackrel{d}{.}`$207. This suggests that the presence of the 4$`\stackrel{d}{.}`$8 and 0$`\stackrel{d}{.}`$142 signals could not be the reason for the appearance of the 0$`\stackrel{d}{.}`$207 peak in the CLEANed spectrum of PX And data.

In order to refine the periods we again performed a non-linear least-squares multi-sinusoidal fit to the data, but this time the periods were also allowed to vary in a narrow interval around the peaks maxima. $`P_{\mathrm{SH}}^{}`$ and $`P_{\mathrm{prec}}`$ are not independent, so in the function used to fit the data they were coupled through Eq. (2). Thus, only $`P_{\mathrm{SH}}^{}`$ and $`P_1`$=0$`\stackrel{d}{.}`$207 were varied. In that way the following periods were determined: $`P_{\mathrm{SH}}^{}`$=0$`\stackrel{d}{.}`$142 $`\pm `$0.002, $`P_1`$=0$`\stackrel{d}{.}`$207 $`\pm `$0.004 and $`P_{\mathrm{prec}}`$=4$`\stackrel{d}{.}`$8 with semi-amplitudes of 0.086, 0.076 and 0.256 mag, respectively. The errors of $`P_{\mathrm{SH}}^{}`$ and $`P_1`$ are estimated from the half-width at half-maximum of the central peak of the spectral window. A reliable error estimate of the 4$`\stackrel{d}{.}`$8 period is not possible mainly because the data set is short and covers only $``$1 cycle. The best fit is shown in Fig. 1 and the mean wave shapes of the ”negative superhumps” and 0$`\stackrel{d}{.}`$207 signal in Fig. 3. The CLEANed spectrum of the residuals (Fig. 2$`c`$) shows no peaks which implies that all peaks in the power spectrum are successfully accounted for by the fit with the detected periods.

Periodogram analysis was performed to the three runs obtained in 2001, after the subtraction of the mean of each run. The Lomb-Scargle periodogram (Scargle scar (1982)) is shown in Fig. 2$`d`$. The sparse distribution of the runs resulted in a complex alias structure of the periodogram. It is clear, however, that the strongest peak corresponds to a “negative superhump” with $`P_{\mathrm{SH}}^{}`$0$`\stackrel{d}{.}`$141. This value of $`P_{\mathrm{SH}}^{}`$ gives $`P_{\mathrm{prec}}`$3$`\stackrel{d}{.}`$8. The fit with these two periods is shown in Fig. 1. The corresponding semi-amplitudes are 0.069 and 0.21 mag. The 0$`\stackrel{d}{.}`$207 period is not detected in 2001. Note, however, that in both runs obtained on 18 and 19 October there is a slight increase of the out-of-eclipse magnitude, which might be a manifestation of a periodicity around 0$`\stackrel{d}{.}`$207.

#### 3.2.2 High-frequency periodicities

To search for high-frequency periodicities we subtracted the low-frequency signals from the data and calculated the power spectra. The periodograms show many peaks, but no coherent oscillations were detected. In double logarithmic scale the power decreases linearly with frequency (Fig. 4) and this is usually interpreted as a result of the flickering (Bruch bruch (1992)). The power spectrum can be described by the equation

$$P(f)=\frac{\alpha }{1+(2\pi \tau f)^\gamma }$$
(3)
with $`\gamma `$ being the slope of the linear part and $`\tau `$ may serve as an estimation of the mean duration of the dominant structures in the light curves. Equation (3) was fitted to the mean power spectrum yielding $`\gamma `$=1.71 and $`\tau `$=6 min. We also determined $`\tau `$ by means of autocorrelation function as described in Kraicheva et al. (1999b ). The mean value of $`\tau `$$`\pm `$1.5 min is slightly lower, but consistent with that determined from the fit of the mean power spectrum. The mean standard deviation in the light curves after the subtraction of the low-frequency signals is 0.061 $`\pm `$0.008 mag. This value is consistent with the standard deviation found in the light curves of TT Ari during the “negative superhumps” regime (Kraicheva et al. 1999a ). After switching to “positive superhumps”, the standard deviation in the TT Ari light curves decreased by a factor of 2. It would be interesting to compare this quantity in the novalikes showing only “negative” or “positive superhumps” to see if there is a systematic difference.

### 3.3 Variable eclipse depth

As first noticed by Thorstensen et al. (th91 (1991)) the eclipse depth in PX And is highly variable. The October 2000 runs leave the impression that the eclipse depth is modulated with the precession cycle (Fig. 5 upper row). This however could be a visual effect: the eclipses look deeper when they coincide with the superhumps minima. To verify this we have subtracted the best fit from the data and show the residuals in relative flux in the lower row of Fig. 5. The eclipse depth modulation with the precession period is evident even after the subtraction. Figure 6 shows the eclipse depth as a function of the mean out-of-eclipse magnitude. We fitted the eclipse depths with a sinusoid with the precession period. This shows that the eclipses are deepest at the minimum of the precession cycle. The mean eclipse depth is $``$0.56 ($``$0.89 mag) and the amplitude of the variation is $``$0.11. Thus, the eclipses we see are on average much deeper than the value of 0.5 mag reported by Thorstensen et al. (th91 (1991)). In Figs. 5 and 6 one can also find hints that in 2001 PX And was $``$0.2 mag brighter than in 2000.

## 4 Discussion

The origin of the “negative superhumps” is still an open question. The most plausible model is based on the assumption of a retrograde precession of a tilted accretion disk. While the ”positive superhumps” have been simulated numerically since the work of Whitehurst (wh (1988)), there are some difficulties to simulate the ”negative superhumps”. In fact, all attempts to simulate “negative superhumps” failed in the sense that they were not able to produce a significant tilt starting from a disk lying in the orbital plane (Murray & Armitage mur (1998); Wood et al. wood (2000)). Once tilted, however, the accretion disk starts precessing in retrograde direction (Larwood et al. lar (1996); Wood et al. wood (2000)) thus giving two additional clocks: $`P_{\mathrm{prec}}`$ in inertial co-ordinate system and $`P_{\mathrm{SH}}^{}`$ in a system co-rotating with the binary. $`P_{\mathrm{prec}}`$ is typically a few days and thus $`P_{\mathrm{SH}}^{}`$ is slightly shorter than the orbital period. The angle at which the accretion disk is seen from the Earth is modulated with the precession period, giving brightness modulations with $`P_{\mathrm{prec}}`$. Taking into account foreshortening and limb-darkening, one can estimate the disk tilt needed to produce the observed amplitude of the 4$`\stackrel{d}{.}`$8 modulation in PX And. With a limb-darkening coefficient $`u=0.6`$ the tilt angle is between 2.5 and 3, depending on the assumed system inclination. For comparison the simulations of Wood et al. (wood (2000)) were performed with a tilt angle of 5.

The mechanism generating the “negative superhump” light itself is however rather uncertain. According to the simulations of Wood et al. (wood (2000)) the two opposite parts of the disk which are most displaced from the orbital plane are tidally heated by the secondary. The heating is maximal when these parts point to the secondary and this happens twice per superhump cycle. Wood et al. (wood (2000)) find that the tidal stress is asymmetric with respect to the disk mid-plane and the disk surface facing the orbital plane is more heated. Thus, if the disk is optically thick the observer sees only one side and a modulation with the superhump period is observed. Let us define precession phase $`\varphi _{\mathrm{prec}}`$ to be zero at the maximum of the 4$`\stackrel{d}{.}`$8 cycle, i.e. when the disk is most nearly face-on. The relative phasing between the periodic signals determined from the fit shows that the maxima of the superhumps coincide with the eclipses at $`\varphi _{\mathrm{prec}}`$0.13. Figure 7 shows a sketch of the system configuration in eclipse at $`\varphi _{\mathrm{prec}}=0.13`$. The expected position of the superhump light source (SLS) according to Wood et al. (wood (2000)) and the line of nodes are also shown. The dashed line marks the part of the disk which lays below the orbital plane. According to the simulations of Wood et al. (wood (2000)) the ”negative superhumps” maxima occur when the region labelled as ”SLS” in Fig. 7 lies on the line connecting the two system components, which means that the superhumps maxima will coincide with the eclipses at $`\varphi _{\mathrm{prec}}=0`$. Our results show that in PX And this is not exactly the case and the superhumps maxima are observed $`0.13P_{\mathrm{SH}}^{}`$ later.

There is another mechanism which could generate “negative superhumps” and at the same time account for the delay. Patterson (patt99 (1999)) noticed that most of the SW Sex novalikes show “negative superhumps” and this could naturally explain the accretion stream overflow thought to be responsible for the SW Sex phenomenon (Hellier & Robinson hel94 (1994)). In a system with a precessing tilted accretion disk, the amount of gas in the overflowing stream will vary with the “negative superhumps” period. Correspondingly, the intensity of the spot (shown in Fig. 7) formed where the overflowing stream re-impacts the disk should be modulated and might be the “negative” SLS. A similar model was proposed by Hessman et al. (hess (1992)) for the origin of ”positive superhumps”. The maximal overflowing is expected to take place when the most displaced part of the disk points to the accretion stream. The intensity of the re-impact spot will however reach its maximum later when the gas in the stream reaches the re-impact zone (the moment shown in Fig. 7). This means that the gas in the accretion stream which at the moment shown in Fig. 7 is near the re-impact zone, has passed over the disk edge earlier, when the system configuration was different and the most displaced part of the disk pointed to the accretion stream.

Let the gas in the overflowing stream travels from the outer disk edge to the re-impact zone in the time $`\mathrm{\Delta }t`$. Then, because in co-ordinate system rotating with the binary the line of nodes precesses with period $`P_{\mathrm{SH}}^{}`$, at the superhumps maxima the displaced part of the disk will not point to the accretion stream but will be rotated with respect to it by an angle $`\theta 360\mathrm{\Delta }t/P_{\mathrm{SH}}^{}`$ deg. Figure 7 shows that in the case of PX And $`\theta 0.15\times 360=54^{}`$. Therefore, in order to observe the maxima of the “negative superhumps” at $`\varphi _{\mathrm{prec}}`$0.13 and $`\varphi _{\mathrm{orb}}`$0 the travel time from the outer disk edge to the re-impact zone has to be $`30`$ min. This is slightly longer than expected (Warner & Peters war (1972)), but we have to take into account that the velocity of the incoming gas is most probably reduced by the first impact with the accretion disk. This model for the origin of the “negative” SLS can be independently checked by tracking the intensity of the high-velocity $`s`$-wave observed in the spectra of SW Sex novalikes.

The origin of the 0$`\stackrel{d}{.}`$207 signal is quite uncertain. The corresponding peak in the power spectrum could not be a result of an amplitude modulation of the “negative superhumps”. If this is the case, the power spectrum should show two additional peaks. The secondary crosses the line of nodes twice per superhump cycle. The time between three consecutive passages coincides within 1.5% with the observed period. It is however not clear how this would produce periodic brightness modulations, moreover with the observed amplitude of $``$0.15 mag.

The most puzzling observation of PX And is its highly variable eclipse depth. Our observations suggest that the eclipse depth is modulated with the precession period (Figs. 5 and 6), indicating that this phenomenon could be related to the disk precession. In the precessing tilted accretion disk model, in some particular combinations of $`q`$ and $`i`$ the eclipsed area of the disk is expected to vary with the precession period and thus to modulate the eclipse depth. We have simulated eclipses of a tilted precessing disk which showed that the eclipse depth indeed varied but with a very low amplitude of 2-3%, which is insufficient to explain the observations.

The relatively shallow eclipses show that the accretion disk in PX And is not totally eclipsed. In this case, if there is a significant constant light source which is totally eclipsed, the eclipse depth will be modulated since different fraction of the total light emitted by the system will be eclipsed at different precession phases. The deepest eclipses are expected at the precession cycle minima, as observed. As the constant light source must be eclipsed, its most likely location is close to the white dwarf. The white dwarf itself could not be the source of the constant light as it does not contribute significantly to the total system light. The most obvious candidate is the inner hot part of the disk. Since the brightness modulation with the precession cycle comes from more or less pure geometrical considerations it is not clear why the emission of the inner parts of the disc could be constant. One possible explanation might be that the inner part of the disk is not tilted or is less tilted than the outer part. In this case the emission from the inner disk will be constant. Another source of the constant light could be the emission from the boundary layer between the disk and the white dwarf surface. This however is not very likely since the temperature of the boundary layer is thought to be rather high to emit significantly in the optical wavelengths. Of course, there exists the possibility that the eclipse depth modulation with the precession cycle is an artifact. If some flickering peaks occur during the eclipse (which is likely because the accretion disk in PX And is not totally eclipsed), this will decrease the observed eclipse depth. Figure 8 shows that some eclipses are indeed badly affected by the flickering. Because of the small number of eclipses covered this might introduce a spurious modulation. But since the correlation between eclipse depth and mean magnitude is good, we suggest that the eclipse depth modulation with the precession cycle is real.

There is another way to modulate the eclipse depth. As the SLS is non-uniformly distributed over the disk surface and the accretion disk is not totally eclipsed, there are two extreme possibilities when the SLS is (i) totally and (ii) never eclipsed. In both cases the eclipse depth will depend on the actual superhump light at the moment of eclipse, but will have opposite behavior. If the SLS is totally eclipsed, then the deepest eclipses will be observed when superhump maxima coincide with the eclipses. If the SLS is never eclipsed, the eclipse depth will be affected in a way analogous to the veiling of the absorption lines in T Tau stars. In these stars the additional emission in the continuum reduces the observed depth of the absorption lines. Then the eclipses will be deepest when superhump minima coincide with the eclipses. Using the parameters of the periodic signals obtained from the multi-sinusoidal fit, we calculated the total contribution of the modulated part of the flux to the unmodulated one, $`F_{\mathrm{mod}}/F_{\mathrm{const}}`$, at the moments of the eclipses. $`F_{\mathrm{const}}`$ is modulated with the precession cycle and hence may be thought to be constant on the superhump time scale. The $``$0$`\stackrel{d}{.}`$207 signal was included in the calculations for 2000. The reason is that the eclipse depth variation observed in the 2000 data is rather high to be explained if only the “negative” SLS is involved. Figure 9 shows the eclipse depth $`d_\mathrm{e}`$ as a function of $`F_{\mathrm{mod}}/F_{\mathrm{const}}`$. Although the scatter is large, one should keep in mind that any flickering near the mid-eclipses will decrease the eclipse depth and will move the points in Fig. 9 downward. If a small number of eclipses are observed then this could easily obscure any correlation. Despite the large scatter, there are some hints of an inverse correlation in Fig. 9. This suggests that the SLS is not eclipsed. In Fig. 9 is also shown an arbitrary curve defined by the equation:

$$d_\mathrm{e}=\frac{d_{\mathrm{e},0}}{1+F_{\mathrm{mod}}/F_{\mathrm{const}}}$$
(4)
which gives the eclipse depth as a function of the additional modulated light in the case when this light is not eclipsed. Here, $`d_{\mathrm{e},0}`$ is the eclipse depth if $`F_{\mathrm{mod}}=0`$.

Although mathematically the veiling could explain the eclipse depth changes, there are many uncertainties with this interpretation. Most of them are related to the unknown place of origin of the periodic signals. If the conclusion of Wood et al. (wood (2000)) about where the “negative superhumps” are generated is correct, then the “negative” SLS should be at least partially eclipsed. This suggests that the stream overflow could be the mechanism generating the “negative superhumps” in PX And. Moreover, in a system with grazing eclipses the re-impact zone is not eclipsed. Apart of this, the origin of the $``$0$`\stackrel{d}{.}`$207 signal is unknown and it is not clear if its source is eclipsed or not. Thus, the question why the eclipse depth in PX And varies remains open and we have only pointed to some of the possible explanations.

Thorstensen et al. (th91 (1991)) have estimated $`q`$, $`i`$ and the disk radius of PX And by fitting a mean eclipse. It should be noted however that the observations of Thorstensen et al. (th91 (1991)) point to an average eclipse depth of $`0.5`$ mag while our observations show much deeper eclipses. It is clear that the highly variable shape and depth of eclipses could affect significantly any attempt to estimate the system parameters of PX And. Because of this we have not attempted such an estimation from the eclipse profiles. It would also not be correct to define so called ”mean” eclipse and to analyze it with the eclipse mapping technique. Moreover, a large part of the accretion disk in PX And is not eclipsed. The eclipse mapping algorithm could be easily modified to allow analysis of tilted accretion disks. To use this technique however one will need to average at least several eclipses at given precession phase in order to reduce the influence of flickering and other noise.

Generally, if a given cataclysmic variable shows “positive superhumps” one could estimate $`q`$ using the existing $`ϵ^+(q)`$ relations ($`ϵ^+=(P_{\mathrm{sh}}^+P_{\mathrm{orb}})/P_{\mathrm{orb}}`$). Patterson (patt99 (1999)) reported that PX And shows “positive superhumps” with $`ϵ^+`$0.09 and this could be used to estimate $`q`$. Recently, Montgomery (mon (2001)) published an analytic $`ϵ^+(q)`$ expression which takes into account the pressure effect. Applying this relation (Eq. (8) in Montgomery mon (2001)) to PX And we obtain $`q`$0.27. Although this value of $`q`$ is more plausible for a system showing ”positive superhumps”, one should keep in mind that the peak in the power spectrum of the CBA data is very weak. The ”negative superhumps” are much more confidently detected, but unfortunately an $`ϵ^{}(q)`$ relation for the ”negative superhumps” has not been established yet. Patterson (patt99 (1999)) noticed that $`ϵ^{}0.5ϵ^+`$. Our periodogram analysis yields $`ϵ^{}0.029`$, hence $`ϵ^+0.058`$. From Montgomery (mon (2001)) relation, one obtains $`q0.18`$. This is not an unrealistic value but we cannot be sure that the relation $`ϵ^{}0.5ϵ^+`$ holds true for PX And.

Based on the analysis of our new photometric observations and on the results of Patterson (patt99 (1999)) we can with no doubt place the novalike PX And among the cataclysmic variables showing permanent ”negative superhumps”. The observations of the star in 2000 show two unusual features: the presence of another superhump with a period of $``$0$`\stackrel{d}{.}`$207 and a modulation of the eclipse depth with the precession period of the accretion disk. Our 2001 observations are however limited in time coverage and it is difficult to say if these two phenomena are typical of the star or are a peculiarity of this particular data set only. As a suggestion for future work we recommend performing a multi-site photometric observations covering several precession cycles. If the modulation of the eclipse depth with the precession period is confirmed, this phenomenon could be further studied by analysis of the eclipse profiles at different precession phases.

###### Acknowledgements.

We are grateful to the referee John Thorstensen for his valuable comments and suggestions. HB wishes to thank Prof. Wilhelm Seggewiss for generously allocating time at Hoher List. The work was partially supported by NFSR under project No. 715/97.


format: text

text:
# Excitonic optical spectrum of semiconductors obtained by time-dependent density-functional theory with the exact-exchange kernel

## Abstract

Applying the novel exact-exchange (EXX) Kohn-Sham method within time-dependent density-functional theory, we obtained the optical absorption spectrum of bulk silicon in good agreement with experiments including excitonic features. Analysis of the EXX kernel shows that inclusion of the Coulomb coupling of electron-hole pairs and the correct long-wavelength behavior in the kernel is crucial for the proper description of excitonic effects in semiconductors.

The calculation of optical spectra of solids has been traditionally one of the most important and challenging areas of first-principles material investigations. The interaction of excited electrons and holes plays a crucial role in optical excitations, however properly incorporating such effects ab initio is theoretically and computationally very difficult. In the past several years there has been a major advance in the field BSE98  based on the solution of the Bethe-Salpeter equation (BSE) Hanke79  for the two-particle Green’s function starting from quasiparticles obtained within the GW approximation (GWA) Hybertsen85 , but unfortunately the calculation scheme is by nature complicated and demanding.

One attractive alternative is extending density-functional theory (DFT) to the time-dependent (TD) case TDDFT . TDDFT is a well-founded theory for the treatment of electronic excitations in general and, unlike the GWA-BSE route which is presently implemented as a post-DFT method, provides a computational scheme based directly on Kohn-Sham (KS) one-particle states. Moreover, once the KS one-particle states are given, the electronic linear response properties are completely determined by the orbital-independent Hartree kernel $`F_\mathrm{H}(𝐫,𝐫^{})1/|𝐫𝐫^{}|`$ AU  and the dynamic exchange-correlation kernel

$$F_{\mathrm{xc}}(𝐫,𝐫^{};tt^{})\frac{\delta v_{\mathrm{xc}}(𝐫;t)}{\delta n(𝐫^{};t^{})}.$$
(1)
Thus, in principle, one only needs to find a good approximation to the exchange-correlation kernel for practical TDDFT applications, which represents a significant conceptual and computational simplification. In practice, however, the conventional adiabatic local-density approximation (LDA) and generalized-gradient approximation (GGA) kernels are known to be inadequate for the study of optical excitations in solids. A primary example is their incorrect long-wavelength behavior for insulators which has been recently emphasized in several contexts Aulbur96 ; Ghosez97 ; Reining02 .

Deficiencies of the LDA and GGA in fact appear already on the level of the ground-state exchange-correlation energy $`E_{\mathrm{xc}}`$ and potential $`v_{\mathrm{xc}}(𝐫)\delta E_{\mathrm{xc}}/\delta n(𝐫)`$. With respect to energetics, for example, the LDA and GGA inherently fail to describe the quasi-two-dimensional electron gas 2D-egas  and the van der Waals interactions vdW . For $`v_{\mathrm{xc}}`$ and the corresponding KS eigenvalue spectrum, LDA and GGA $`v_{\mathrm{xc}}`$ incorrectly decay exponentially rather than as $`1/r`$ in finite systems and accordingly their highest occupied orbital energies are too high and unoccupied orbitals do not exhibit Rydberg series, and for solids their band gaps are too small. The deficiencies of the LDA and GGA $`v_{\mathrm{xc}}`$ in particular pose significant practical difficulties for the study of electronic excitations within the DFT formalism.

In this regard, recent realizations of multi-dimensional KS exact-exchange (EXX) Stadele97 ; Gorling99 ; Ivanov99  and effective EXX Bylander95 ; Grabo97 ; KLI99 ; DellaSala01  methods provide an interesting opportunity. The orbital-based self-interaction-free EXX methods not only provide realistic local multiplicative KS exchange potentials and KS eigenvalue spectra for molecules Grabo97 ; KLI99 ; Gorling99 ; Ivanov99 ; DellaSala01  but also give band gaps of semiconductors in good agreement with experiments Stadele97 ; Bylander95 . We have also recently shown DFT01  that the EXX spectrum at the one-particle level without any previously applied post-DFT modification such as a quasiparticle shift Levine89  gives a very good description of the absorption spectrum of semiconductors excluding excitonic features resulting from two-particle interactions, which supports the notion that KS eigenvalue differences represent well-defined approximations for excitation energies Gorling96 ; Umrigar98 . However, the remaining excitonic character in the spectrum is not properly treated by the LDA or GGA kernel, so a complete set of DFT methods for the investigation of electronic optical excitations in solids is still missing.

Given the importance of proper inclusion of excitonic effects in the study of electronic optical excitations and the encouraging performance of the EXX method for band structures of semiconductors, we consider in this work a TDDFT scheme based on the EXX kernel. At the formal level, in contrast to LDA and GGA kernels, the EXX kernel is nonlocal in real space, depends explicitly on the frequency, and has the correct long-wavelength limit behavior Ghosez97 ; FXEXX02 . Taking bulk silicon as the representative semiconductor system, we show that the EXX kernel indeed provides optical spectra of semiconductors in good agreement with experiments including excitonic effects and further analyze the origin of the excitonic peak and the validity of locality approximations in space and time. We additionally make comparisons between EXX-TDEXX and GWA-BSE, which renders new insight into the TDDFT approach in general Gonze99 .

To obtain the optical spectrum of a solid, we computed the frequency $`\omega `$-dependent macroscopic dielectric function $`ϵ_\mathrm{M}(\omega )`$ which can be written in terms of the modified full linear response matrix $`\widehat{\chi }`$ as Hanke79

$$ϵ_\mathrm{M}(\omega )=1\underset{q0}{lim}F_\mathrm{H}(𝐪)\widehat{\chi }(𝐪;\omega )|_{𝐆=𝐆^{}=0},$$
(2)
where $`𝐪`$ is the photon momentum, and $`F_\mathrm{H}(𝐪)F_\mathrm{H}(𝐆,𝐆^{},𝐪)=\delta _{𝐆,𝐆^{}}4\pi /|𝐪+𝐆|^2`$. Within TDDFT, $`\widehat{\chi }`$ is completely determined by the KS linear response matrix $`\chi _0`$ and the kernel matrix $`F_{\mathrm{xc}}`$ according to

$$\widehat{\chi }(𝐪;\omega )=\left[1\chi _0(𝐪;\omega )\{\widehat{F}_\mathrm{H}+F_{\mathrm{xc}}(𝐪;\omega )\}\right]^1\chi _0(𝐪;\omega ),$$
(3)
where $`\widehat{F}_\mathrm{H}=0`$ for $`𝐆=𝐆^{}=0`$ and $`\widehat{F}_\mathrm{H}=F_\mathrm{H}`$ otherwise. At the independent-particle level, $`\widehat{\chi }=\chi _0`$, the two-particle interaction effects due to $`\widehat{F}_\mathrm{H}`$ and $`F_{\mathrm{xc}}`$ are both ignored, while only the $`F_{\mathrm{xc}}`$ part is ignored at the time-dependent Hartree level. As is apparent from Eqs. (2) and (3), the long-wavelength $`q0`$ behavior of the “head” ($`𝐆=𝐆^{}=0`$) and “wing” ($`𝐆=0`$ and $`𝐆^{}0`$ or vice versa) elements of $`F_{\mathrm{xc}}`$ (in insulators, $`𝒪(1/q^2)`$ for the head and $`𝒪(1/q)`$ for the wings Ghosez97 ; FXEXX02 ) plays an important role.

The computation of $`\chi _0`$ is straightforwardly done with KS one-particle states and energies, and consequently $`F_{\mathrm{xc}}`$ remains as the only important quantity to be determined. So far, for $`F_{\mathrm{xc}}`$, the adiabatic LDA kernel $`F_{\mathrm{xc}}^{LDA}=\delta ^2E_{\mathrm{xc}}^{LDA}/(\delta n\delta n^{})`$ has been almost exclusively employed. It is local in real space and frequency-independent, which results in a reciprocal-space representation independent of $`q`$ and $`\omega `$, $`F_{\mathrm{xc}}^{LDA}(𝐆,𝐆^{},𝐪;\omega )=F_{\mathrm{xc}}^{LDA}(𝐆𝐆^{})`$. While numerical advantageous, this drastic simplification is known to be problematic for the study of optical responses in solids. For example, the head and wings of the LDA kernel are incorrectly nondivergent in the $`q0`$ limit and this defect cannot be corrected by a semilocal GGA kernel.

Instead, we adopted in this work the EXX kernel,

$$\begin{array}{cc}\hfill F_\mathrm{x}^{EXX}& (𝐆,𝐆^{},𝐪;\omega )=\underset{𝐆_1,𝐆_2}{}\chi _0^1(𝐆,𝐆_1,𝐪;\omega )\hfill \\ & \times H_\mathrm{x}(𝐆_1,𝐆_2,𝐪;\omega )\chi _0^1(𝐆_2,𝐆^{},𝐪;\omega ),\hfill \end{array}$$
(4)
which is nonlocal in real space and explicitly depends on the frequency. The expression of the EXX kernel “core” $`H_\mathrm{x}`$, which was interpreted in the many-body diagrammatic picture as the first-order self-energy and vertex corrections to the irreducible polarizability $`\stackrel{~}{\chi }=\chi _0+\chi _0F_{\mathrm{xc}}\stackrel{~}{\chi }`$ Tokatly01 , has been presented in Ref. FXEXX02  for the case of periodic insulators. We have further shown that $`F_\mathrm{x}^{EXX}`$ has the $`q0`$ behavior of the exact $`F_{\mathrm{xc}}`$ and thus rectifies one serious deficiency of the LDA and GGA kernels. The full expression for $`H_\mathrm{x}`$ is rather lengthy and we list here only the three resonant terms convention :

$$\begin{array}{cc}\hfill H_\mathrm{x}^{A\text{res}}(𝐪;\omega )& \frac{2}{\mathrm{\Omega }}\underset{as𝐤}{}\underset{bt𝐤^{}}{}\left[\frac{a𝐤|e^{i(𝐪+𝐆)𝐫}|s𝐤+𝐪s𝐤+𝐪;b𝐤^{}|\widehat{w}_C|t𝐤^{}+𝐪;a𝐤t𝐤^{}+𝐪|e^{i(𝐪+𝐆^{})𝐫}|b𝐤^{}}{(ϵ_{a𝐤}ϵ_{s𝐤+𝐪}+\omega +i\delta )(ϵ_{b𝐤^{}}ϵ_{t𝐤^{}+𝐪}+\omega +i\delta )}\right];\hfill \\ \hfill H_\mathrm{x}^{B\text{res}}(𝐪;\omega )& \frac{2}{\mathrm{\Omega }}\underset{abs𝐤}{}\left[\frac{a𝐤|e^{i(𝐪+𝐆)𝐫}|s𝐤+𝐪b𝐤|\widehat{\mathrm{\Sigma }}_x\widehat{v}_x|a𝐤s𝐤+𝐪|e^{i(𝐪+𝐆^{})𝐫}|b𝐤}{(ϵ_{a𝐤}ϵ_{s𝐤+𝐪}+\omega +i\delta )(ϵ_{b𝐤}ϵ_{s𝐤+𝐪}+\omega +i\delta )}\right]\hfill \\ & +\frac{2}{\mathrm{\Omega }}\underset{ast𝐤}{}\left[\frac{a𝐤|e^{i(𝐪+𝐆)𝐫}|s𝐤+𝐪s𝐤+𝐪|\widehat{\mathrm{\Sigma }}_x\widehat{v}_x|t𝐤+𝐪t𝐤+𝐪|e^{i(𝐪+𝐆^{})𝐫}|a𝐤}{(ϵ_{a𝐤}ϵ_{s𝐤+𝐪}+\omega +i\delta )(ϵ_{a𝐤}ϵ_{t𝐤+𝐪}+\omega +i\delta )}\right],\hfill \end{array}$$
(5)
where two is the spin factor, $`\mathrm{\Omega }`$ is the crystal volume, $`\{a,b\}`$ are valence bands, $`\{s,t\}`$ are conduction bands, the matrix elements $`i;j|\widehat{w}_C|l;m`$ are four-index integrals defined as

$$\begin{array}{cc}\hfill i𝐤+𝐪;j𝐤^{}|\widehat{w}_C|l𝐤^{}+𝐪;m𝐤& 𝑑𝐫𝑑𝐫^{}\varphi _{i𝐤+𝐪}^{}(𝐫)\varphi _{j𝐤^{}}^{}(𝐫^{})\hfill \\ & \times w_C(𝐫,𝐫^{})\varphi _{l𝐤^{}+𝐪}(𝐫)\varphi _{m𝐤}(𝐫^{}),\hfill \end{array}$$
(6)
$`\widehat{\mathrm{\Sigma }}_x`$ denotes an orbital-dependent exchange operator of the form of the Hartree-Fock-like exchange operator but constructed from the KS orbitals, i.e.,

$$i𝐤+𝐪|\widehat{\mathrm{\Sigma }}_x|j𝐤+𝐪\underset{a𝐤^{}}{}i𝐤+𝐪;a𝐤^{}|\widehat{w}_C|a𝐤^{};j𝐤+𝐪,$$
(7)
and $`\widehat{v}_x`$ is generated by the orbital-independent EXX KS potential $`v_x(𝐫)`$. Here $`w_C(𝐫,𝐫^{})`$ yielding $`\widehat{w}_C`$ denotes the generalized Coulomb interaction, e.g., $`1/|𝐫𝐫^{}|`$ for the bare Coulomb interaction.

Compared with the standard BSE approach, although we also have to perform the four-index Coulomb integrals of Eqs. (6) and (7), we can calculate optical spectra without diagonalizations in the space of occupied and unoccupied single-particle states. Due to the huge number of $`𝐤`$-points involved in the optical spectrum calculations, the size of the matrix to be diagonalized can be very large, and the fact that we can avoid the diagonalization process altogether in principle represents a significant numerical advantage. In addition, note that $`F_\mathrm{x}^{EXX}`$ is free of Coulomb singularities FXEXX02  and thus the formal validity of the EXX TDDFT approach does not depend on the consideration of the thermodynamic limit. In the BSE equation, on the other hand, integrable singularities appear which require special numerical care BSE98 .

We have implemented the full dynamic EXX kernel in reciprocal space employing the plane wave basis. The accuracy of the code has been carefully checked by numerically testing if the calculated exchange kernel acts as a functional derivative of the exchange potential with respect to the electron density at $`\omega =0`$, $`\delta v_\mathrm{x}^{EXX}=F_\mathrm{x}^{EXX}\delta n`$. To obtain $`F_\mathrm{x}^{EXX}`$ we first carried out self-consistent EXX ground-state calculations at the experimental lattice constant of Si, $`5.43\AA `$, and generated the KS potential Stadele97 ; DFT01 . Ten special-$`𝐤`$ points and an orbital kinetic energy cutoff of 12.5 Ha have been used. In the response calculation step, we solved the KS equations once more at a larger number of $`𝐤`$-points and obtained KS orbitals and eigenvalues as the input for the construction of $`\chi _0`$ and $`F_\mathrm{x}^{EXX}`$. We adopted a shifted uniform $`𝐤`$-mesh, for which we employed up to $`22\times 22\times 22`$ $`𝐤`$-points for $`\chi _0`$ and up to $`9\times 9\times 9`$ $`𝐤`$-points for $`F_\mathrm{x}^{EXX}`$. The kinetic energy cutoff of 10 Ha and $`\delta `$ = 0.15eV have been used, and ten conduction bands were included.

The first question we addressed with small $`𝐤`$-points sets (e.g. $`5\times 5\times 5`$) was whether the adiabatic approximation, defined as the $`\omega =0`$ limit of Eqs. (4) and (5), is justified. Although the nature of the approximation within TDEXX is different from that of the BSE approach in which the static screened interaction is taken, we found that the adiabatic EXX kernel generates an overall similar spectrum as the dynamic one. Computing $`F_\mathrm{x}^{EXX}`$ only at a single $`\omega `$ results in a significant reduction of computational workload, and thus the converged calculation was performed in the adiabatic approximation.

Figure 1 shows the absorption spectrum of Si obtained at the single-particle EXX level and by taking into account two-particle interaction effects via $`\widehat{F}_\mathrm{H}`$ plus $`F_{\mathrm{xc}}^{LDA}`$ (EXX+TDLDA) and via $`\widehat{F}_\mathrm{H}`$ plus $`F_\mathrm{x}^{EXX}`$ (EXX+TDEXX). Compared with the much discussed LDA (and LDA+TDLDA) spectrum, which is incorrectly shifted to the lower frequency region by about 1 eV due to its well-known band gap underestimation, the single-particle EXX absorption edge and second ($`E_2`$) peak are in excellent agreement with the experimental curve Herzinger98  due to the realistic EXX band structure DFT01 . However, the first ($`E_1`$) peak originating from electron-hole attractions is much underestimated at the single-particle EXX level, and it is not recovered by taking into account electron-hole interactions via TDLDA, demonstrating the failure of $`F_{\mathrm{xc}}^{LDA}`$ to describe excitonic effects. The EXX+TDEXX spectrum, on the other hand, is in excellent agreement with the measured data: The absorption strength at the $`E_1`$ peak region is correctly enhanced while that of the higher frequency region is reduced, which shows that $`F_\mathrm{x}^{EXX}`$ indeed provides the correct description of the important excitonic effects.

To obtain the EXX+TDEXX spectrum in Fig. 1, we used a slightly modified bare Coulomb interaction as $`w_C`$. Employing the bare Coulomb interaction resulted in a collapse of the spectrum due to too strong long-range Coulomb interaction of electron-hole pairs at different $`𝐤`$-points. We therefore cut off the long-range Coulomb interaction for these pairs, or set to zero the contributions to the integrals of Eq. (6) for which $`(𝐆+𝐤𝐤^{})`$ lies within the first Brillouin zone knekp . Devising a more systematic strategy might be desirable in the future.

Now, to understand the underlying mechanism of the encouraging EXX+TDEXX result, we analyze $`F_\mathrm{x}^{EXX}`$ by investigating the role of different contributions to $`H_\mathrm{x}`$ as shown in Fig. 2. First, by taking only $`H_\mathrm{x}^B`$, we obtain a spectrum resembling that of the one-particle EXX but shifted upward by about 0.3 eV. No excitonic feature arises at this level. This spectrum is similar to that of the GWA, which is in agreement with previous predictions Gonze99 ; Tokatly01  that GWA results should be recovered with the $`H_\mathrm{x}^B`$ part of the screened EXX kernel. However, note that we have used an only slightly modified bare Coulomb interaction, and, while the GWA spectrum results from independent quasiparticles that involve ($`N\pm 1`$)-electron excitations, the DFT “quasiparticles” spectrum already represents electron-number conserving ($`N`$-electron) one-particle excitations at the time-independent level Gorling96 ; Umrigar98 ; DFT01 . Second, if we take only $`H_\mathrm{x}^A`$, the spectrum is strongly enhanced at the lower frequency region ($``$ 4 eV) while it is much reduced at the higher frequency region ($``$ 4 eV). The excitonic $`E_1`$ peak appears predominantly while the $`E_2`$ peak is slightly red-shifted by about 0.2 eV. This shows that the Coulomb coupling of electron-hole pairs is precisely the origin of the excitonic $`E_1`$ peak. Finally, we set the head and wings of $`F_\mathrm{x}^{EXX}`$ to zero and thus make them $`𝐪`$-independent as in the case of the LDA or GGA. By doing so, we obtain a spectrum similar to that of the EXX+TDLDA without excitonic peaks or whatsoever. This demonstrates that taking into account the correct $`q0`$ behavior of $`F_{\mathrm{xc}}`$ is crucial.

We finally comment on other recent TDDFT works that have also obtained optical spectra of semiconductors in good agreement with experiments. First, de Boeij et al. have left the domain of conventional DFT and resorted to TD current-DFT deBoeij01 . Albeit they had to introduce an empirical energy shift and a prefactor for the polarization functional, their approach may represent a different route to include the space nonlocality discussed above via a macroscopic functional. However, our result shows that employing the current response derived from a macroscopic functional is not necessary if one properly incorporates the space nonlocality in the microscopic density response. In this context, recent work of Reining et al. which was based on the idea of mapping the BSE to the TDDFT is in agreement with ours Reining02 . But in their scheme they first invoked the GWA quasiparticle shift and second employed an empirical static scaled Coulomb kernel as $`F_{\mathrm{xc}}`$. Although their approach is not rigorous, we note that their simple static kernel nevertheless exhibits a correct $`q0`$ divergence.

In summary, we have reported calculations of the optical absorption spectrum of Si with the novel TDEXX approach. Of particular interest was the nature of the dynamic and nonlocal orbital-based EXX kernel which exhibits the $`q0`$ behavior of the exact $`F_{\mathrm{xc}}`$ and successfully generated excitonic features in the spectrum. We showed that including Coulomb coupling of KS electron-hole pairs at different $`𝐤`$-points and the correct long-wavelength behavior of the exchange-correlation kernel is crucial, while the adiabatic approximation can be justified.

We thank L. Reining, L. X. Benedict, and E. L. Shirley for helpful discussions. This work has been supported by the Humboldt Foundation (Y.-H.K.), the Deutsche Forschungsgemeinschaft, and the Fonds der Chemischen Industrie (A.G.).


text:
# High-frequency fluctuational conductivity in YBa2Cu3O7-δ

We present measurements of the fluctuational excess conductivity at 48 GHz and 24 GHz in YBa2Cu3O7-δ films above the transition temperature. The measurements depart from the gaussian prediction for finite-frequency fluctuational conductivity. We focus on the region not too close to $`T_c`$, where the real part of the excess conductivity drops much faster than the prediction of the gaussian model. We calculate the dynamic excess conductivity within a Ginzburg-Landau approach. In the calculation we insert a short-wavelength cutoff of the order of the inverse coherence length, in order to suppress high-momentum modes. The excess conductivity of all samples measured can be described very well by the modified model.
PACS numbers: 74.40.+k, 74.76.-w, 78.70.Gq
J. Low Temp. Phys., to be published.

High-$`T_c`$ superconductors are model systems for the study of superconducting fluctuations. High transition temperatures, short coherence lengths and anisotropy contribute to unveil fluctuative regimes that are hindered in conventional superconductors. In fact, in the latter the superconducting fluctuations are rather well described in terms of relatively simple models, such as (focussing on the electrical transport) the Aslamazov and Larkin theory (in a microscopic framework) or gaussian models.
A large part of the experiments on transport properties in the fluctuational region have been performed by measuring the dc electrical conductivity. While this approach benefits from the higher precision of dc measurements, it is intrisically unable to give information on the characteristic time of the fluctuations, as opposed to the dynamic conductivity. It is then interesting to examine appropriate tools for the analysis of the dynamic conductivity.
Most of such models have been developed to hold in a temperature region neither too close to $`T_c`$ (that is, outside of the critical region), nor too far from $`T_c`$ (where the single $`q`$=0 fluctuation mode is supposed to fully describe the system). Due to the wide temperature region where fluctuation effects are observed in high-$`T_c`$ cuprates, these simple models have been recently reexamined. In this paper we focus on the analysis of the temperature region not too close to $`T_c`$, introducing a simple Ginzburg-Landau (GL) model for the finite-frequency excess conductivity. We then discuss our microwave data for the excess conductivity in light of the proposed extension to the gaussian theory.
We briefly sketch the theoretical approach. We derive the (a,b)-plane complex conductivity within the correlation function formalism, using the GL functional written for a three-dimensional uniaxial anisotropic superconductor, and limited to terms containing $`\psi ^2`$. We do not attempt here to include the quartic term since its effect becomes relevant only very close to $`T_c`$. The dynamic of the order parameter, introduced through the time-dependent GL theory, is governed by $`w=\omega \tau _0/ϵ`$, where $`\omega /2\pi `$ is the measuring frequency, $`ϵ`$=ln$`(T/T_c)`$ is the reduced temperature, and $`\tau _0`$ is the temperature-independent GL relaxation time. The complex conductivity is obtained by averaging the current operator with respect to the noise as $`J_{ab}\left(t\right)=\frac{\mathrm{}e^{}}{m_{ab}}\frac{d^3q}{\left(2\pi \right)^3}q_{ab}C\left[𝐤=𝐪\frac{e^{}}{\mathrm{}}𝐀\left(t\right);t,t\right]`$, where $`C`$ is the order-parameter correlator, $`𝐀`$ is the harmonic vector potential, $`e^{}=2e`$ is twice the electronic charge, and $`m_{ab}`$ is the in-plane mass of the pair. Since $`\tau _0`$20 fs, a small-$`w`$ expansion is appropriate. The common approach involves integration of the expression for the current over all $`q`$ values, since the $`q`$=0 mode is the most diverging close to $`T_c`$. The result is the well-known gaussian result by Schmidt:

$$\mathrm{\Delta }\sigma _g=\mathrm{\Delta }\sigma _{1,g}+\mathrm{i}\mathrm{\Delta }\sigma _{2,g}\mathrm{\Delta }\sigma _0\left[\left(1\frac{w^2}{16}\right)+\mathrm{i}\frac{w}{6}\right]$$
(1)
where $`\mathrm{\Delta }\sigma _0=\frac{e^2}{32\mathrm{}\xi _c(0)ϵ^{1/2}}`$ is the gaussian, dc excess conductivity, and $`\xi _c(0)`$ is the zero temperature out-of-plane GL coherence length. When temperature raises above $`T_c`$, the relative weight of the $`q`$ 0 modes increases. However modes with $`q`$ higher than, e.g., the inverse coherence length should not be allowed. The customary way to take into account this feature is to put a cutoff in the integration of the current operator (see an extended discussion in Ref.). Using a single cutoff such that $`\sqrt{\underset{j=a,b,c}{}[q_j^{max}\xi _j(0)]^2}<\mathrm{\Lambda }`$, the dynamic (cutoffed) excess conductivity becomes (for small $`w`$):

$`\mathrm{\Delta }\sigma _{1,c}\mathrm{\Delta }\sigma _0\left\{{\displaystyle \frac{2}{\pi }}\left[atnK{\displaystyle \frac{K}{\left(1+K^2\right)^2}}\left(1+{\displaystyle \frac{5}{3}}K^2\right)\right]+o(w^2)\right\}`$ (2)
$`\mathrm{\Delta }\sigma _{2,c}\mathrm{\Delta }\sigma _0{\displaystyle \frac{w}{6}}{\displaystyle \frac{2}{\pi }}\left[atnK+{\displaystyle \frac{K}{\left(1+K^2\right)^3}}\left(K^4{\displaystyle \frac{8}{3}}K^21\right)\right]`$ (3)
where $`K=\mathrm{\Lambda }/ϵ^{\frac{1}{2}}`$. We note that the frequency dependence remains unaltered with respect to the gaussian expression, while the temperature dependence changes due to the cutoff effects. In Fig.1a we plot a simulation of $`\mathrm{\Delta }\sigma _{1,c}/\mathrm{\Delta }\sigma _{1,g}`$ and $`\mathrm{\Delta }\sigma _{2,c}/\mathrm{\Delta }\sigma _{2,g}`$ as a function of $`ϵ`$, for cutoff number $`\mathrm{\Lambda }`$=0.5, and neglecting all but the leading term in the expressions for $`\mathrm{\Delta }\sigma `$. It is immediately seen that the cutoff effects, if present, are best visible in $`\mathrm{\Delta }\sigma _1`$.
We have measured the microwave response at 23.9 and 48.2 GHz in two YBa2Cu3O7-δ films, grown by different techniques. In all cases the film was the end wall of a cylindrical metal cavity, resonating in the TE011 mode. Details of the technique have been given elsewhere. Measurements of the quality factor yielded the effective surface resistance $`R`$. The measured films were thinner than 2500 Å, so that in the transition region $`R\rho _1/d`$, where $`\rho _1`$ is the real part of the resistivity and $`d`$ is the film thickness.
In one sample (measured at 48.2 GHz) we have carefully measured the frequency shift near the transition, in order to obtain the effective surface reactance $`X\rho _2/d`$, where $`\rho _2/d`$ is the imaginary part of the resistivity. The data are reported in Fig.2a. We discuss with some length these data, in order to illustrate the main aspects of the analysis.
The complex excess conductivity $`\mathrm{\Delta }\sigma `$ is defined as follows:

$$\rho _1+\mathrm{i}\rho _2=\frac{1}{\sigma _n+\mathrm{\Delta }\sigma }$$
(4)
where $`\rho _n=\frac{1}{\sigma _n}`$ is the normal state resistivity, that we take as the linear extrapolation from high temperatures. The inversion of Eq.4 gives $`\mathrm{\Delta }\sigma _1`$ and $`\mathrm{\Delta }\sigma _2`$, that are reported as a function of $`T`$ in Fig.2b. Such data are first used to evaluate $`T_c`$, using the scaling theories prediction $`\mathrm{\Delta }\sigma _1(T_c)=\mathrm{\Delta }\sigma _2(T_c)`$. Since $`\mathrm{\Delta }\sigma _2`$ is vanishingly small $``$ 1 K above $`T_c`$, we can approximate $`\mathrm{\Delta }\sigma _1\frac{1}{\rho _1}\frac{1}{\rho _n}`$ not too close to $`T_c`$. The result is reported in Fig.2c, and compared to the exact inversion of the $`\rho `$ data. As can be seen, above 1.01 $`T_c`$ the two procedures yield the same result. This feature can be put on quantitative grounds by plotting the percentage deviation of the approximate from the exact derivation of $`\mathrm{\Delta }\sigma _1`$, as in Fig.2d. It is seen that the deviation is largely below 5% for $`TT_c>`$0.5 K.
From this discussion of these data, we can assert that the $`w`$ 0 regime for $`\mathrm{\Delta }\sigma _1`$ is reached in a wide temperature range, even at our highest measuring frequency of 48.2 GHz.
On the basis of the previous discussion, we have extracted in the second sample $`\mathrm{\Delta }\sigma _1`$ with the approximate expression, and we have used the inflection point of $`\rho _1`$ to estimate $`T_c`$ (we note that the precise choice of $`T_c`$ does not significantly affect the data analysis in the range of interest here, $`ϵ>10^2`$). In Fig.1b we report the ratio $`\mathrm{\Delta }\sigma _1/\mathrm{\Delta }\sigma _0`$, having used $`\xi _c(0)`$ as a scale factor in order to get the correct limit at low $`ϵ`$. These data show a clear evidence for the failure of the simple gaussian approximation, which predicts a constant value for all $`ϵ`$. In the same Fig.1b, continuous lines are fits performed through Eq.2, in good agreement with the data. We can conclude that our measurements of the excess conductivity, performed at microwave frequencies, show clear departures from the gaussian prediction, and that the correction is well described by the introduction of a cutoff in the spectrum of accessible momenta.
We are indebted to T.Mishonov, D.Neri, R.Raimondi for helpful discussions.


text:
# The relative ages of LMC old clusters, and the case of NGC 1841

## 1. Introduction

The LMC is the closest disk galaxy (at $`50`$ kpc) with a large ensamble of star clusters. At least $`11`$ of them can be considered  *bona fide* counterparts (e.g. Olszewski et al. 1996) of Galactic globular clusters (GGC). Their absolute ages are comparable to those of the Milky Way metal-poor ones (Brocato et al. 1996; Olsen et al. 1998; Johnson et al. 1999), however, their kinematics is strikingly different, since they show a definite disk-like rotation of $`80\mathrm{km}\mathrm{s}^1`$ (Freeman et al. 1983). It is of great interest to understand the reasons of the simultaneous formation of the two systems in such dynamically different protogalaxies. A key input to this problem is the age dispersion of LMC globulars.

## 2. The Relative ages of LMC old clusters, and the case of NGC 1841

The CMDs are based on archival HST/WFPC2 data (GO 5897). Data reduction and calibration will be presented in a forthcoming paper. In Fig. 1 (left panel) the trend of $`\delta (VI)_{\mathrm{@}2.5}`$ vs. \[Fe/H\] is compared to a set of theoretical isochrones (see R99 for details), and the six LMC clusters fall inside the $`2`$ Gyr strip defined by metal-poor GGCs. We thus confirm that most LMC clusters are coeval with the Milky Way globulars, and more importantly, for the first time, we can constrain the age dispersion of the LMC flattened halo, which results to be as low as $`\sigma _t0.5`$ Gyr. According to Fig. 1 (left panel), NGC 1841 is $`2`$ Gyr younger than the mean Galactic/LMC sample, a fact that is confirmed by a direct comparison of its CMD to that of Hodge 11 (right panel of Fig. 1). Indeed, one can see that the young age is revealed by the HB-TO difference as well, since the NGC 1841 turnoff is brighter than that of Hodge 11.

## 3. Conclusions

We have found that both the mean age and age dispersion of the LMC GCs is comparable to that of the metal-poor Galactic GCs, and clear evidence has been obtained in favor of a younger age for NGC 1841. Since it is the most metal-poor LMC cluster, it means that the interstellar medium (ISM) where it formed, has had a lower degree of chemical evolution. NGC 1841 could have been formed in a relatively isolated fragment of the proto-LMC, or it could have been part of an independent system now disrupted. Further support to the latter hypothesis could be the fact that NGC 1841 is the farthest cluster from the LMC center ($`10`$ kpc) and that its radial velocity is incompatible with the rotation of the old LMC halo (Freeman et al. 1983). The fact that LMC globulars formed in a rotating protogalaxy with perfect synchronization with the GGCs, seems to require a joint formation scenario, possibly induced by the Milky Way early SF (e.g. Taniguchi et al. 1999).

## References

Brocato, E., et al. 1996, MNRAS, 282, 614

Freeman, K.C., Illingworth, G., & Oemler, A. 1983, ApJ, 272, 488

Johnson, J. A., et al. 1999, ApJ, 527, 199

Olsen, K. A. G., et al. 1998, MNRAS, 300, 665

Olszewski, E.W., Suntzeff, N.B., & Mateo, M. 1996, ARA&A, 34, 511

Palmieri, R., et al. 2002, A&A, in press

Rosenberg, A., Saviane, I., Piotto, G., & Aparicio, A. 1999, AJ, 118, 2306 (R99)

Taniguchi, Y., Trentham, N., & Ikeuchi, S. 1999, ApJ, 526, 13

VandenBerg, D. A., et al. 2000, ApJ, 532, 430


text:
# Early-Phase Spectra of ”Hypernova” SN 2002ap

## 1 Introduction

The very powerful supernova (SN) explosions called “hypernovae” (Iwamoto et al., 1998) are perhaps the most energetic events in the current universe. Their kinetic energy is $`10^{52}`$ ergs, which is about 10 times larger than of normal SNe. SN 1998bw, probably associated with GRB 980425 (Galama et al., 1998), was the first recognized hypernova. Although classified as a Type Ic SN (SN Ic), the spectrum of SN 1998bw showed unusual broad spectral features in the early phase. The derived velocities were as large as 30,000 km s-1. The kinetic energy of SN 1998bw was deduced to be 3 – 5 $`\times `$ 1052 ergs (Iwamoto et al., 1998; Nakamura et al., 2001). Other Type Ic hypernovae were subsequently recognized based on their broad absorption features. These include SNe 1997ef (e.g., Iwamoto et al. 2000), 1998ey (Garnavich, Jha, & Kirshner, 1998), and the recently discovered 2002bl (Filippenko et al., 2002). Moreover, SN 1997dq was suggested to be a hypernova from its similarity to SN 1997ef at later phases (Matheson et al., 2001). The kinetic energy deduced for SN 1997ef was about $`1.8\times 10^{52}`$ ergs (Mazzali, Iwamoto, & Nomoto 2000), which was smaller than that of SN 1998bw, but still significantly larger than normal SNe. Hypernovae also occur among SNe IIn, for example, SN 1997cy, one of the brightest SNe ever discovered and characterized by a large kinetic energy, $`3\times 10^{52}`$ ergs (Turatto et al., 2000), and SN 1999E (Cappellaro, Turatto, & Mazzali 1999). The range of properties of hypernovae, and their possible association with GRBs are also topics of great interest (e.g., Nomoto et al. 2001).

SN 2002ap was discovered in M74 on January 29.4 UT by Y.Hirose (Nakano et al., 2002). Spectra obtained on January 30 – 31 showed very broad absorption features, resembling those in SNe 1998bw and 1997ef. Therefore, SN 2002ap was suggested to be a hypernova (Kinugasa et al., 2002; Meikle et al., 2002; Gal-Yam & Shemmer, 2002). SN 2002ap is the nearest hypernova discovered to date, hence it was a good target even for small telescopes. Various multi-wavelength observations of SN 2002ap have been performed (summarized in Gal-Yam, Ofek, & Shemmer 2002). NIR photometry was also performed frequently at Gunma Astronomical Observatory (GAO) (E. Nishihara et al. 2002, in preparation).

We report here on the optical spectroscopy of SN 2002ap performed at GAO. We compare the spectra of SN 2002ap with those of other hypernovae. We also used the parameterized SN synthetic-spectrum code SYNOW to perform line identification and deduce velocity information from the spectra in the early phase, when very broad features resulting from line-blending were present.

## 2 Observations and Data Reduction

Observations of SN 2002ap started on the night of the discovery report. Nightly early phase spectral observations were carried out at GAO and Bisei Astronomical Observatory (BAO) for about 40 days, from discovery to solar conjunction. Spectra were obtained with Gunma Compact Spectrograph (GCS) attached to the classical-Cassegrain focus of the GAO 65cm telescope (F/12) and with a spectrograph attached to the BAO 1.01-m telescope (F/12). The spectral resolving power ($`\lambda /\mathrm{\Delta }\lambda `$) of GCS is about 500 at 5500 Å and the wavelength coverage is from 3800 to 7600 Å. The BAO spectrograph covers from 4700 to 7000 Å with a resolving power of about 1000. Observational parameters are listed in Table 1. Since all the spectra were obtained with narrow slits (2′′ on GCS and 2.6′′ at BAO, compared to a typical seeing of 2-3′′) and the sky conditions were generally non-photometric, the flux was not on an absolute flux scale. The slit was aligned in E-W direction at GAO and in the N-S direction at BAO, which in both cases was not the parallactic angle. Consequently, the blue end of the spectra could be significantly affected at advance epochs, as the observations were done at low altitude. In those cases, we cut the blue end of the spectra (see Figure 1).

Data reduction was performed with IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation.. Calibration for the relative intensity was performed based on the observations of spectrophotometric standard stars. A Fe-Ne-Ar lamp for GAO and a Fe-Ne lamp for BAO were used for the wavelength calibration.

## 3 Results and Discussion

The time-sequence of the calibrated spectra of SN 2002ap is shown in Figure 1. The evolution of the spectra is apparent. As time progresses, most absorption features shift to the red, and the continuum becomes redder.

The first two spectra, taken on the same night, show a blue continuum, with neither deep absorptions nor strong emissions (see also the top of Figure 2). Broad and shallow depressions are seen near 4700, 5700, and 6200 Å. These broad features resemble those of SNe 1998bw and 1997ef in the earliest phases (Patat et al., 2001; Iwamoto et al., 2000). However, SN 2002ap is bluer, possibly because it was caught at an earlier epoch. The fact that the first observations must have been very close to the time of explosion is exemplified by the rapid development of spectral features. The first two spectra also show a broad peak near 5000Å. But in the third spectrum taken the following night, this peak has shifted significantly to the red ($`5100`$Å), indicating a rapid decrease of the material velocity sampled by the spectra, as the photosphere moves inwards.

Photometric observation at Wise Observatory established that the peak of the light curve occurred on February 7.1 and 8.8 UT in the $`B`$ and $`V`$-band, respectively (Gal-Yam et al., 2002). We designate the epochs using the number of days from $`B`$ maximum. Figure 2 also shows the spectrum of SN 2002ap taken on day $`+`$2 (February 9.4 UT), near the time of maximum. We compared this with the near-maximum spectra of other SNe Ic: SNe 1998bw, 1997ef (hypernovae) and 1994I (a normal SN Ic; Nomoto et al. 1994). Based on the line width, the spectra can easily be separated into two groups. One is that of “hypernovae”, including SNe 1998bw, 1997ef, and 2002ap, while the other includes only the normal Type Ic SN 1994I. Therefore, SN 2002ap is easily identified as a hypernova. Moreover, the spectrum of SN 2002ap shows emission-like features around 4600 Å and 5300 Å, resembling the spectrum of SN 1997ef more than that of SN 1998bw.

In Figure 3 we compare similar-looking spectra of SNe 2002ap and 1997ef. We can see that the spectral evolution of SN 2002ap is about 1.5 times faster. The earliest available spectrum of SN 1997ef, on day $``$6, appears to have properties intermediate between the spectrum of SN 2002ap on day $``$1 and that on day $``$6, judging from the progressive redshifting of line features. The near-maximum spectrum of SN 1997ef resembles that of SN 2002ap on day $`+`$2 in all major features. The similarity between the spectrum of SN 1997ef on day $`+`$19 and that of SN 2002ap on day $`+`$15 is striking, both showing minor features of comparable intensities at about 4900, 5900, 6200, 7100 and 7300 Å. Most of the new features visible in the spectrum of SN 1997ef on day $`+`$52, some of which are probably net emissions, are already visible in SN 2002ap on day $`+`$25, although they are not as developed. Both of these spectra suggest that the SNe are making a transition to the nebular phase, while the spectrum of SN 1997ef on day $`+`$27 is still essentially photospheric in nature (Mazzali et al., 2000).

Especially in the earliest spectra, most features are so broad because of line blending that it is not easy to identify which lines contribute to them. Thus we used the parameterized SN spectrum synthesis code SYNOW, described by Fisher (2000), to perform line identification and deduce velocity information. It is simpler and faster than the Monte-Carlo code of Mazzali et al. (2002), and better suited for our objective, as it does not require the construction of realistic SN models. In this paper, the radial-dependent line optical depth above the photosphere is assumed to be proportional to $`(v/v_{\mathrm{ph}})^n`$, where $`v`$ is the expansion velocity and $`v_{\mathrm{ph}}`$ is the photospheric velocity. A reasonable fit to the observations could be obtained using just a few ions: Ca II, O I, Si II, Fe II, Co II and Ni II, consistent with the line identification in SNe 1997ef and 1998bw (Branch, 2001). Synthetic spectra and line identification for day $``$7 and $`+`$2 are shown in Figure 2.

Either a high $`v_{\mathrm{ph}}`$ or a small $`n`$ can make the lines broad. Guided by Branch’s (2001) work, we found that both are necessary to account for the extreme line widths. Similar results were obtained by Mazzali et al. (2002), who used a Monte-Carlo code to synthesize some selected spectra of SN 2002ap. Our best fit value of the power-law index was 3 for the first spectrum. This index value was used for the radial dependence of the opacity for all lines and at all epochs. As shown in Figure 4, the photospheric velocity obtained from the fitting decrease from 35,000 $`\mathrm{km}\mathrm{s}^1`$ on day $`7`$ to 13,000 $`\mathrm{km}\mathrm{s}^1`$ on day $`+`$2, with the uncertainty being about 1000 $`\mathrm{km}\mathrm{s}^1`$ empirically (Branch et al., 2002).

Our fits are reasonable for these epochs, allowing us to identify most absorption lines. However, some problems remain. One is the unwanted “peak” around 4100 Å, which is actually the remainder of continuum. We can remove it by introducing both Ti II and Sc II lines. However, there is no other argument for Sc II in SNe Ic. Such difficulties may show the limitations of SYNOW. They are probably due to the assumption of a pure resonant scattering source function and on the simple power-law dependence of the optical depth on velocity applied to all ions. The hump around 5200 Å is also incorrectly reproduced, probably the consequence of not transferring sufficient flux from the blue ($`4000`$Å) in SYNOW. Finally, we introduced Na I to fit the absorption around 5700 Å. He I $`\lambda `$5876 may also contributes there. The NIR spectra show a significant feature, possibly due to He I $`\lambda `$10830 (Motohara et al., 2002). However, that feature in SN 1994I might have been C I, O I, and/or Si I lines. (Baron et al., 1999; Millard et al., 1999).

Among those absorption lines identified, the Si II $`\lambda \lambda `$6347, 6371 (collectively called $`\lambda `$6355) absorption minimum was easily distinguishable during almost the entire observation period. Therefore, we measured the minimum of the Si II absorption as an alternative method of estimating the photospheric velocity. In Figure 1, we marked with an arrow the minimum of the Si II absorption. The evolution of the photospheric velocity, as deduced from this method, is also shown in Figure 4. At the earliest phases, the typical statistical error was around $`\pm `$ 3000 km s-1. Moreover, strong line blending makes the estimate of the velocity somewhat uncertain at those epochs. Near maximum light, although the minima are clearly visible (the statistical error was about $`\pm `$ 1500 km s-1), the systematic error on the velocity is mostly due to the ambiguity introduced by the uncertain continuum. At late phases, the error (around $`\pm `$ 2500 km s-1) was mainly due to poor statistics. The difference between the velocities from Si II and those from SYNOW fitting may be mostly caused by systematic errors introduced by the uncertainty in the continuum and the heavy line blending.

In Figure 4, we compare the evolution of the photospheric velocity of SN 2002ap, 1998bw and 1997ef (Patat et al., 2001). We adopted for SN 2002ap an explosion date of January 28.9 UT (day $``$9.2), as proposed by Mazzali et al. (2002). The photospheric velocity of SN 2002ap declines faster than for the other two hypernovae, suggesting that the ejecta mass of SN 2002ap is smaller. Although the velocity evolution of the three SNe is initially well fitted with exponential functions, both SNe 1997ef and 1998bw eventually seem to flatten out, SN 1998bw at about 7000 km s-1 and SN 1997ef at about 2000 km s-1. While the measurement for SN 1998bw can be done directly on the spectra, the value for SN 1997ef is derived from model fitting, since, as for SN 2002ap, the Si II line becomes indistinguishable at advanced epochs (see spectra on day $`+`$52 for SN 1997ef and on day $`+`$25 for SN 2002ap in Figure 3). No sign of the flattening for SN 2002ap was seen over the period of our observations, but it might occur later, which may have to be established through detailed spectral modeling.

The rapid evolution of the photospheric velocity and of the spectra of SN 2002ap suggests that SN 2002ap will develop a nebular spectrum earlier than the other hypernovae presented here. In SN 1997ef nebular emission was observed in the Ca II IR triplet on day $`+`$27. Our spectra of SN 2002ap unfortunately do not extend to that region. However, the spectrum on day $`+`$25, which is intermediate between the day $`+`$27 and the day $`+`$52 spectra of SN 1997ef, may show the onset of weak \[O I\] $`\lambda \lambda `$6300, 6364 net emission. E. Nishihara et al. (2002, in preparation) show that a break occurs in the evolution of the IR colors around February 25, suggesting the onset of nebular emission. The first spectra of SN 2002ap after solar conjunction, taken on day 131 and 140, show that it has already entered the nebular phase (Leonard et al., 2002).

The broad spectral features, suggesting high kinetic energy, distinguish hypernovae from normal SNe Ic. However, the observations of SN 2002ap suggest that there is a variety of hypernovae, based on the observational features (spectra, spectral evolution, luminosity, light curves, the association with a GRB, etc.). More observations of hypernovae are necessary in order to understand what physical parameters are responsible for the variations among hypernovae and whether there is a gap in properties between hypernovae and normal SNe Ic (Figure 2 of Nomoto et al. 2002).


text:
# Commensurate-incommensurate transition of cold atoms in an optical lattice

## Abstract

An atomic gas subject to a commensurate periodic potential generated by an optical lattice undergoes a superfluid–Mott insulator transition. Confining a strongly interacting gas to one dimension generates an instability where an arbitrary weak potential is sufficient to pin the atoms into the Mott state; here, we derive the corresponding phase diagram. The commensurate pinned state may be detected via its finite excitation gap and the Bragg peaks in the static structure factor.

Atomic gases are developing into the ultimately tunable laboratory system allowing to study complex quantum phenomena Anglin and Ketterle (2002). Recently, subjecting an atomic Bose-Einstein condensate to an optical lattice, Greiner et al. Greiner et al. (2002) have succeeded in tuning the system through a quantum phase transition separating a superfluid (S) from a Mott insulating (MI) phase. This 3D bulk transition involves weakly interacting bosons and is well understood within the Bose-Hubbard description Fisher et al. (1989); Jaksch et al. (1998): the system turns insulating when the on-site interaction energy $`U`$ becomes of the order of the hopping energy $`J`$. This strong coupling transition is a result of quenching the kinetic energy by a strong lattice potential. Amazingly, by confining the atomic gas to one dimension, the strong coupling limit can be reached without the optical lattice: in 1D, the ratio $`\gamma `$ between the interaction- and kinetic energies per particle scales inversely with the density $`n`$ and thus it is the low-density limit which is interacting strongly (Tonks gas) Petrov et al. (2000). A new instability then appears in the strongly interacting 1D quantum gas: the superfluid groundstate in the homogeneous system turns insulating in the presence of an arbitrarily weak commensurate optical lattice and the S–MI transition changes to a transition of the incommensurate–commensurate type. In this letter, we analyze this new instability and derive the complete phase diagram for the S–MI transition in the limits of both weakly and strongly interacting gases. Remarkably, this goal can be achieved by a mapping to two classic problems, the Bose-Hubbard model (introduced in Ref. Jaksch et al., 1998) and the sine-Gordon problem describing the weakly- and strongly interacting limits of the atomic gas, respectively. Below, we first summarize the main results providing us with the phase diagram shown in Fig. 1; we then proceed with a detailed analysis of the strongly interacting Bose gas subject to a weak optical potential, leading us to a proper understanding of the new instability.

A weakly interacting atomic gas subject to an optical lattice is well described by the Bose-Hubbard model, which starts from a tight-binding model and takes the interaction between bosons into account perturbatively; the hopping amplitude $`J(V)`$ and the on-site interaction energy $`U(V,\gamma )`$ follow from the underlying parameters of the atomic gas, the dimensionless interaction parameter $`\gamma `$ and the amplitude $`V`$ of the optical lattice. The phase diagram of the Bose-Hubbard model is well known Fisher et al. (1989) and involves insulating Mott-lobes embedded in a superfluid phase. In 3D, the mean-field analysis for densities commensurate with the lattice provides the critical parameter $`U/J|_{\mathrm{S}\mathrm{MI}}5.8z`$, in good agreement with the experimental findings of Greiner et al. Greiner et al. (2002) (here, $`z`$ denotes the number of nearest neighbors). Going to 1D, fluctuations become important and appreciably modify the mean-field result: numerical simulations Kühner and Monien (1998); Rapsch et al. (1999) place the transition at the critical value $`U/J|_{\mathrm{S}\mathrm{MI}}=2C3.84`$. This result is easily transformed into the $`\gamma `$$`V`$ phase diagram of the weakly interacting atomic gas, once the relations $`J(V)`$ and $`U(V,\gamma )`$ to the Mott-Hubbard parameters are known; combining $`J\mathrm{exp}[2(V/E_r)^{1/2}]`$ as derived from a WKB-calculation, with the on-site repulsion $`U\gamma `$, the transition line $`V_c(\gamma )`$ (see Fig. 1) then derives from the implicit equation

$$4V/E_r=\mathrm{ln}^2\left[4\sqrt{2}\pi C\left(V/E_r\right)^{1/2}/\gamma \right].$$
(1)
Here, $`E_r=\mathrm{}^2k^2/2m`$ is the recoil energy with $`m`$ the boson mass and $`k=2\pi /\lambda `$ the wave vector of the light generating the optical lattice $`V(x)=V\mathrm{sin}^2(kx)`$. The dimensionless interaction strength $`\gamma `$ is defined via $`\gamma =mg/\mathrm{}^2n`$, with $`n`$ the density and $`g`$ the strength of the $`\delta `$-function interaction potential; $`g`$ is related to the 3D scattering length $`a_s`$ and the transverse confining frequency $`\omega _{}`$ via $`g=2\mathrm{}\omega _{}a_s`$ Olshanii (1998).

Increasing the interaction strength $`\gamma `$, the critical amplitude $`V_c`$ of the optical lattice triggering the S–MI transition decreases, see Fig. 1; the description of the atom gas in terms of the Bose-Hubbard model breaks down and we have to look for a new starting point. For a weak optical potential, a natural choice is the 1D Bose gas with $`\delta `$-function interaction, which resides in the strong coupling regime at small densities $`n`$ Petrov et al. (2000); the presence of the optical lattice is taken into account perturbatively. The homogeneous 1D Bose gas with $`\delta `$-function interaction has been solved exactly by Lieb and Liniger Lieb and Liniger (1963); the corresponding low-energy physics is properly described in terms of a Luttinger liquid with a parameter $`K(\gamma )`$ derived from the exact solution, see below. Adding the optical potential $`V`$, we arrive at the sine-Gordon model; the critical value $`K(\gamma ,V)=2(1+V/4E_r)`$ separating the Mott insulating phase from the superfluid one determines the phase line $`V_c/E_r=(\gamma ^1\gamma _c^1)/5.5`$, with the critical value $`\gamma _c3.5`$. Combining the results of the Bose-Hubbard- and sine-Gordon models we can complete the phase diagram for the commensurable situation as shown in Fig. 1: most remarkable is the appearance of a critical interaction strength $`\gamma _c`$ above which an arbitrary weak optical lattice is able to pin the system into a Mott insulator state. The presence of this instability is due to the closeness of the dilute 1D Bose liquid to Wigner-crystallization. Tuning the system away from commensurability with $`Q2\pi (n2/\lambda )0`$ the Mott insulator survives up to a critical misfit $`Q_c(V,\gamma )`$, see Fig. 1. The corresponding physics is similar to that of the commensurate-incommensurate transition of adsorbates on a periodic substrate as studied by Pokrovsky and Talapov Pokrovsky and Talapov (1979).

For a qualitative understanding of the transition in the strongly interacting 1D Bose gas, it is useful to consider the limit $`\gamma 1`$, where the behavior is essentially that of an ideal Fermi gas Girardeau (1960). A weak periodic potential $`V(x)=V\mathrm{sin}^2(kx)`$, with a lattice constant such that an integer number $`i=1,2,\mathrm{}`$ of particles will fit into one unit cell, gives rise to a single-particle band structure in which the $`i`$ lowest bands are completely filled. The ground state for noninteracting Fermions is then a trivial band insulator, separated from excited states by an energy gap which scales like $`(V/E_r)^i`$ in the limit $`VE_r`$. Similar to the Mott phase in the Bose-Hubbard model, the insulating state has a fixed integer density, commensurate with the lattice. It remains locked in a finite regime of the chemical potential, characterizing an incompressible state Fisher et al. (1989). Clearly, for a weak periodic potential, the lowest energy gap with $`i=1`$ is much larger than the higher-order ones. In the following, we will thus confine ourselves to studying the commensurate-incommensurate transition near integer filling $`i=1`$, where the commensurate phase has maximal stability.

For a quantitative theoretical analysis, it is convenient to use Haldane’s description of the interacting 1D Bose gas in terms of its long wave length density oscillations Haldane (1981). Introducing the two fields $`\varphi (x)`$ and $`\theta (x)`$ describing phase and density fluctuations, the Hamiltonian of the homogeneous gas (without longitudinal confining trap and periodic potential) is a quadratic form involving kinetic and interaction energies,

$$H_0=\frac{\mathrm{}}{2\pi }𝑑x\left[v_J(_x\varphi )^2+v_N(_x\theta \pi n)^2\right].$$
(2)
Here, $`v_J=\pi \mathrm{}n/m`$ is the analog of the bare ‘Fermi’ velocity at given average density $`n`$ while $`v_N=_n\mu /\pi \mathrm{}`$ is determined by the inverse compressibility, giving a sound velocity $`v_s=\sqrt{v_Jv_N}`$ consistent with the standard thermodynamic relation $`mv_s^2=n_n\mu `$. For short-range repulsive interactions, the dimensionless ratio $`K=\beta ^2/4\pi =v_J/v_s`$ is larger than one, approaching $`\mathrm{}`$ in the noninteracting case and unity in the hard core limit Haldane (1981). In the idealized model with $`\delta `$-function interactions, the exact solution by Lieb and Liniger shows that $`K`$ is a monotonically decreasing function of the ratio $`\gamma =mg/\mathrm{}^2n`$ between the interaction and kinetic energies; the limiting behavior for small values of $`\gamma `$

$$K(\gamma 0)=\pi [\gamma (1/2\pi )\gamma ^{3/2}]^{1/2}$$
(3)
follows from the Bogoliubov approximation in 1D. Surprisingly, this result remains quantitatively correct for $`\gamma `$ values up to $`10`$ Lieb and Liniger (1963). At large $`\gamma >10`$, the asymptotic behavior is $`K(\gamma \mathrm{})=(1+2/\gamma )^2`$.

For our subsequent analysis, it is convenient to introduce the conjugate fields $`\mathrm{\Theta }(x)=(2/\beta )\left[\theta (x)\pi nx\right]`$ and $`\mathrm{\Pi }(x)=\mathrm{}(\beta /2\pi )_x\varphi (x)`$ obeying the commutation relation $`[\mathrm{\Theta }(x),\mathrm{\Pi }(x^{})]=i\mathrm{}\delta (xx^{})`$. Using these fields,

$$H_0=\frac{\mathrm{}v_s}{2}𝑑x\left[\left(\mathrm{\Pi }/\mathrm{}\right)^2+\left(_x\mathrm{\Theta }\right)^2\right]$$
(4)
takes the form of the Hamiltonian for a 1D harmonic string with a linear spectrum $`\omega =v_sq`$ describing the long wave length density modulations of the interacting Bose gas. The assumption of a linear spectrum inherent in the simple form (4) is valid only below a momentum cutoff $`1/a\pi n`$ Lieb and Liniger (1963); the choice of the length scale $`a`$ fixes the energy scale of $`H_0`$. The Boson density operator $`n(x)`$ is related to the field $`\mathrm{\Theta }(x)`$ by Haldane (1981)

$$n(x)=\left[n+\frac{\beta }{2\pi }_x\mathrm{\Theta }\right]\left[1+2\underset{l=1}{\overset{\mathrm{}}{}}\mathrm{cos}\left(\frac{l\beta \mathrm{\Theta }}{2}+l\pi nx\right)\right];$$
(5)
the last factor accounts for the discrete nature of the particles. Adding an external periodic potential with amplitude $`V/2`$ and period $`\lambda /2`$ gives rise to the perturbation

$$H_\mathrm{V}=\frac{V}{2}𝑑xn(x)\mathrm{cos}\frac{4\pi x}{\lambda }.$$
(6)
As noted already by Haldane, insertion of the Fourier expansion (5) generates terms of the type appearing in the quantum (1+1)-dimensional sine-Gordon theory Coleman (1975); Gogolin et al. (1998). Close to commensurability the dominant term arising from the lowest harmonic in (5) has the conventional sine-Gordon form Fisher et al. (1989)

$$H_\mathrm{V}=\frac{Vn}{2}𝑑x\mathrm{cos}\left[\beta \mathrm{\Theta }+Qx\right]$$
(7)
with coupling parameter $`\beta =2(\pi K)^{1/2}`$ and a twist $`Q=2\pi (n2/\lambda )`$. The strength of the nonlinear $`\mathrm{cos}\beta \mathrm{\Theta }`$\- perturbation is conveniently expressed through the dimensionless parameter $`u=\pi a^2nV/2\mathrm{}v_s`$ which naturally involves the cutoff parameter $`a`$ Kehrein (1999). The twist $`Q`$ vanishes at commensurability; away from commensurability the finite twist $`Q`$ acts as a chemical potential for excitations and is preferably incorporated into the free Hamiltonian (4) via the replacement $`_x\mathrm{\Theta }_x\mathrm{\Theta }Q/\beta `$.

At fixed potential $`V`$, the quantum sine-Gordon model describes the competition between the preferred average inter-particle distance at given density due to the repulsive interaction and the period imposed by the external potential. A perturbative calculation (see Ref. Gogolin et al., 1998 for a review) tells that for $`\beta ^2/4\pi =K>2`$ a weak periodic potential is unable to pin the density; hence for $`K>K_c=2`$ ($`\gamma <\gamma _c3.5`$) the ground state remains gapless and superfluid in the presence of a small-amplitude lattice potential. In the strong coupling regime $`K<2`$, however, the atoms are locked even to a weak lattice, as long as the twist $`Q`$ is less than a critical value $`Q_c`$. Beyond that, there is a finite density of ‘solitons’ (or domain walls for adsorbates on a periodic substrate, see Ref. Pokrovsky and Talapov, 1979), which interpolate between minima of the external potential, relieving the frustration present at incommensurate densities $`Q0`$. The ‘solitons’ behave like relativistic particles with energy $`E_q=\mathrm{}v_s\sqrt{q^2+M^2}`$ and reestablish the superfluid response. The ‘mass’ $`M`$ determines the excitation gap in the Mott insulating state, which translates into a jump $`\mathrm{\Delta }\mu `$ in the chemical potential at the commensurate density: given that an additional/missing atom involves $`K`$ solitonic excitations with energy $`E_{q=0}=\mathrm{}v_sM`$ Japaridze and Nersesyan (1984) one obtains

$$\mathrm{\Delta }\mu =\frac{2\pi \mathrm{}^2n}{m}M;$$
(8)
furthermore, this mass is also simply related to the critical twist via $`Q_c=2K^2M`$. The precise numerical value of $`M`$ depends on the high momentum cutoff $`1/a`$ via the dimensionless amplitude $`u`$ of the lattice potential. The free-fermion limit $`K=1`$ fixes this cutoff at $`1/a=\pi n`$, resulting in the simple form $`u=V/4E_r`$; we ignore small corrections arising due to a possible modification in the cutoff away from $`K=1`$. The dependence of the mass $`M`$ on $`u`$ can be obtained from a recent nonperturbative renormalization group analysis of the quantum sine-Gordon model by Kehrein Kehrein (1999); for small values $`VE_r`$ and $`K`$ away from $`K_c`$, the gap in the chemical potential takes the form

$$\mathrm{\Delta }\mu =2E_r\left[\frac{V}{(2K)4E_r}\right]^{1/(2K)}.$$
(9)
In the limit $`K1`$, the size of the gap approaches $`V/2`$, in agreement with the above fermionic picture of the Tonks gas limit, where the appearance of an insulating ground state, even for a weak periodic potential, is due to the opening of a single-particle band gap at the Fermi energy. In the practically accessible regime of large but finite $`\gamma >\gamma _c`$ the gap depends on $`V`$ in the more complicated manner as given by (9) and vanishes exponentially as $`K`$ approaches $`K_c=2`$ Kehrein (1999), see Fig. 2(a). Similarly, the density range $`n2/\lambda =\pm Q_c/2\pi `$ over which the ground state remains locked approaches zero as $`KK_c`$. Finally, for $`K>2`$, the dependence of the critical interaction parameter $`K_c`$ (or $`\gamma _c`$) on the lattice amplitude $`V`$ follows easily from the Kosterlitz-Thouless nature of the scaling flow near $`K_c=2`$: to lowest order in $`u`$, $`K_c(u)=2(1+u)`$. Combining this result with (3) it is straightforward to determine the line $`V_c(\gamma )`$ separating the gapped insulating regime from the superfluid at small but finite values of $`u`$, see Fig. 1.

In order to analyze the consequences of the commensurate-incommensurate transition for cold atoms in a trap, we consider a 1D Bose gas in the presence of a weak longitudinal confining potential $`V(x)=m\omega ^2x^2/2`$. Provided the associated oscillator length $`l=(\mathrm{}/m\omega )^{1/2}`$ is much larger than the inter-particle distance, the density profile in this inhomogeneous situation may be obtained from the Thomas-Fermi approximation Dunjko et al. (2001)

$$\mu \left[n(x)\right]+V(x)=\mu [n(0)]$$
(10)
where $`\mu [n]`$ is the chemical potential of the homogeneous system. The central density $`n(0)`$ and the associated radius $`R=\left(2\mu [n(0)]/m\omega ^2\right)^{1/2}`$ of the cloud is obtained from the normalization condition

$$N=_R^R𝑑xn(x)=2R_0^{n(0)}𝑑n\sqrt{1\frac{\mu [n]}{\mu [n(0)]}}$$
(11)
with $`N`$ the particle number, provided the relation $`\mu [n]`$ is known explicitly.

In the absence of an optical lattice the density profiles are smooth functions of the coupling $`\gamma `$ Dunjko et al. (2001). In the limit $`\gamma 1`$, $`\mu [n]\mu _\mathrm{F}[n]=(\mathrm{}\pi n)^2/2m`$ approaches the Fermi energy of an ideal Fermi gas with density $`n`$, resulting in a profile $`n(x)=(2N/\pi R)\sqrt{1(x/R)^2}`$ with radius $`R=(2N)^{1/2}l`$. Adding now a periodic potential which is nearly commensurate with the density in the trap center opens a gap $`\mathrm{\Delta }\mu `$ in the chemical potential. The density profile will then exhibit a flat, incompressible regime in the trap center. In order to find the detailed shape of the density profile, we approximate the sine-Gordon model in its strong coupling, gapped phase by a gas of noninteracting relativistic Fermions of mass $`M`$; in the relevant regime $`1<K<2`$, this approximation is known to work extremely well Kehrein (1999). The chemical potential as a function of density can then be calculated explicitly,

$$\mu [n]\mu _\mathrm{F}[n]+\frac{\mathrm{\Delta }\mu }{2}f\left(\frac{4KE_r}{\mathrm{\Delta }\mu }\frac{\delta n}{n_c}\right)$$
(12)
with $`\delta n=nn_c`$ ($`n_c=2/\lambda `$). The dimensionless function

$$f(z)=\pm \left(1+z^2\right)^{1/2}z$$
(13)
is discontinous at $`z=0`$ and incorporates, via (12), the jump $`\mathrm{\Delta }\mu `$ in the chemical potential at $`n_c`$ associated with the incompressible commensurate state. Using this approximation, the density profile and the fraction of particles participating in the commensurate phase derive from an integration of (11) with (12), see Fig. 2(b). Knowledge of the locked fraction $`\mathrm{\Delta }N/N`$ plays an important role in the experimental detection of a commensurate Mott phase. This can be achieved by measuring the excitation gap through a phase gradient method as done previously for the Bose-Hubbard transition Greiner et al. (2002). Alternatively, it should be possible to directly observe the increase in the long-range translational order in the Mott phase via Bragg diffraction Birkl et al. (1995); Weidemüller et al. (1995); in either case the fraction $`\mathrm{\Delta }N/N`$ determines the experimentally available signal. The latter can be further enhanced by generating an array of parallel ‘atom wires’ with the help of a strong 2D optical transverse lattice. Using numbers similar to those in the recent experiment by Greiner et al. Greiner et al. (2002), it is possible to generate several thousand parallel 1D wires with a transverse confining frequency $`\nu _{}=20`$ kHz. A longitudinal harmonic trap with frequency $`\nu =40`$ Hz then encloses $`N50`$ atoms in each 1D wire; the associated central density in the absence of a longitudinal periodic potential is $`n(0)=2\mu \text{m}^1`$ for $`\gamma 1`$, commensurate with the lattice constant $`\lambda /20.5\mu `$m of a typical optical lattice Greiner et al. (2002). A weak periodic potential will then lead to an incompressible Mott state in the center of the cloud, provided the parameter $`\gamma =2a_s/n(0)l_{}^2`$ is larger than the critical value $`\gamma _c3.5`$. For 87Rb with a scattering length $`a_s5`$ nm, the resulting $`\gamma `$ is equal to one, i.e., not quite in the required range. As noted already by Petrov et al. Petrov et al. (2000), however, larger and in particular tunable values of $`\gamma `$ may be realized by changing $`a_s`$ via a Feshbach resonance as present, e.g., in 85Rb.

In conclusion we have shown that a commensurate Mott state can be realized in dilute 1D BEC’s already with an arbitrary weak lattice potential, provided that the ratio $`\gamma `$ between the interaction and kinetic energies is larger than a critical value $`\gamma _c3.5`$. This instability provides a new and experimentally accessible tool for the quantitative characterization of 1D atomic gases in the strongly correlated ‘Tonks gas’ limit. Also, the observation of a Mott state in a regime where the atoms are not confined to discrete lattice sites would give direct evidence for the granularity of matter in strongly interacting dilute gases.

It is a pleasure to acknowledge fruitful discussions with I. Bloch, T. Esslinger, M. Greiner and S. Kehrein. This work was supported by the DFG Schwerpunkt Ultrakalte Quantengase.