text:
# ON THE CORRELATED X-RAY AND OPTICAL EVOLUTION OF SS CYGNI

## 1 INTRODUCTION

SS Cygni is a dwarf nova (DN) cataclysmic variable. The main characteristic of DN are outbursts which have recurrence timescales of days to years, and last for days to weeks. Outbursts occur in SS Cygni every $`40`$ d, when the source rises from a quiescent magnitude of $`V12`$ to an outburst magnitude of $`V8.5`$ Szkody & Mattei (1984). The outbursts of SS Cygni exhibit a bi-modal distribution, with wide and narrow outbursts that last $`20`$ and $`12`$ d, respectively Szkody & Mattei (1984).

SS Cygni was first detected in X-rays by Rappaport et al. (1974). The main site of X-ray emission in DN is believed to be the boundary layer between the white dwarf and the accretion disk. Accretion material in the innermost orbit of the disk is decelerated in the thin boundary layer before coming to rest on the surface of the white dwarf. During this process the material should release $``$ half of its accretion luminosity.

SS Cygni shows both hard and soft X-ray spectral components Córdova et al. (1980). The two spectral components are usually attributed to different states of the boundary layer flow, which depend on $`\dot{M}`$ Pringle & Savonije (1979). For $`\dot{M}\dot{M}_{\mathrm{crit}}`$, $`\dot{M}_{\mathrm{crit}}10^{16}`$ g s-1, the heated gas expands adiabatically and forms a hot ($`T10^8`$ K) corona which cools by thermal bremsstrahlung, radiating a hard spectrum. For $`\dot{M}\dot{M}_{\mathrm{crit}}`$ the heated gas cools without expanding out of the optically thick boundary layers, giving rise to an approximate blackbody spectrum (at $`T10^5`$ K). During this high-accretion state, tenuous regions of the boundary layer flow, above and below the center plane, remain in the hot, optically thin state, yielding a residual hard component.

Comparison of X-ray data from EXOSAT with the optical light curve of SS Cygni by Jones & Watson (1992) confirmed this X-ray behavior. In quiescence, the hard component is strong. The outburst is first seen in the optical band. After a short delay the hard X-ray component flares, then declines to a level below that of quiescence. The ultrasoft X-ray outburst begins as the hard X-ray flare ends. The ultrasoft X-ray emission tracks the optical emission closely during the peak of the optical burst, but the ultrasoft X-rays decline to quiescence more quickly. Once they start to decline, the hard X-ray component increases again, to a level brighter than quiescence, then declines more slowly than the optical in the final decline. Studies of detailed RXTE/PCA observations of single outbursts of SS Cygni exhibit this X-ray/optical behavior (Wheatley, Mauche & Mattei 2000; Priedhorsky, Trudolyubov & McGowan 2002; McGowan, Priedhorsky & Trudolyubov 2002; Wheatley, Mauche & Mattei 2003).

We analyze several outbursts of SS Cygni from archival data, to investigate the optical and X-ray evolution of the source. In particular, we want to confirm the optical/X-ray outburst phenomenology found in previous works, which were based on studies of only a few outbursts.

## 2 OBSERVATIONS AND DATA ANALYSIS

Its brightness and frequent outbursts makes SS Cygni an ideal candidate for observations by amateur astronomers. Since its discovery in 1896, SS Cygni has been monitored by the American Association of Variable Star Observers (AAVSO***http://www.aavso.org). We have obtained the 1996–2001 light curve of SS Cygni from the AAVSO.

We have also obtained the X-ray light curve of SS Cygni from the All Sky Monitor (ASM; Doty 1988; Levine et al. 1996) that is onboard the Rossi X-ray Timing Explorer (RXTE) satellitehttp://xte.mit.edu (Bradt, Rothschild & Swank 1993). The data contained in the archive includes background subtraction. The ASM has three X-ray channels spanning $`1.312.2`$ keVftp://legacy.gsfc.nasa.gov, from which we derive $`1.33.0`$ keV (band1) and $`3.012.2`$ keV (band2) fluxes for our analysis. While the ASM is small, with a $`90`$ cm2 effective area, it is persistent in its coverage.

We also use the much more sensitive, but sporadic, pointed observations taken by the Proportional Counter Array (PCA; Jahoda et al. 1996) instrument on RXTE in the $`320`$ keV band. For our analysis we used publicly available RXTE/PCA data collected during $`19962000`$ (a total of $`262`$, observations $`100020000`$ s duration). The observation log for the data used in our analysis is shown in Table 1.

The PCA data were reduced with the standard RXTE FTOOLS package, version 5.2. We used PCA data collected in the $`320`$ keV energy range for the spectral analysis. The PCA response matrices for individual observations were constructed using the pcarmf task and the background estimation was performed by applying very large events (VLE) or faint source models, depending on the level of the source flux §§§See http://heasarc.gsfc.nasa.gov/docs/xte/recipes/pcabackest.html. The standard dead-time correction procedure was applied to the PCA data. In order to account for the uncertainties of the response matrix, a $`0.5\%`$ systematic error was added in quadrature to the statistical error for each PCA energy channel.

We generated energy spectra of SS Cygni averaging the data over whole individual observations using the Standard 2 mode data (binned data; 129-channel spectra accumulated every 16 s). For the temporal analysis of SS Cygni we used the PCA data in Standard 2 and Good Xenon modes (unbinned data; time-stamped events with 256-channel resolution).

## 3 COMPARISON OF AAVSO DATA AND ASM DATA

The contemporaneous AAVSO and ASM data cover the period from 1996 January 4 to 2001 August 30. The data are shown in Fig. 1; while the optical outbursts are obvious, it is difficult to see individual outbursts in the X-rays.

In order to investigate the optical and X-ray variability we have adopted a stacking technique to average the optical and X-ray data and improve signal to noise. The analysis of the outburst and quiescent (inter-outburst) data are presented in the following sections.

### 3.1 Outburst

We separate the optical data into wide (Fig. 1, e.g. JD 2450320), narrow (Fig. 1, e.g. JD 2450370), and anomalous or mini outbursts (Fig. 1, e.g. JD 2451050). In our data we find 22 wide, 27 narrow and 7 anomalous outbursts. For rescaling purposes we define start ($`T_{start}`$) and end points ($`T_{end}`$) to occur at $`V=10.5`$ for the wide and narrow optical outbursts, and at $`V=11.7`$ for the anomalous outbursts. By calculating the average duration of the wide, narrow, and anomalous outbursts we determine a characteristic length $`\mathrm{\Delta }T`$ for the three different outbursts ($`\mathrm{\Delta }T=16`$ d for the wide, 9 d for the narrow and 3 d for the anomalous outbursts). To ensure we do not ignore any data which may be part of an outburst we include $`30`$ d prior to $`T_{start}`$ and $`40`$ d after $`T_{start}`$ in each cut of data for the wide and narrow outbursts, and $`5`$ d before and $`10`$ d after $`T_{start}`$ for the anomalous outbursts. Using the start and end times of the sections of optical data, we separate the band1 ($`1.33.0`$ keV) and band2 ($`3.012.2`$ keV) X-ray data. The optical and X-ray data are rescaled on the characteristic outburst timescales such that $`T_0=T_{start}`$ and $`T_{\mathrm{\Delta }T}=T_{end}`$. We can then calculate the average time history for each type of outburst, in the optical, band1 and band2 X-rays, by stacking and averaging the rescaled data (Fig. 2).

For both the wide and narrow outbursts (Fig. 2, left and middle panels), the band2 X-rays ($`312.2`$ keV) are anti-correlated with the optical flux, i.e. the band2 X-ray flux falls when the optical flux rises during outburst. This is consistent with previous reports that the hard X-rays are suppressed during optical outburst Jones & Watson (1992); Priedhorsky et al. (2002); McGowan et al. (2002); Wheatley et al. (2003). We have confirmed that this behavior holds over a large ensemble of outbursts. The previously-reported X-ray maximum at the decline of an outburst is also visible. The band2 X-rays for the narrow outburst shows an increase in flux which could be associated with the previously-reported X-ray maximum at the beginning of the optical outburst (Fig. 2, middle panel). However, as this is based on only one point, and a similar behavior is not seen in the wide outburst, this result is inconclusive. The averaged anomalous outbursts (Fig. 2, right panel) suggest a positive correlation between the optical and band2 X-ray flux.

### 3.2 Quiescence

We studied the optical and X-ray correlation during quiescence using a similar method as for the outburst data. We use the start and end times of the optical outbursts determined above to separate the inter-outburst data. The inter-outburst was therefore the interval between the $`V=10.5`$ crossing on the decline from outburst, and the $`V=10.5`$ crossing on the rise to outburst. We note that this is an arbitrary definition of quiescence, which we use only for our stacking analysis. We removed all data from times of anomalous outburst. We calculate the average duration of quiescence, and define a characteristic length for the quiescent interval of 38 d which we used to rescale the data. The optical, band1 and band2 X-ray data are then averaged by stacking all the quiescent intervals (Fig. 3).

A linear fit including a decline to the inter-outburst $`1.312.2`$ keV X-ray flux of SS Cygni results in a $`\chi ^2`$ of 36.8 for 29 d.o.f. We chose the interval from rescaled days 5 to 36 as the largest span that is clearly free from the optical outburst. During this interval the flux seems to fade, with a decline significant at the 3.6-sigma level. The error bars on the ASM X-ray fluxes reflect all currently known systematic effects. The best linear fit to the $`1.312.2`$ keV X-ray flux between rescaled days 5 and 36 declines at the rate of 1.3% per scaled day, with a decrease of $`40`$% over the day 5 – day 36 interval (see Fig. 3, fourth panel). In other words, the flux between days 5 and 36 was fit by a mean flux of 0.354 counts s-1 and a slope of $`5.7\times 10^3`$ counts s-1 day-1, with an error in the slope of $`\pm 1.6\times 10^3`$ counts s-1 day-1, based on an increase in the $`\chi ^2`$ of 1 for one interesting parameter.

#### 3.2.1 Folding on the Orbital Period

To test for a modulation at the orbital period, we combined the band1 and band2 X-ray data from the same quiescent intervals. We folded the $`1.312.2`$ keV data on $`P_{orb}=0.2751297`$ d using $`T_0=`$ HJD 2447403.6295, where $`T_0`$ corresponds to the time of inferior conjunction, i.e. the time of zero velocity of the absorption line star Friend et al. (1990). The resulting light curve is shown in Fig. 4.

We fit the data with a constant linear fit and a sinusoid. The best linear fit results in $`\chi ^2=7.1`$ for 19 d.o.f., the best-fit sinusoid gives $`\chi ^2=5.4`$ for 16 d.o.f. The semi-amplitude of the sinusoid is 0.024 counts which corresponds to a $`7`$% modulation. While the sinusoid gives a better fit to the data, employing the $`F`$-test we find a probability of 78%. This indicates that the variation is random, thereby failing to rule out a constant flux.

## 4 COMPARISON OF SHORT-TERM OPTICAL AND X-RAY TIME HISTORIES

### 4.1 Light curve phenomenology. Correlated hard X-ray and optical evolution.

In five years, six optical outbursts of SS Cygni were covered with RXTE/PCA pointed observations. Two of the outbursts were poorly sampled, therefore we present results from only four of the outbursts here. In spite of the differences in the duration of the four outbursts, the morphology of the hard X-ray evolution ($`320`$ keV) is essentially the same in each (Fig. 5, upper and middle panels). Unlike the soft X-ray emission below $`1`$ keV Jones & Watson (1992), the hard X-ray emission of SS Cygni shows a general anticorrelation with optical flux, except for “brief” X-ray/optical outbursts (’spikes’). The positions of the spikes, which last for $`<2`$ d, are marked in Fig. 5 (right panel) at MJDs 51610, 51622, 51644 and 51680. The rise of the optical outburst coincides with a short spike in the hard X-rays: the flux increases by several times, then fades, all in $`<1`$ d. At the maximum of the optical outburst, the hard X-ray flux is at its lowest. The decline of the optical outburst corresponds to a gradual increase of the hard X-ray flux, followed by a brief $`24`$ d maximum. Periods of optical quiescence are generally characterized by a relatively high level of hard X-ray flux. Rare anomalous optical outbursts are accompanied by short X-ray spikes (i.e. MJDs 51610 and 51644 in Fig. 5, right panel).

We studied the characteristic timescales of the hard X-ray luminosity transitions. The average duration of the X-ray spike coincident with an anomalous optical outburst is $`<2`$ d, with a nearly symmetrical time history. Exponential fits to the rising/decaying parts of the X-ray spikes (MJDs 51610 and 51644) give an e-folding time of $`0.5`$ d (Fig. 5).

The X-ray spikes corresponding to the beginning of the optical outburst have a shorter duration of $`<1`$ d. The observations of 1996 October 9–11 allow a detailed study of the source evolution during such a spike (Fig. 5, left panel). The light curve of SS Cygni during the spike, based on four consecutive PCA observations Wheatley et al. (2000, 2003), is shown in Fig. 6.

The spike has an asymmetric shape with an initial quasi-exponential rise with $`\tau _{rise}25000`$ s followed by an abrupt decline in $`5000`$ s. The total duration of the hard X-ray spike is about half a day. The source $`320`$ keV flux rises by $`5`$ times from the quiescent level to the spike maximum and then drops by a factor of $`18`$ during the decline (Fig. 5, left panel). According to the AAVSO data, the optical outburst starts $`1`$ day prior to the hard X-ray spike (Wheatley et al. 2000), while the UV rise is coincident with sharp decline in hard X-rays.

The hard X-ray recovery at the end of an outburst has a consistent pattern from one outburst to another (Fig. 5). An exponential fit to the optical decline and X-ray rise is shown in Table 2. The X-rays appear to evolve almost twice as quickly as the optical light.

### 4.2 Energy spectra

Previous X-ray studies of SS Cygni Swank (1979); Córdova et al. (1980); Jones & Watson (1992); Yoshida et al. (1992) have demonstrated the presence of two distinct spectral components: a hard component with a cut-off at $`1020`$ keV (bremsstrahlung model), and an ultrasoft optically thick component that is cut-off at $`60`$ eV. These components behave very differently during optical outburst and quiescence. The energy range of the PCA detector allows the direct study of the hard spectral component without contamination from an ultrasoft component. The ultrasoft component is invisible to the ASM, also.

The RXTE/PCA spectral data were approximated with one of the simplest models that fit the data: a two-component XSPEC (version 11) model that includes an absorbed bremsstrahlung emission model and a Gaussian emission line at $`6.7`$ keV. Line emission from the source can be present at 6.4, 6.7 and/or 6.9 keV, however the energy resolution of the PCA is not sensitive enough to be able to distinguish between these lines, and we are most likely seeing a blend. While this two-component model is unphysical, it allows us to to compare our spectral results directly with previously published results. The width of the Gaussian line was fixed at 0.1 keV due to the limited energy resolution of PCA detector ($`1`$ keV at 6 keV), while the line centroid energy was left as a free parameter. Typical energy spectra of SS Cygni corresponding to the optical quiescence (upper histogram) and outburst (lower histogram) are shown in Fig. 7.

The combination of an absorbed bremsstrahlung radiation model and a Gaussian emission line at $`6.7`$ keV provides a satisfactory description for the PCA spectral data in the $`320`$ keV energy band (Fig. 7). The best-fit value of the characteristic temperature, $`\mathrm{kT}_{\mathrm{bremss}}`$, of the bremsstrahlung model varies between $`6`$ and $`26`$ keV, and the resulting model $`320`$ keV energy flux changes from $`10^{11}`$ to $`7\times 10^{10}`$ ergs s-1 cm-2. The evolution of best-fit spectral parameters (X-ray flux and $`\mathrm{kT}_{\mathrm{bremss}}`$) during the 1996, 1999 and 2000 RXTE/PCA observations is shown in Fig. 5 (middle and lower panels). There is a unique correlation between both the hard X-ray flux and hardness of the X-ray spectrum, and the optical flux level. The X-ray flux is high and the X-ray spectrum is hard during optical quiescence, and at short periods at the beginning and the end of optical outburst. On the other hand, during the maximum of optical outburst the X-ray flux is extremely low and the X-ray spectrum is soft (corresponding to the lowest level of cut-off energy). The relative strength of the 6.7 keV emission line also is correlated with optical flux: the equivalent width of the emission line rises from $`350`$ eV during optical quiescence to $`800`$ eV during the maximum of the optical outburst.

The relation between the level of X-ray flux and the best-fit value of $`\mathrm{kT}_{\mathrm{bremss}}`$ (hardness-intensity diagram) is shown in Fig. 8. Each point represents the data averaged over one individual observation. Fig. 8 clearly shows that the observations can be subdivided into three distinct regions according to their positions on the hardness-intensity diagram, illustrating the correlation between X-ray and optical evolution.

The first, densely populated region is concentrated around a temperature of $`20`$ keV and an X-ray flux of $`2\times 10^{10}`$ ergs s-1 cm-2. This region corresponds to the periods of optical quiescence, with a relatively high level of hard X-ray flux (Fig. 5). The second, densely populated region is concentrated around a temperature of $`8`$ keV and an X-ray flux of $`3\times 10^{11}`$ ergs s-1 cm-2. This group represents the periods of optical/soft X-ray outburst, characterized by an extremely low level of hard X-ray flux (Fig. 5). The third, sparsely populated group comes mostly from periods of transition between the optical outburst and quiescence, characterized by the rise of hard X-ray luminosity and general hardening of X-ray spectrum (Fig. 5). A small fraction of the observations from the third group correspond to the fast transitions from optical quiescence to the outburst.

## 5 DISCUSSION

It is generally believed that outbursts in DN are caused by a sudden brightening of the optically thick accretion disk due to an increased accretion rate (see Osaki 1996 for review). Two models were proposed to account for the intermittent accretion, the mass-transfer instability model (MTI; Bath 1973) and the disk instability model (DIM; Osaki 1974). It is now generally accepted that the DIM is a strong factor for DN outbursts Osaki (1996). In the DIM the secondary star’s mass transfer rate is continuous, but during quiescence the mass flow from the secondary is not fully accreted onto the white dwarf but is mostly held up in the outer parts of the accretion disk. Once the stored mass reaches a critical amount, a thermal instability Meyer & Meyer-Hofmeister (1981) within the disk causes this matter to flow through the disk and onto the white dwarf, resulting in optical and UV/soft X-ray outburst.

The general pattern of dwarf novae optical and X-ray behavior during outburst is traditionally explained in the framework of the DIM Lasota (2001). Two observed features, however, appear to pose significant problems for the simple single component disk models Meyer & Meyer-Hofmeister (1994) which are unable to reproduce: i) a sizable lag in rise-time of the UV light with respect to optical and hard X-rays rise at the beginning of the outburst and ii) the appearance and gradual decline of hard X-rays and UV light when the system falls back into quiescence.

The observed optical/hard X-ray/UV delay and the post-outburst evolution of the dwarf novae can be explained by considering models in which the inner edge of the accretion disk does not extend down to the surface of the white dwarf. Truncation of the disk can be due to evaporation of the inner part of the accretion disk Meyer & Meyer-Hofmeister (1994), magnetic fields (e.g. Livio & Pringle 1992), or to an advection-dominated accretion flow (ADAF; Menou 2001). Here we examine the observed multiwavelength behavior of SS Cygni in the framework of the first model, evaporation of the inner disk. Let us assume that the quiescent accretion flow in SS Cygni consists of an outer optically thick cool disk truncated by evaporation at a certain radius and, inside that, an optically thin quasi-spherical boundary layer/corona structure extending down to the surface of the white dwarf.

The relative timing of the optical, hard X-ray, and soft X-ray/UV components in an SS Cygni outburst are now well established Wheatley et al. (2003). We suggest that this sequence is an important clue for understanding the movement of material through the disk. To summarize the outburst phenomenology, after the initial optical decline, there is a short delay ($`0.5`$ days) before the flaring of the hard X-ray component and an additional $`1`$ day delay before the flare in the soft X-ray and UV bands Swank (1979); Wheatley et al. (2000, 2003). The key features of the hard X-ray flare are the fast quasi-exponential rise of the flux ($`\tau _{rise}25000`$ s) and an abrupt decline (Fig. 6) coincident with the rise of the UV flux (Wheatley et al. 2000).

These delays are conventionally interpreted as that passage of a heating wave through the accretion disk (Wheatley et al. 2003 and references therein). The wave starts in the outer part of the disk, causing the optical rise, then reaches the innermost region, where it initiates an increase in the mass flow from the inner edge of the disk to the boundary layer and white dwarf. This gradual increase in flow first yields an increase in the optically thin, hard X-ray emission; at a certain critical point, the flow switches its mode into an optically thick, soft X-ray emitting plasma.

We propose an alternative explanation for the optical/hard X-ray/ soft X-ray sequence, inspired by the fact that the optical/hard X-ray delay is comparable to the dynamical timescale to fall in from the outer region of the accretion disk. In our picture, when the outburst starts, the outer disk transitions from a cool to a hot state. The increased mass accretion rate gives rise to the optical flare. The increase of this mass transfer rate, either in the outer part of the accretion disk, or in the vicinity of the inner Lagrangian point, somehow leads to an increase of the mass flow through the optically thin inner region of the accretion flow on a dynamical timescale that is comparable to the orbital period of the SS Cygni system. This requires a path for the accreting material that bypasses diffusive flow through the disk, perhaps via the sub-Keplerian flow above the disk. Once this material arrives at the inner region, it causes the initial fast rise of the hard X-ray flux (Fig. 6). The accretion rate through the optically thick accretion disk increases on a diffusive timescale, causing the rise of UV and soft X-ray flux. At some point, the rise of the UV/soft X-ray flux, combined with the increase of the accretion rate through the optically thin region, will lead to its catastrophic cooling and an abrupt softening of X-ray spectrum and abrupt drop of hard X-ray flux, as seen in the RXTE light curve (Fig. 6).

In this alternative picture, the optically thick disk/boundary layer will extend down to the white dwarf surface when the outburst is at maximum, and radiation driven disk winds will probably strip the coronal layers from the disk Meyer & Meyer-Hofmeister (1994). The decline of the accretion rate during late stages of the outburst allows disk evaporation to resume, restoring the optically thin disk corona. As the corona returns, the hard X-ray flux recovers and the X-ray spectrum hardens. At the end of the optical outburst we expect that the hard X-ray flux should decay more slowly than the optical, as we observe, because the variation in $`\dot{M}`$ will be “smeared out” as it travels through the disk Bath & Pringle (1981); Jones & Watson (1992). The hot/cool transition seen directly in the optical is thereby delayed and spread out as it reaches the innermost disk.

This picture also explains the phenomenon of short hard X-ray spikes associated with optical anomalous outbursts. The short-term increase of the mass transfer rate in the outer optically thick disk or in the vicinity of the inner Lagrangian point creates an increase of the mass flow in the inner optically thin region, by the same bypass that takes place in a normal outburst, and causes an increase in the hard X-ray flux. However, the smaller short spikes decays as the small mass spike is accreted, rather than being quenched by the optically thick flow through the diffusive disk as in a normal outburst.

The series of pointed RXTE/PCA observations allow us to trace the evolution of the X-ray spectrum of SS Cygni in unprecedented detail. The results of spectral analysis of RXTE/PCA observations of SS Cygni provide additional support for the two-component model of the accretion flow mentioned above. The hard X-ray emission from non-magnetic CVs can be interpreted as due to: i) optically thin radiation from strong shocks formed in the boundary layer when the mass transfer rate was low; ii) an X-ray emitting corona-like structure above the plane of the disk/boundary layer (Meyer & Meyer-Hofmeister 1994; Frank, King & Raine 2002). In quiescence and at low disk accretion rates, the optically thin boundary layer gas is likely to be thermally unstable to turbulent viscous heating, rapidly achieving hard X-ray temperatures ($`10^8`$ K). At the same time, this gas is likely to expand out of the disk and boundary layer, and form a hot corona-like structure around most of the white dwarf. At quiescence (low $`\dot{M}`$) SS Cygni shows relatively bright, highly variable hard X-ray emission (see Wheatley et al. 2003). The energy spectrum of the source is extremely hard with a cut-off at $`2025`$ keV (corresponding to the “QUIESCENCE” and upper parts of the “TRANSITION” region in Fig. 8) indicating the temperature of the optically thin emitting gas is $`10^7`$ K. At high disk accretion rates (the maximum of the optical outburst) the instability will be largely suppressed, and the boundary layer luminosity mostly emitted as soft X-rays, while the corona-like structure will be completely cooled down by rising flux of soft photons from the disk/boundary layer. After this stage we can expect the radiation driven winds to remove the coronal layer, and the hard X-rays observed during optical outburst will originate from shocks connected with such winds Lucy & Solomon (1970); Mauche & Raymond (1987). The X-ray spectrum of SS Cygni during the optical outburst is much softer than in quiescence, with a cut-off energy of $`68`$ keV (Fig. 5, 8). The abrupt softening of the X-ray spectrum observed with RXTE/PCA at the end of the hard X-ray spikes (Fig. 5, 6) also supports the above interpretation. The recovery of the hard X-ray flux during the decline of the optical outburst is accompanied by gradual hardening of the source energy spectrum (Fig. 5, 8 (the “TRANSITION” region)), as expected from the recovery of the inner hot optically thin region. The subsequent decline of the optical and hard X-ray luminosity can be attributed to the further decrease of the mass accretion flow through the optically thick and thin regions of the accretion disk.

Our analysis of the RXTE/ASM and AAVSO observations of SS Cygni over a period of 6 years confirms the general picture of hard X-ray evolution deduced from more sensitive pointed observations of a handful of individual outbursts Swank (1979); Córdova et al. (1980); Jones & Watson (1992); Yoshida et al. (1992); Wheatley et al. (2000, 2003). The average ASM light curves of SS Cygni obtained by stacking a large number of individual outbursts still show the suppression of the hard ($`312.2`$ keV) X-rays during both wide and narrow optical outbursts (Fig. 2).

The RXTE/ASM observations of SS Cygni during quiescence show some evidence for a decrease in the hard X-ray flux during the inter-outburst periods. The observations of the post-outburst evolution of other dwarf novae often show a decrease of the optical and UV flux between outbursts Smak (2000). The observed decrease in the X-ray, UV and optical poses additional problems for the simple DIM models, which imply a gradual increase of the luminosity (in all bands) during quiescence. On the other hand, the two-component model of the accretion flow involving disk evaporation Meyer & Meyer-Hofmeister (1994) naturally solves this problem. As the outburst approaches its end, the inner edge of the optically thick accretion disk progressively moves outwards with time. The accretion through the coronal flow thereby decreases, causing a decline in the hard X-ray flux Meyer & Meyer-Hofmeister (1994). This decrease in X-rays reduces the level of X-ray flux illuminating the outer parts of the accretion disk. Hence, less X-rays are reprocessed which leads to a decrease in the optical flux.

## 6 CONCLUSIONS

We have analyzed the variability and spectral evolution of the prototype dwarf nova system SS Cygni using simultaneous RXTE data (hard X-rays) and AAVSO (optical) observations. The series of pointed RXTE/PCA observations allow us to trace the evolution of the X-ray properties of SS Cygni in unprecedented detail. We find that our results are in general agreement with predictions of modified disk instability models that imply a two-component accretion flow, formed through the evaporation of the inner part of the optically thick accretion disk Meyer & Meyer-Hofmeister (1994) or other process. We find that the optical and X-ray outburst morphology of SS Cygni predicted by the modified DIM and observed in a series of pointed observations is confirmed by stacking the marginal detections of individual outbursts from RXTE/ASM observations. The averaged light curves for the wide and narrow outbursts display a suppression of the hard ($`312.2`$ keV) X-ray flux during optical outbursts.

Based on the hard X-ray (RXTE/PCA/ASM) data and optical (AAVSO) observations, we quantify the time lags between the optical and X-ray outbursts, and the time scales of the X-ray spectral recovery from outburst. We also find that anomalous optical outbursts of SS Cygni are accompanied by short hard X-ray outbursts lasting $`12`$ days.

Using the RXTE/ASM data, we find evidence for a decrease in the hard X-ray flux during quiescence. While this is in disagreement with a simple DIM which predicts an increase in the disk mass, constant or decreasing flux has been observed in optical and UV studies of other dwarf novae (e.g. Smak 2000). It should be mentioned that a modified DIM involving a two-component accretion flow Meyer & Meyer-Hofmeister (1994) naturally solves this problem, by predicting the decrease of hard X-ray and optical/UV flux due to the outward motion of the inner boundary of optically thick accretion disk due to its evaporation into a corona-like structure.

We folded the quiescent $`1.312.2`$ keV X-ray data on the orbital period, but do not find any significant modulation. An all-sky monitor a few times more sensitive than the RXTE/ASM would allow us to study individual outbursts and the quiescent behavior of SS Cygni with much greater signal to noise, and to extend these studies to a large number of dwarf nova (Priedhorsky, Peele & Nugent 1996).

## 7 ACKNOWLEDGMENTS

We thank Peter Wheatley for useful discussions. This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center. In this research, we have used, and acknowledge with thanks, data from the AAVSO International Database, based on observations submitted to the AAVSO by variable star observers worldwide.


format: text

text:
# A new model for simulating colloidal dynamics

## 1 Introduction

Understanding the dynamics of colloidal dispersions has been a long-standing problem in condensed matter physics. Despite the continuous need for accurate theoretical predictions concerning particle mobilities or sedimentation rates, the development of the field meets a nearly unpassable barrier. The many-body character of the hydrodynamic interactions (HI) between the colloidal particles (i. e. their correlated motion as a result of solvent-mediated fast momentum transport) poses a serious challenge for analytic theory. On the other hand, the rapid growth of available computer power in the last decades has made it more and more feasible to descend to a more basic description of colloidal systems rather than trying to develop an adequate macroscopic description. Such a step was made recently for the statics of charged colloids, where simulations of the primitive electrolyte model, taking the particle character of the charges and fluctuations in the charge distribution fully into acount, resolved some important strong correlation issues and thus have resulted in a significant refinement of the existing theories of screening . In a similar fashion, many-body effects in the HI of suspensions of only moderate density are quite important. Although their analytical form is, in principle, known in terms of a multipole expansion , it is numerically very cumbersome to take these higher-order many-body terms into account in a Brownian or Stokesian dynamics  simulation, where only the positions of the colloidal particles enter. For this reason, it is more convenient and for large particle numbers also more efficient, to take these effects into account by simulating the solvent degrees of freedom, and, in particular, the momentum transport through the solvent explicitly. An additional bonus is that retardation effects, if present, are automatically taken into account.

Standard Molecular Dynamics (MD), where just Newton’s equations of motion are integrated numerically, is of course able to simulate HI correctly. However, it is very time-consuming to simulate the motion of a huge number of particles explicitly in such detail. The time step is governed by the local oscillations of the solvent particles in their temporary “cages”, which is much faster than the motion of the colloidal particles. This information is however not needed for correctly reproducing HI. One just needs an appropriate mechanism for momentum transfer on somewhat larger time scales (however still small compared to the motion of the colloidal particles). This can be done in various ways: Dissipative Particle Dynamics (DPD)  is essentially MD, where however the particles are made very soft in order to increase the time step, and where a momentum-conserving Langevin thermostat is added . Another possibility is multi-particle collision dynamics , where the solvent is modeled as an ideal gas, and collisions are modeled as simple stochastic updating rules which conserve energy and momentum. Grid-based methods can be either the direct solution of the Navier-Stokes equation by a finite difference method, or the Lattice Boltzmann (LB)  approach, where essentially a linearized Boltzmann equation is solved in a fully discretized version (i. e. space, time, and velocities are discretized).

In our opinion, all of these approaches have their advantages and disadvantages. Nevertheless, it seems that all of them are roughly equivalent with respect to their computational efficiency, and with respect to the level of accuracy of describing the physical phenomena. Therefore, we believe that the choice of method is mainly a matter of convenience.

For our model we have chosen the LB route, mainly because this algorithm is particularly simple to implement and parallelize. Moreover, we were inspired by previous highly successful simulations of hard-sphere colloidal suspensions . A further advantage of a grid-based method is that thermal fluctuations can be both turned on and off, depending on the physical situation under consideration (in a particle method they are always present). This approach is essentially a hybrid method, where the solvent is modeled via LB, while the colloidal particles are run with MD.

The original approach by Ladd  models the colloidal particles as extended hollow spheres, while stick boundary conditions at the surface are implemented (roughly spoken) via bounce-back collision rules. One should note that for moving spheres this involves some minor fluctuations in the fluid mass, which is included within a sphere: The moving shell incorporates some new fluid at the front, while it releases mass at the rear. The fluctuations are a natural effect of the thermal density fluctuations of the fluid, and of the lattice discretization. Using this method, the translational and the rotational dynamics of colloidal spheres were accurately described.

This model has been recently extended to the case of charged systems , where the colloidal particles carry a charge, while the counterions and salt ions are taken into account as LB populations, such that the electrostatics is essentially taken into account on the level of the Poisson-Boltzmann equation. This method has two disadvantages: firstly, the discrete nature of the ions and correlations beyond the Poisson-Boltzmann level are not taken into account; secondly, one cannot avoid a leakage of charge into (and out of) the sphere (as in the case of mass), such that it is hard (if not impossible) to maintain a well-defined Debye layer of charges around it . For these reasons, it is advisable to take the counterions explicitly into account. The disadvantage, however, is that only rather small size ratios (size of colloidal particle vs. size of ion) are easily accessible. Nevertheless, we believe it is useful to study such a system, with respect to both its equilibrium and its nonequilibrium properties. The purpose of the present paper is to describe first steps in the development and validation of the corresponding model.

For the coupling of the small ions to the LB hydrodynamics, one would like to simply model the former as point particles. Fortunately, such a coupling has been recently developed by Ahlrichs and Dünweg  with the purpose of studying the dynamics of polymer solutions . Each Brownian point particle is assigned a phenomenological friction coefficient $`\zeta `$, and the coupling is implemented just as a dissipative force $`\stackrel{}{F}=\zeta (\stackrel{}{V}\stackrel{}{u})`$ acting on the particle, where $`\stackrel{}{V}`$ is the particle velocity, while $`\stackrel{}{u}`$ is the solvent flow velocity at the particle’s position, obtained via linear interpolation from the surrounding lattice sites. Furthermore, the LB variables at these sites are adjusted to ensure momentum balance, and thermal fluctuations are added to both types of degrees of freedom. For more details on implementation and validation of this method, see Refs. .

It is then convenient to use this coupling not only for the ions, but also for the colloidal particles, whose larger size is taken into account by modeling them as some arrangement of interaction sites. For reasons of efficiency, we take only points at the surface of the colloid. Furthermore, the coupling constant $`\zeta `$ is chosen rather large, which approximates stick boundary conditions. It should be noted that the dissipative nature of the coupling does not model any squeezing-out of solvent, hence the sphere is filled with the same amount of solvent regardless if it is hollow or if it were filled with particles. The solvent within the sphere follows the motion of the surrounding shell (with some minor time lag, see below), and therefore one can view the solvent inside the sphere as just belonging to the colloidal particle. Adding further particles in the sphere’s volume would have no effect except coupling the fluid within the sphere even tighter to its motion. Within our desired level of accuracy, this turned out not to be necessary. The long-time motion of the sphere is the same as that of a corresponding hard sphere, as we will show in the present paper.

We analyze basic dynamic properties of such a model colloid, where we restrict ourselves to the case of a single sphere. A detailed analysis is done for the neutral case, while some preliminary data for the case of a charged sphere are presented. We show that the (neutral) model exhibits the essential dynamical features of a spherical colloidal particle in a liquid. In our view, our method provides comparable efficiency to Ladd’s approach, while being quite straightforward to implement. Furthermore, it provides substantial flexibility with respect to the properties of the colloidal surface, namely, deformable, permeable, and non-stick surfaces can be easily simulated.

The remainder of this article is organized as follows: In Sec. 2, we describe our simulation model, while Sec. 3 contains the numerical results on translational and rotational diffusion. Finally, Sec. 4 concludes with a brief summary.

## 2 Model

Our hybrid simulation method involves two subsystems: the solvent that is modeled via LB with fluctuating stress tensor (i. e. we run a constant-temperature version of the LB method) and a Langevin MD simulation for the particles immersed in the solvent. The LB simulation is performed using the 18-velocity model , using the protocol described in . The fluid simulation consists of collision and propagation steps, the former being performed with inclusion of the momentum transfer from the solute particles (surface beads, and, for charged systems, ions).

The colloidal particle is represented by a two-dimensional tethered bead-spring network consisting of $`100`$ beads, which is wrapped around a ball of a radius $`\sigma _{cs}`$ (for notation, see below), so that the whole construction resembles a raspberry (see Fig. 1). The network connectivity is maintained via finitely extendable nonlinear elastic (FENE) springs,

$$V_{FENE}(r)=\frac{kR_0^2}{2}\mathrm{ln}\left(1\left(\frac{r}{R_0}\right)^2\right),$$
(1)
where $`k`$ is the spring constant, and $`R_0`$ the maximum bond extension. Furthermore, the beads repel each other by a modified Lennard-Jones (LJ) potential

$$V_{LJ}(r)=\{\begin{array}{cc}4ϵ_{ij}\left(\left(\frac{\sigma _{ij}}{r}\right)^{12}\left(\frac{\sigma _{ij}}{r}\right)^6+\frac{1}{4}\right)\hfill & r<2^{1/6}\sigma _{ij}\hfill \\ 0\hfill & r2^{1/6}\sigma _{ij}.\hfill \end{array}$$
(2)
An additional repulsive LJ bead is introduced at the center of the sphere in order to maintain its shape. In Eq. 2, $`i,j`$ denote either a central (“c”) or a surface (“s”) bead. The unit system is completely defined by the surface bead parameters by setting $`ϵ_{ss}`$, $`\sigma _{ss}`$, and the surface bead mass $`m_s`$ to unity. The interaction between the central bead and the surface beads is described by $`\sigma _{cs}=3`$, which is thus the sphere radius, and $`ϵ_{cs}=8`$. Furthermore, the FENE spring constant for the surface beads is $`k=300`$ and the maximum bond extension is $`R_0=1.25`$. To simulate a charged colloidal particle, we place the charge at the central bead, and add an appropriate number of counterions (LJ beads with “s” properties) outside the sphere. The electrostatic interaction is taken into account via the Coulomb potential

$$V_{el}(r)=\lambda _Bk_BT\frac{q_iq_j}{r}$$
(3)
between the various charges, where the standard Ewald summation technique  is applied. In Eq. 3, $`\lambda _B=e^2/\left(4\pi \epsilon _0\epsilon _rk_BT\right)`$ is the Bjerrum length, $`k_B`$ the Boltzmann constant, $`q_i`$ the charge of species $`i`$ in units of the elementary charge $`e`$, and $`T`$ the temperature.

The LB lattice constant is chosen as one (in our LJ unit system), and the fluid is simulated in a cubic box with periodic boundary conditions. The force between the LB fluid and the surface beads (or ions) is given by

$$\stackrel{}{F}=\zeta \left(\stackrel{}{V}\stackrel{}{u}\right)+\stackrel{}{f}.$$
(4)
Here, $`\zeta `$ is the “bare”  friction coefficient, $`\stackrel{}{V}`$ and $`\stackrel{}{u}`$ are the velocities of the bead and the fluid (at the position of the bead), respectively, while $`\stackrel{}{f}`$ is a Gaussian white noise force with zero mean, whose strength is given via the standard fluctuation-dissipation theorem  to keep the surface beads and ions at the same temperature as the solvent. The central bead is not coupled to the solvent (here $`\zeta =\stackrel{}{f}=0`$); as discussed in the Introduction, the behavior of the model would change only marginally if such a coupling were included. In our simulation we used a friction constant $`\zeta =20`$, a temperature $`k_BT=1`$, a fluid mass density $`\rho =0.85`$, and a kinematic viscosity $`\nu =3`$, resulting in a dynamic viscosity $`\eta =2.55`$. At least 20000 MD steps were performed to equilibrate the initial random bead configuration before the interaction with the LB solvent was turned on. Further details on the method can be found in .

## 3 Results

We test the simulation method against basic relations for an isolated sphere in solvent. First, we look at the center of mass’ velocity relaxation. The simplest experiment to perform is a “kick”. The sphere is placed in a LB fluid at rest (i. e. without thermal noise), and at time $`t=0`$ all particles of the sphere are assigned an identical velocity $`\stackrel{}{V}=1`$ in $`x`$ direction. Fig. 2 monitors the time behavior of the sphere’s center of mass velocity, normalized by the initial value. According to linear response theory, this relaxation function must be identical to the normalized center-of-mass velocity autocorrelation function for Brownian motion in thermal equilibrium, if the initial kick is weak enough. This is indeed satisfied, as a comparison of the two curves in Fig. 2 shows. For the experiment in thermal equilibrium, we performed 10 runs with different random number generator initializations, in order to reduce the statistical uncertainty.

It is well-known that simulations of Brownian motion in a hydrodynamic solvent are always strongly affected by finite size effects. The diffusion constant, and therefore also the relaxation function, depend on the linear system size $`L`$ due to hydrodynamic interactions with the periodic images. The diffusion constant exhibits a finite-size correction of order $`R/L`$ . Asymptotic behavior can therefore only be expected for $`R/L1`$, and this is why we performed the experiment in a rather large box of size $`L=80`$. For the same reason, the equivalence between “kick” experiment and Brownian motion will only hold if the comparison is done for the same box sizes.

One clearly sees that the initial decay of the relaxation function is characterized by two relaxation processes, one initial fast decay followed by a somewhat slower relaxation. Qualitatively, this may be explained as follows: A compact sphere of radius $`R`$ and mass $`M`$ should exhibit a velocity relaxation which, after transient ballistic motion, is initially characterized by an exponential decay $`\mathrm{exp}\left(t/\tau \right)`$, with a relaxation time $`\tau =M/\zeta _{tot}`$. Here, $`\zeta _{tot}`$ is the total friction coefficient, which we estimate via Stokes’ law for stick boundary conditions as $`\zeta _{tot}6\pi \eta R144`$. However, the effective mass is time-dependent: While initially only the mass of the beads $`M=101`$ contributes, at later times the fluid within the sphere is dragged as well, such that then the mass is roughly estimated as $`M101+4\pi \rho R^3/3214`$. This gives rise to the initial and final relaxation times $`\tau _{in}0.7`$ and $`\tau _{fin}1.5`$. However, the initial relaxation time cannot be observed, since in the extreme short-time regime ballistic effects play a role. Conversely, the decay with $`\tau _{fin}`$ is clearly visible (see Fig. 2).

After $`t1`$, the famous long-time tail  (normalized relaxation function $`V(t)/V(0)=Bt^{3/2}`$) sets in. The physical mechanism of this slow relaxation is the fact that the momentum is conserved and hence transported away diffusively from the particle . For this reason, it can only be observed in a sufficiently large box. The prefactor of the power law is known in the colloidal limit, $`B=(1/12)(Nm/\rho )(\pi \nu )^{3/2}`$, where $`N`$ is the number of beads, and $`m`$ the bead mass . It should be noted that the prefactor of the power-law decay of the unnormalized velocity autocorrelation function does not depend on the properties of the sphere at all, but only on the temperature, and the hydrodynamic properties of the solvent . The mass dependence of $`B`$ results only from the normalization, i. e. the value of $`V^2`$ according to the equipartition theorem, $`V^2=3k_BT/(Nm)`$. From this consideration it is clear that the short-time value of the sphere’s mass ($`Nm`$) enters (at $`t=0`$ the shell and the inner fluid are not yet coupled). For our simulation parameters, we find $`B=0.342`$. As Fig. 2 shows, the data exhibit the expected behavior very nicely.

Figure 3 illustrates the finite size effect in the diffusive properties by plotting the mean square displacement of the colloidal sphere as a function of time, for two different box sizes $`L=6.33R`$ and $`L=20R`$, for Brownian motion at $`k_BT=1`$. In the long-time regime, the slope is given by $`6D^S`$, where $`D^S`$ is the self-diffusion coefficient. The figure clearly shows that $`D^S`$ increases with box size. This can be explained in terms of hydrodynamic interactions with the periodic images, or, equivalently, in terms of suppression of long-wavelength hydrodynamic modes . We have therefore performed a systematic finite-size analysis, and measured $`D^S`$ for various box sizes by integrating the velocity autocorrelation function over time . In Fig. 4 we plot $`D^S/D_0`$ ($`D_0`$ denoting the asymptotic Stokes-Einstein value $`D_0=k_BT/(6\pi \eta R)=6.9\times 10^3`$) as a function of $`1/L`$; the expected linear behavior  is clearly seen.

We now look at the relaxation of the rotational motion. We performed a similar kick experiment as described above, where now an initial angular velocity $`\omega _0=1`$ was provided to the sphere. The data for the normalized decay function $`\omega (t)/\omega (0)`$ are presented in Fig. 5. The sphere dynamics shows a characteristic “raw egg” (damped) rotation pattern with a fast initial decay and a subsequent slower one. For the rotational relaxation of a hard sphere with moment of inertia $`I`$ and rotational friction coefficient $`8\pi \eta R^3`$, theory predicts a decay according to the Debye law $`\omega (t)/\omega (0)=\mathrm{exp}\left(t/\tau _r\right)`$, where the relaxation time is given by $`\tau _r=I/(8\pi \eta R^3)`$ . Similarly to the effective mass, the effective moment of inertia is expected to increase as a function of time, due to the time-delayed dragging of the fluid. Initially we expect a hollow-sphere value of roughly $`I_{in}(2/3)MR^2600`$, where we used the mass of the outer shell $`M=100`$. A more accurate value is obtained by direct summation over the contributing beads, yielding $`I_{in}=546`$. In the late-time regime the moment of inertia is expected to be the compact-sphere value $`I_{fin}(2/5)MR^2770`$, where the total mass $`M214`$ has been used. This results in relaxation times $`\tau _{r,in}=0.34`$ and $`\tau _{r,fin}=0.44`$. These values are roughly consistent with a fit to the data in the interval $`0.1<t<1.0`$, which yields a somewhat larger relaxation time $`\tau _r=0.68`$ (see Fig. 5). Given the general inaccuracy of the estimates, and the crossover of the Debye relaxation function both to short-time ballistic behavior, and long-time hydrodynamic behavior, this deviation is not too surprising.

Similar to the translational motion, the exponential decay is then followed by a power-law long-time tail. Theory predicts $`\omega (t)/\omega (0)=(\pi I/\rho )(4\pi \nu t)^{5/2}`$ . As we have seen before, the long-time tail in the translational case is governed by the short-time mass as a result of normalization. Similarly, the rotational tail must be controlled by the short-time moment of inertia $`I_{in}`$, which also governs the mean square fluctuations of the angular velocity via the equipartition theorem $`\omega ^2=3k_BT/I_{in}`$. For this reason, we can determine $`I_{in}`$ by fitting a $`t^{5/2}`$ law to the data, which is much more accurate than our rough geometrical estimate. The fit results in $`I_{in}=533`$, which is reasonably consistent. If we insert this value into the relaxation time expression, we obtain $`\tau _r=0.64`$, in quite good agreement with the data. It hence seems that fluid dragging effects are not yet very important for the initial Debye relaxation.

At longer times, one can again notice a significant finite-size effect: the curves obtained in the smaller simulation box depart from the asymptotic power-law line earlier, i. e. the long-time diffusion is hindered at the small system sizes. In fact, the power-law regime is inaccessible at $`L=10`$ while for $`L=100`$ it extends up to $`t=200`$. This value characterizes the interval after which the particle starts feeling its own periodic images. The rotational diffusion constant $`D^R`$ is given by the Green-Kubo integral

$$D^R=\frac{1}{3}_0^{\mathrm{}}𝑑t\stackrel{}{\omega }(t)\stackrel{}{\omega }(0).$$
(5)
We can evaluate this by again making use of linear response theory: The correlation function $`\stackrel{}{\omega }(t)\stackrel{}{\omega }(0)`$ is identical to the relaxation function presented in Fig. 5, multiplied with the initial value $`\omega ^2`$, which we know from the equipartition theorem (see above). Using this approach, we have calculated $`D^R`$ for different box sizes $`L`$. In an infinite fluid, $`D^R`$ has the Stokes-Einstein-Debye value $`D^R=k_BT/(8\pi \eta R^3)=0.58\times 10^3`$, towards which the results indeed converge for $`L\mathrm{}`$. Again a $`1/L`$ behavior (but weaker than for translational diffusion) is observed, as shown in Fig. 4. It should be noted that the accuracy of the data is slightly hampered by the fact that we had to use a somewhat arbitrary criterion for cutting off the integral function, which does not arrive at a constant value but rather ends up with a linear increase originating from a small nonzero constant angular velocity in the long-time limit (as a manifestation of the finite-size effect).

Finally, we have started to study the effect of charge on the self-diffusion of the colloidal particle. We performed simulations of a sphere with central charge $`Z_C=10`$ in an LB box of size $`L=40`$. Ten counterions of charge $`1`$ were also added. The Bjerrum length was set to $`2`$. Technically, the simulation ran without any problems just as well as for the neutral system. In Fig. 6 we compare the decay of the velocity autocorrelation function for the neutral and charged spheres. The difference between the two is not detectable within our error bars at this charge. We cannot say yet much about the influence of the charges on the dynamic properties; we do however expect that for strongly charged systems the Debye layer will effectively increase the hydrodynamic radius and slow down the diffusion. More work needs to be done to resolve this issue.

## 4 Summary

In summary, we introduced and tested a new model for simulating colloidal dynamics with inclusion of the hydrodynamic effects on the level of the lattice-Boltzmann equation. The suggested “raspberry”-like colloidal object exhibits the essential features of the diffusion of a spherical colloidal particle. For the first time we combined a model containing explicit charges with accurate treatment of the hydrodynamics. We expect this method to have some advantages over the previously applied scheme  which builds upon the model of Ladd  and has the problem of charge leakage. In our method: (i) the stick boundary conditions are realized not strictly, i. e. the particle surface is softly bound to the fluid; (ii) the method is easily adaptable to the multiple timestep method where the ionic action on the colloid can be averaged out somewhat within the integration step for the solvent. We should note that such a hybrid simulation requires still some computational effort — two million hybrid steps for the colloid consisting of 100 beads in a box of $`20\times 20\times 20`$ LB cells take about 12 hours on a 2 GHz Pentium 4 processor, while it is expected to run a few times faster on more sophisticated hardware. Here, one step is a MD step of the total system of beads, plus one LB update of the whole lattice, where the latter part completely dominates the CPU effort. Some improvement is possible by not updating the LB fluid at every MD step, as it was done here. For dense systems, where the MD part is no longer negligible, one should also optimize the interaction potentials. Finally, we would like to mention that the algorithm can be fully parallelized and we are working on the parallel version.

We thank Christian Holm, Jürgen Horbach and Kurt Kremer for stimulating discussions, and the latter also for critical reading of the manuscript. This work was funded by the SFB TR 6 of the Deutsche Forschungsgemeinschaft.

## 5 References


text:
# Heavy–light decay constants from the step scaling methodtalk given at Lattice 2003

## 1 Introduction

The decay constants $`f_h\mathrm{}`$ of the heavy–light pseudoscalar mesons depend upon two different energy scales. A first energy scale is the one associated with the heavy quark mass $`m_h`$, and a second one is related to the light quark mass $`m_{\mathrm{}}`$, or $`\mathrm{\Lambda }_{QCD}`$. As a consequence, numerical simulations of heavy–light decay constants require tiny lattice spacings, in order to properly describe the dynamics of the heavy quark (highly localized), and large physical volumes, to accommodate the light quark (widely spread). In practice, a direct reliable simulation of $`f_h\mathrm{}`$ should demand a lattice size of around $`L/aO(100)`$, which is hardly affordable nowadays in quenched simulations as well as in future unquenched ones. Different approaches are available on the market to face this problem . I describe a new method which has been developed by the APE group at the University of Rome “Tor Vergata” , called the step scaling method at Lattice 2002 conference .

## 2 Step scaling method

The step scaling method is a two step procedure where the heavy–light decay constants $`f_h\mathrm{}`$ are first computed on a small unphysical volume with $`L_0=0.4`$ fm, and then evolved towards larger volumes in order to remove finite size effects. The evolution is simply provided by the identity

$$f_h\mathrm{}(L_{\mathrm{}})=f_h\mathrm{}(L_0)\frac{f_h\mathrm{}(L_1)}{f_h\mathrm{}(L_0)}\frac{f_h\mathrm{}(L_2)}{f_h\mathrm{}(L_1)}\mathrm{}$$
(1)
where $`L_0<L_1<L_2\mathrm{}`$ The basic ingredient of the procedure is the computation of the step scaling function

$$\sigma (m_{\mathrm{}},m_h,L_{k1})=\frac{f_h\mathrm{}(m_{\mathrm{}},m_h,L_k)}{f_h\mathrm{}(m_{\mathrm{}},m_h,L_{k1})}|_{L_k=sL_{k1}}$$
(2)
i.e. the ratio of the decay constants at two different volumes and same quark masses. We fix $`s=2`$ as the ratio of subsequent sizes and, after two evolution steps, we end up on $`L_2=1.6`$ fm, which we take as the “infinite” volume, i.e. free of finite size effects. At fixed lattice spacing, both $`f_h\mathrm{}(L_0)`$ and the $`\sigma `$’s are simulated for a set of different quark masses, measured in a RGI quark mass scheme. Then, they are separately extrapolated to the continuum limit. Units are fixed through the $`r_0`$ scale . The main advantage of the method is that $`f_h\mathrm{}(L_0)`$ can be simulated at physical heavy quark masses without big lattice artifacts while the $`\sigma `$’s preserve in all steps a soft dependence upon the heavy quark mass, due to cancellations among volume depending terms. This makes the extrapolation to the physical values of $`m_c^{RGI}`$ and $`m_b^{RGI}`$ numerically safe. Moreover, such values can be obtained by a self–consistent application of the step scaling method to the HL–meson masses .

## 3 Computational framework

The calculation is set up on the Schrödinger Functional scheme, with topological parameters

$$T=2L,C=C^{}=0,\theta =0$$
(3)
$`C`$ and $`C^{}`$ being the gauge fields on time boundaries, and $`\theta `$ being a phase which affects the periodicity of the fermion spatial boundary conditions. We simulate the basic SF correlation functions, i.e. the axial, pseudoscalar and boundary normalization bilinears

$$F_A^I(x_0)=\frac{a^6}{2}\underset{\mathrm{𝐲𝐳}}{}\overline{\zeta }_j(𝐲)\gamma _5\zeta _i(𝐳)A_0^I(x)$$
(4)
$$F_P(x_0)=\frac{a^6}{2}\underset{\mathrm{𝐲𝐳}}{}\overline{\zeta }_j(𝐲)\gamma _5\zeta _i(𝐳)P(x)$$
(5)
$$F_1=\frac{a^{12}}{3L^6}\underset{\mathrm{𝐲𝐳𝐲}^{}𝐳^{}}{}\overline{\zeta }_j(𝐲)\gamma _5\zeta _i(𝐳)\overline{\zeta }_i^{}(𝐲^{})\gamma _5\zeta _j^{}(𝐳^{})$$
(6)
In terms of these, a finite size definition of the decay constants is given by taking the axial correlation function in the middle of the lattice

$$f_h\mathrm{}=\frac{2}{\sqrt{L^3M_A(T/2)}}Z_A(1+b_Aam)\frac{F_A^I(T/2)}{\sqrt{F_1}}$$
(7)
and the normalization is provided in terms of the effective meson mass, still measured in a finite size fashion

$$aM_A(T/2)=\frac{1}{2}\mathrm{ln}\left[\frac{F_A^I(T/2a)}{F_A^I(T/2+a)}\right]$$
(8)
Quark flavors are defined according to their RGI masses (scale/finite–size independent), defined as follows. We start from the bare quark masses

$$am_i=\frac{1}{2}\left[\frac{1}{\kappa _i}\frac{1}{\kappa _c}\right]$$
(9)
$$m_{ij}^{WI}=\frac{(1/2)(_0+_0^{})F_A+ac_A_0^{}_0F_P}{2F_P}|_{T/2}$$
(10)
and connect them to the RGI quark masses through proper renormalisation factors known in literature

$$m_{ij}^{RGI}=Z_M\left[1+(b_Ab_P)\frac{am_i+am_j}{2}\right]m_{ij}^{WI}$$
(11)
$$\widehat{m}_i^{RGI}=Z_MZ[1+b_mam_i]m_i,Z=\frac{Z_mZ_P}{Z_A}$$
(12)
In this way, we get a plethora of definitions, all with the same continuum limit and differing in $`O(a^2)`$ terms:

* $`m_i^{RGI}=m_{ii}^{RGI}`$
* $`m_{i_{(j)}}^{RGI}=2m_{ij}^{RGI}m_{jj}^{RGI},ji`$
* $`\widehat{m}_i^{RGI}`$

After choosing a target pair $`(m_i^{RGI},m_j^{RGI})`$, we tune the hopping parameters $`(\kappa _i,\kappa _j)`$ such that the values of the masses coming from the various definitions equate the chosen pair. Of course, each choice of $`(\kappa _i,\kappa _j)`$ produce a different value of $`f_h\mathrm{}`$, all the values differing in $`O(a^2)`$ terms, as it can be seen in Fig.1.

## 4 Results at $`L_0=0.4`$ fm

Simulations of the decay constants on the smallest volume $`L_0=0.4`$ fm have been performed at five different lattice spacings, using the geometries $`24\times 12^3`$, $`32\times 16^3`$, $`40\times 20^3`$, $`48\times 24^3`$ and $`64\times 32^3`$. For each lattice spacing a set of eight RGI quark masses have been simulated, two around $`m_b^{RGI}`$, two around $`m_c^{RGI}`$, one additional at 4.0 GeV, and three light quark masses at 140, 100 and 60 MeV. We fix the light quark mass at the physical value of $`m_s^{RGI}`$ by interpolation and we get

$$f_{B_s}(L_0)=475(2)\mathrm{MeV}$$
(13)
$$f_{D_s}(L_0)=644(3)\mathrm{MeV}$$
(14)
The errors quoted at this stage are statistical only, and they have been evaluated by a jackknife procedure. The continuum extrapolation has been obtained by performing a combined linear fit in $`(a/r_0)^2`$ over the last three sets of points of Fig.1. We estimate a systematic error of 1%, related to this extrapolation, by the discrepancy observed after restricting the fit to the two sets of points nearest to the continuum. The results obtained are clearly unphysical, as they have to be corrected through the evolution factors given by the $`\sigma `$’s.

## 5 Evolution steps

The first evolution step connects $`L_0`$ to $`L_1=2L_0=0.8`$ fm. Three different lattice spacings have been simulated in order to perform the continuum extrapolation, corresponding to the geometries with $`L_0/a=8,12,16`$ and consequently $`L_1/a=16,24,32`$. The values of the simulated quark masses have been halved in order to leave the lattice artifacts unchanged. In Fig. 2 the continuum step scaling function $`\sigma (L_0)`$ is plotted vs. the inverse of the heavy quark mass with the light quark mass fixed at $`m_s^{RGI}`$.

At this stage we get

$$\sigma _{B_s}(L_0)=0.417(3)$$
(15)
$$\sigma _{D_s}(L_0)=0.414(3)$$
(16)
Analogously, the second evolution step connects $`L_1`$ to $`L_2=2L_1=1.6`$ fm. Also here, three lattice spacings have been simulated with the same lattice geometries as before. Again the values of the simulated quark masses have been halved. Fig. 3 shows the continuum step scaling function $`\sigma (L_1)`$ vs. the inverse of the heavy quark mass at $`m_l^{RGI}=m_s^{RGI}`$. We find

$$\sigma _{B_s}(L_1)=0.97(3)$$
(17)
$$\sigma _{D_s}(L_1)=0.90(2)$$
(18)
## 6 Results

Combining the results of the small volume $`L_0`$ with the evolution factors, we obtain

$$f_{B_s}=192(6)(4)\mathrm{MeV}$$
(19)
$$f_{D_s}=240(5)(5)\mathrm{MeV}$$
(20)
The first error is statistical while the second one is systematic and due to the uncertainties on the continuum extrapolation, on the scale and on the renormalisation factors.


text:
# Coulomb Correlations and Orbital Polarization in the Metal Insulator Transition of VO2

## I Introduction

The metal insulator transition in VO2 has been intensively studied for a long time. At 340 K the resistivity changes by several orders of magnitude.morin  The high-temperature metallic phase has a rutile structure, while the low-temperature insulating phase is monoclinic ($`M_1`$), with a zigzag-type pairing of V atoms along the $`c`$ axis. Both phases are non-magnetic. Although this transition is widely regarded as a Mott-Hubbard transition, pouget ; zilbersztejn ; mott ; rice  the role of the Peierls distortion of the crystal structure in the insulating phase has been the topic of intense debate.goodenough ; caruthers ; wentzcovitch ; eyert ; huang ; korotin ; continenza ; laad  The discovery of a second insulating phase ($`M_2`$), pouget  in which only half of the V atoms form pairs while the other evenly spaced chains behave as magnetic insulators, suggested that both low temperature phases should be regarded as Mott-Hubbard and not as Peierls band insulators.zilbersztejn ; rice  Although the role of the Coulomb interaction in the metal insulator transition of VO2 has been studied previously, huang ; korotin ; continenza ; laad  a consistent description of the rutile and monoclinic phases is not yet available.

The aim of the present work is to elucidate the interplay of Coulomb correlations and orbital polarization in the quasi-particle spectra of VO2. Analyzing recent photoemission dataokazaki ; tjeng  we demonstrate that the metallic and insulating phases show evidence of strong local interactions which manifest themselves in distinctly different ways because of the different degree of orbital polarization in the rutile and monoclinic structures. The key aspect of the metallic phase is the small orbital polarization among V $`3d`$ $`t_{2g}`$ states, implying weak static but strong dynamical correlations. Thus, the spectra reveal band narrowing and reduced weight of the coherent peak near the Fermi level, and an incoherent satellite feature associated with the lower Hubbard band. The monoclinic phase, in contrast, is characterized by a pronounced $`t_{2g}`$ orbital polarization induced by the V–V Peierls distortion. The local Coulomb interaction therefore leads to strong static correlations and to the opening of an excitation gap. The influence of non-diagonal coupling among $`t_{2g}`$ orbitals on the size of the gap is investigated within the static limit and is found to be small.

In the following sections we discuss results obtained within several theoretical approaches, such as the dynamical mean field theory (DMFT)DMFT , the local density approximation plus Hubbard $`U`$ (LDA+U)LDA+U  and $`GW`$GW  methods. We argue that none is presently capable of explaining all of the observed phenomena in a consistent manner. Instead, we focus on the merits and failures of these approaches with the aim of highlighting the different roles of Coulomb correlations and orbital polarization in the metallic and insulating phases of VO2.

## II Results and Discussion

Fig. 1 shows photoemission spectra for VO2 films (150 Å thick) grown on TiO2(001) at binding energies corresponding to the V $`3d`$ bands (21.2 eV photon energy).okazaki  Emission from O $`2p`$ states extends from $`2`$ eV to about $`8`$ eV. As a result of compressive stress from the substrate muraoka  the transition temperature is lowered from 340 K to 290 K. The high temperature spectrum shows a Fermi edge characteristic of metallic behavior which is absent in the low temperature spectrum. While this trend is in agreement with previous photoemission measurements,PES ; shin  the data in Fig. 1 reveal two spectral features in the metallic phase: a peak close to $`E_F`$ and a second one near $`1`$ eV. In the insulating phase, the peak near $`E_F`$ disappears and the feature near $`1`$ eV becomes more intense. Recent VO2 photoemission spectra taken at 520 eV photon energy tjeng  agree with the ones shown in Fig. 1, except for a considerably greater relative weight of the coherent peak near $`E_F`$ in the metallic phase. These spectral changes are consistent with the observation that satellite peaks in low photon energy spectra tend to be more pronounced as a result of a surface induced enhancement of Coulomb correlations. maiti ; sekiyama ; liebsch03

Before analyzing the photoemission data we discuss the single particle properties of VO2 obtained using density functional theory. We have carried out full potential linearized augmented plane wave (LAPW) calculations for the rutile and monoclinic structures using the experimental lattice parameters and treating exchange correlation within the generalized gradient approximation (GGA).perdew  Due to the octahedral crystal field, the states near $`E_F`$ have V $`3d`$ $`t_{2g}`$ character. They are separated by a small gap from the empty $`e_g`$ states, and from the O $`2p`$ states by a gap of about 1.0 eV. The occupancy of the $`t_{2g}`$ manifold is $`3d^1`$. Our results qualitatively confirm previous LDA calculations. wentzcovitch ; eyert ; huang ; korotin  Fig. 2(a) shows a comparison of the V total $`t_{2g}`$ density of states for the rutile and monoclinic phases of VO2. Although there are differences in detail, the overall width of these distributions and the shape of the occupied region are similar for both structures. Evidently, the GGA/LDA does not predict the insulating nature of the monoclinic phase. Moreover, on the basis of these densities one would not expect correlations to play very different roles in the two phases. However, if we plot the subband contributions to the $`t_{2g}`$ density, the two structures are very different, as shown in Figs. 2(b) and (c). Whereas in the rutile phase the $`t_{2g}`$ bands have similar occupation numbers, in the monoclinic phase the $`d_{x^2y^2}`$ band is significantly more occupied than the $`d_{xz,yz}`$ bands. (We adopt the local coordinate system of Ref. goodenough, , i.e., $`x`$ denotes the $`c`$ axis, while $`y`$ and $`z`$ point along the diagonals of the $`a,b`$ plane.) The origin of this orbital polarization is the V–V dimerization along the $`c`$ axis. The $`d_{x^2y^2}`$ band splits into bonding and anti-bonding components, while the $`d_{xz,yz}`$ bands are pushed upwards due to shortening of V-O distances. Investigation of the energy bands shows that, in agreement with earlier work,wentzcovitch ; eyert  the top of the bonding $`d_{x^2y^2}`$ bands is separated by a slight negative gap from the bottom of the $`d_{xz,yz}`$ bands. In the following subsection we show that different degree of orbital polarization in the rutile and monoclinic phases has a pronounced effect on the quasi-particle spectra of VO2.

### II.1 Metallic Rutile Phase

Let us first discuss the rutile phase. Comparing the $`t_{2g}`$ density of states with the photoemission spectra it is plausible to associate the feature close to $`E_F`$ with emission from metallic V $`3d`$ states. The peak near $`1`$ eV, however, lies in the gap between V $`3d`$ and O $`2p`$ states and cannot be understood within the single particle picture. To describe the spectra in the rutile phase it is clearly necessary to account for dynamical Coulomb correlations. For the evaluation of the quasi-particle distributions we use the Dynamical Mean Field Theory combined with the multiband Quantum Monte Carlo (QMC) method. DMFT  Since hybridization among $`t_{2g}`$ states is weak the local self-energy is taken as diagonal in orbital space. The $`t_{2g}`$ density of states components shown in Fig. 2(b) then serve as input quantities accounting for the single particle properties of the rutile structure. The local Coulomb interaction defining the quantum impurity problem is characterized by intra- and inter-orbital matrix elements $`U`$ and $`U^{}=U2J`$, where $`J`$ is the Hund’s rule exchange integral. According to constrained LDA calculations, $`U4.2`$ eV and $`J0.8`$ eV.korotin ; UJ  The calculations are performed at $`T500`$ K with up to $`10^6`$ sweeps. The quasi-particle distributions are obtained via maximum entropy reconstruction. jarrell

Fig. 3(a) shows that, in contrast to the single particle density of states, the calculated $`t_{2g}`$ quasi-particle spectra for the rutile phase of VO2 exhibit two spectral features, a coherent peak near $`E_F`$ and a lower Hubbard band near $`1`$ eV, in agreement with experiment. The peak near $`E_F`$ accounts for the band narrowing and lifetime broadening of the metallic states whereas the Hubbard band is associated with incoherent emission. Since the $`t_{2g}`$ subbands have comparable single particle distributions their quasi-particle spectra reveal similar correlation features. Moreover, because of the weak orbital polarization in the rutile structure static correlations are negligible. Thus, in the metallic phase the spectral weight transfer between coherent and incoherent contributions to the spectrum is primarily the result of dynamical correlations.

We point out that, in view of the approximate nature of the model underlying the DMFT, quantitative agreement with photoemission data cannot be expected. On the theoretical side, the consideration of purely on-site Coulomb interactions and the neglect of the momentum variation of the self-energy permit only a qualitative analysis of the spectra. In addition, there exists some uncertainty regarding the precise values of the Coulomb and exchange energies. Finally, the DMFT results depend on the temperature used in the QMC calculation. The comparison with results obtained for slightly different values of $`U`$, $`J`$ and $`T`$, however, gives us confidence that in the metallic phase the transfer of spectral weight from the coherent peak near $`E_F`$ to the satellite region near $`1`$ eV is qualitatively reliable and consistent with analogous dynamical correlation effects in other $`3d^1`$ transition metal oxides, such as SrVO3. sekiyama ; liebsch03 ; pavarini  As we discuss below, the local DMFT treatment predicts the monoclinic phase to be also metallic. The metal insulator transition in VO2 is therefore not achieved simply by lowering the temperature in the rutile phase. The lattice transformation from the rutile to monoclinic structure must be taken into account. Therefore, the DMFT results shown in Fig. 3(a) for the rutile structure at $`T=500`$ K can be regarded as representative of correlation induced behavior in the metallic domain.

On the experimental side, as pointed out above, it is important to recall that photoemission data taken at low photon energies represent a superposition of bulk and surface contributions. Since correlation effects are observed to be more enhanced at surfaces, maiti ; sekiyama ; liebsch03  the relative intensity of the satellite peak near $`1`$ eV in the 21.2 eV spectra shown in Fig. 1 is considerably more pronounced than at high photon energies tjeng  at which primarily bulk-like valence states are detected.

Dynamical effects in the metallic phase may also be evaluated within the $`GW`$ approach GW  which treats long range Coulomb interactions in the random phase approximation (RPA) and which has proven rather useful to describe excitation spectra of weakly correlated systems. GWreview  Because of the neglect of multiple electron-electron and hole-hole scattering processes, this scheme fails when local Coulomb interactions are important.gunnarsson  Presumably, therefore, for VO2 the GW method does not reproduce the lower Hubbard band. The satellite in the rutile phase is, of course, also beyond the scope of the static LDA+U approach.comment  The qualitative agreement between the measured high-temperature spectra and the theoretical results shown in Fig. 2(a) suggests that the DMFT captures the key spectral weight rearrangement induced by dynamical correlations.

### II.2 Insulating Monoclinic Phase

Turning now to the monoclinic phase we first calculated the $`t_{2g}`$ quasi-particle spectra within the DMFT. Because of the orbital polarization caused by the V–V Peierls distortion, correlation effects in the $`d_{x^2y^2}`$ band are now stronger so that the intensity of the lower Hubbard band is enhanced compared to the rutile phase (not shown). While this trend agrees with the photoemission data, the coherent peak near $`E_F`$ persists. Thus, the single-site DMFT based on a diagonal self-energy does not reproduce the insulating nature of the monoclinic phase.heldVO2  This failure is in striking contrast to the results obtained within the LDA+U and $`GW`$ methods.huang ; continenza

Let us discuss first the LDA+U approach. In this scheme different orbital occupations give rise to orbital dependent potential terms which shift different $`t_{2g}`$ bands in opposite directions. Thus, orbital polarization leads to strong static correlation effects. Using the same Coulomb parameters as for the rutile phase we find that the local density of states calculated via the LAPW/LDA+U exhibits a gap of about 0.7 eV (see Fig. 3(b)), in agreement with optical data shin  and previous LDA+U results. huang  The bonding $`d_{x^2y^2}`$ bands are now completely filled and the $`d_{xz,yz}`$ and anti-bonding $`d_{x^2y^2}`$ bands are shifted upwards. Thus, the LDA orbital polarization in the monoclinic phase is further increased as a result of static correlations. In turn, this reduction of orbital degeneracy enhances the trend towards insulating behavior.koch

We emphasize that the LDA+U does not provide a complete description of correlation effects in the insulating phase. Essentially, it amounts to a fully self-consistent, non-local treatment of static screening within the dimerized structure. Genuine dynamical effects, such as the spectral weight transfer from the coherent to the incoherent peak as observed in the metallic phase are ignored. Neverteless, since static correlations appear to be the origin of the excitation gap in VO2, it is instructive to inquire which features of the LDA+U contribute to the opening of the gap. Of particular interest is the role of the non-diagonal occupation matrix $`n_{\alpha \beta }`$ which is the key input quantity in the orbital dependent perturbation potential used in the LDA+U. Recent work by Pavarini et al. pavarini  on $`3d^1`$ perovskite materials exhibiting non-diagonal $`t_{2g}`$ orbital coupling caused by octahedral distortions suggests that this mechanism tends to suppress orbital fluctuations and to enhance insulating behavior. Fig. 3(c) compares the VO2 LDA+U density of states with results of an approximate LDA+U treatment in which at each iteration only the diagonal elements $`n_{\alpha \alpha }`$ are retained. The size of the gap is seen to be only slightly reduced. In the case of VO2, therefore, non-diagonal coupling among $`t_{2g}`$ orbitals is evidently not the main reason for the existence of the gap.

To understand the gap formation obtained within the LDA+U, it is useful to formally express the $`t_{2g}`$ self-energy matrix as $`\mathrm{\Sigma }(\omega ,k)=\mathrm{\Sigma }^{\mathrm{HF}}(k)+\mathrm{\Delta }\mathrm{\Sigma }(\omega ,k)`$, where $`\mathrm{\Sigma }^{\mathrm{HF}}(k)=H^{\mathrm{LDA}+\mathrm{U}}(k)H^{\mathrm{LDA}}(k)`$ is real and accounts for spectral changes associated with the LDA+U.HF-limit  $`\mathrm{\Delta }\mathrm{\Sigma }(\omega ,k)`$ is complex and describes purely dynamical effects. The important point is that the LDA+U includes the full momentum variation of $`\mathrm{\Sigma }^{\mathrm{HF}}(k)`$ within the true lattice geometry. Thus, in a site representation, static screening processes generate finite inter-site elements $`\mathrm{\Sigma }_{ij}^{\mathrm{HF}}`$ (each element is a matrix in orbital space) even if the bare LDA+U perturbation potential is site-diagonal. Thus, $`\mathrm{\Sigma }_{ij}^{\mathrm{HF}}(k)`$ is more accurate than the static limit of the single-site DMFT which neglects the $`k`$ dependence and assumes the impurity environment to be isotropic, i.e., $`\mathrm{\Sigma }_{ij}^{\mathrm{DMFT}}(\omega )\delta _{ij}`$. The results shown in Fig. 3(b,c) suggest that the proper evaluation of $`\mathrm{\Sigma }^{\mathrm{HF}}(k)`$ within the Brillouin Zone of the dimerized structure is the crucial ingredient to an adequate description of the insulating behavior in VO2. The LDA+U amounts to a self-consistent treatment of $`\mathrm{\Sigma }^{\mathrm{HF}}(k)`$ since the solution of the Schrödinger equation imposes no restrictions on how the wave functions adjust to the LDA+U potential.

In order to go beyond static screening and include dynamical correlations in the monoclinic phase a cluster extension of the DMFT is most likely required. Such an extension is beyond the scope of the present work. A cluster DMFT would include the crucial inter-site elements $`\mathrm{\Sigma }_{ij}^{\mathrm{DMFT}}(\omega )`$ which arise naturally in a cluster representation of the lattice and of the impurity Green’s functions $`G(\omega )`$ and $`G_0(\omega )`$, even for a purely on-site Coulomb interaction. Preliminary results for VO2 within a cluster DMFTsascha  show that dynamical screening processes beyond the static correlations included in the LDA+U cause a broadening of the LDA+U density distribution and a shift of the main spectral peak towards the Hubbard bands. Possibly, a multi-site extension of the DMFT might also identify the true origin of the metal insulator transition in VO2, i.e., whether it is primarily caused by the lattice reconstruction or by Coulomb correlations, or whether these mechanisms mutually enhance each other. All one can say at present is that, given the orbital polarization induced by the lattice transition, Coulomb correlations have a pronounced effect on the quasi-particle spectra.

As noted above, the $`GW`$ approach applied to the monoclinic phase of VO2 also yields an excitation gap of the correct magnitude. continenza  These calculations utilize a model self-energy gygi  consisting of an approximate short-range contribution given by the local exchange correlation potential, and a correction due to incomplete screening of the Coulomb interaction. Essentially, in this simplified $`GW`$ scheme the self-energy correction consists of a “scissor” operator and additional, non-rigid shifts of energy eigenvalues.gygi  Presumably, the reason why this model self-energy yields a gap is that $`\mathrm{\Sigma }(q,\omega )`$ is non-local and non-diagonal in site space, i.e., it includes the important static correlations in the dimerized structure in a similar fashion as the LDA+U. The approximate nature of the $`GW`$ method, in particular, the neglect of electron-electron and hole-hole ladder type interactions, affects mainly the remaining dynamical corrections caused by the strong local Coulomb energy. An adequate treatment of these corrections would require going beyond the RPA and should lead to a more realistic description of the position and width of the Hubbard bands.

An interesting additional feature observed by Okazaki et al. okazaki  is the temperature dependence of the photoemission spectra below the metal insulator transition. Essentially, towards lower $`T`$ the excitation gap becomes more clearly defined and the main peak near $`1`$ eV gets slightly sharper. These changes, however, are very small on the scale of the main discrepancy still existing between the LDA+U or $`GW`$ results and the experimental data. At present it is not clear whether cluster DMFT calculations in the monoclinic phase as a function of temperature will be able to explain the observed trend or whether an explicit treatment of electron-phonon coupling is required.

## III Summary

The metal insulator transition in VO2 appears to be remarkably complex and its origins are not yet fully understood. In the present work we focussed on the important role of two aspects, local Coulomb correlations and orbital polarization, in the low and high temperature photoemission spectra of VO2. Whereas the metallic phase exhibits weak static and strong dynamical correlations, the monoclinic phase is dominated by static Coulomb correlations. Accordingly, the rutile spectra reveal a double-peak structure where the feature close to $`E_F`$ is identified with metallic V $`3d`$ states and the peak near $`1`$ eV with the lower Hubbard band. Since the $`t_{2g}`$ states in the metallic phase are roughly equally occupied orbital polarization is negligible. The fundamental difference in the monoclinic phase is the large orbital polarization induced by the symmetry breaking due to V–V dimerization. The preferential occupation of the $`d_{x^2y^2}`$ bonding states implies strong static correlations which are the main origin of the excitation gap. It would be of great interest to study the additional dynamical correlation effects in this phase within an extension of the single-site DMFT approach.

A. L. likes to thank K. Okazaki, A. Fujimori, and L.H. Tjeng for sharing their photoemission data prior to publication. He also thanks F. Aryasetiawan, A. Bringer, K. Held, G. Kotliar, A. I. Lichtenstein, A. Poteryaev, and D. Vollhardt for fruitful discussions.

email: a.liebsch@fz-juelich.de; ishida@chs.nihon-u.ac.jp; g.bihlmayer@fz-juelich.de