text:
# Modeling of the Lyman 𝛾 satellites in FUSE spectra of DA white dwarfsThis paper is dedicated in memory of J. L. Greenstein, discoverer of the quasi-molecular lines in white dwarfs, who passed away on October 21, 2002.

## 1 Introduction

Quasi-molecular satellites are absorption features due to transitions that take place during close collisions of a radiating hydrogen atom with a perturbing atom or proton. These spectral features, which are present in the red wings of the Lyman series lines, provide an important source of opacity in the atmospheres of hydrogen-rich white dwarfs.

The IUE spectra of the DA white dwarf 40 Eri B were the first ones which revealed a strong, unexpected absorption feature near 1400 Å (Greenstein 1980). Absorption features near 1600 Å were thereafter detected in spectra of cooler DA white dwarfs (Holm et al. 1985). These features remained unexplained until Koester et al. (1985) and Nelan & Wegner (1985) simultaneously identified them as quasi-molecular satellites of the Lyman $`\alpha `$ line. The H2 and H$`{}_{}{}^{+}{}_{2}{}^{}`$ satellites at 1600 Å and 1400 Å were thereafter observed in other DA white dwarfs, as well as in laboratory plasmas (Kielkopf & Allard 1995). The H$`{}_{}{}^{+}{}_{2}{}^{}`$ satellite absorption features at 1058 Å and 1076 Å were first identified in the spectrum of the DA white dwarf Wolf 1346, as observed with HUT (Koester et al. 1996). ORFEUS (Koester et al. 1998) and FUSE (Wolff et al. 2001; Hébrard et al. 2002b) observations allow these Lyman $`\beta `$ satellites to be observed in other targets.

New calculations for the absorption profiles of Lyman $`\gamma `$ line of atomic hydrogen perturbed by protons have been used to predict synthetic spectra of hot hydrogen-rich white dwarfs and allow the identification of a strong feature near 995 Å. This feature was present in the HUT spectrum of Wolf$ 1346$ (Koester et al. 1996) but it is only in Koester et al. (1998) that a similar feature present in ORFEUS spectra of WD 1031$``$114 and WD 0644$`+`$375 was suspected to be a new satellite in the wing of Lyman $`\gamma `$. It was also distinctly detected in FUSE spectra of CD$ -38^∘10980$ (Wolff et al. 2001) and of Sirius B (Holberg 2002).

These last two recent observations provided motivation to investigate the theory of the Lyman $`\gamma `$ line profile using accurate theoretical molecular potentials, instead of using the simple Stark broadening approximation. Our new calculations of Lyman $`\gamma `$ satellites are summarized in Sect. 2. We present in Sect. 3 a comparison between synthetic spectra including Lyman $`\gamma `$ satellites and the FUSE spectra of four DA white dwarfs. The results are discussed in Sect. 4.

## 2 Theoretical line profiles

The theoretical approach is based on the theory of pressure broadening due to Baranger (1958), developed in an “adiabatic representation” that does not exclude degeneracy of atomic levels. A detailed description of our semi-classical approach as applied to the shape of the Lyman lines has been given by Allard et al. (1994) and Allard et al. (1999).

The Lyman $`\gamma `$ profile and its satellites are calculated for physical conditions encountered in the atmospheres of white dwarfs. Because the density of the hydrogen atoms is low (1015 to 1017 cm-3), we computed the line profile by using the low density approximation as described by Allard et al. (1994). This approximation uses the expansion of the autocorrelation function in powers of density. We also used accurate theoretical potentials describing the binary interactions of one hydrogen atom with a ionized hydrogen atom. The H$`{}_{}{}^{+}{}_{2}{}^{}`$ molecular potentials are available for the H ($`n3`$) states (Madsen & Peek 1971). We computed all those related to $`n=4`$ using a code kindly provided by J. Peek.

Our theoretical approach allows the variation of the radiative dipole moment during a collision. But through lack of appropriate data for Lyman $`\gamma `$, we assumed that the dipole moment for the transitions is constant. Such approximation has been applied to the calculations of the quasi-molecular satellites observed in the IUE and HST spectra of white dwarfs (see, e.g., Koester & Allard 1993; Koester et al. 1994; Bergeron et al. 1995). In the last section we discuss the possible effects of a variable dipole moment on the line profile calculations.

An H$`{}_{}{}^{+}{}_{2}{}^{}`$ correlation diagram was constructed for Lyman $`\gamma `$ (Allard et al. 2003), to which 14 transitions contribute. Two of these transitions, $`1s\sigma _g6h\pi _u`$ and $`2p\sigma _u7i\sigma _g`$, yield satellites in the red portion of Lyman $`\gamma `$ at 992 Å and 996 Å. Fig. 1 shows the profiles of both individual satellites and their effect on the total profile of Lyman $`\gamma `$. We obtain a blend of the line satellites with a maximum near 992 Å. The profile gives a shape similar to the 995 Å feature observed in white dwarf spectra.

## 3 Astronomical applications to synthetic spectra for DA white dwarfs

Atmosphere models have been calculated using the programs TLUSTY (Hubeny 1988; Hubeny & Lanz 1992, 1995). We considered LTE model atmospheres with pure hydrogen composition that explicitly include the Lyman $`\alpha `$ and Lyman $`\beta `$ quasi-molecular opacities. The resulting spectra were computed by using the spectral synthesis code SYNSPEC that incorporates the quasi-molecular satellites of Lyman $`\alpha `$, Lyman $`\beta `$, and Lyman $`\gamma `$. Models predict that the Lyman $`\gamma `$ satellites are visible roughly for effective temperatures between $`\mathrm{15\hspace{0.17em}000}`$ K and $`\mathrm{30\hspace{0.17em}000}`$ K. We verified that non-LTE effects are unimportant in the range of effective temperatures considered in this study.

We used observations of four DA white dwarfs performed with FUSE (Moos et al. 2000). The spectra were obtained between 2000 and 2002 and were retrieved from the FUSE public database (see Table 1). The one-dimensional spectra were extracted from the two-dimensional detector images and calibrated using the CalFUSE pipeline. When possible, we did not use exposures with strong airglow emission. We also excluded regions of the spectra near the edges of the detector segments, which provide poor quality data. The eight FUSE detector segments of the different exposures were co-added together and projected on a $`0.16`$ Å-pixel base, so pixels about 25 times larger than the original FUSE detector pixels. This degradation of the FUSE spectral resolution (typically $`\lambda /\mathrm{\Delta }\lambda 15000`$ for this kind of target with the slit LWRS or MDRS; see Hébrard et al. 2002a or Wood et al. 2002) has no effect on the shapes of the large stellar features which we study and allows the signal-to-noise ratio to be increased.

In Fig. 2 we compare the FUSE spectra with the models that we obtained for these targets. The effective temperatures are within the intermediate temperature range for the visibility of both the Lyman $`\beta `$ and Lyman $`\gamma `$ satellites. The agreement indicates, by the coincidence of position and approximate shape, that the large feature near 995 Å is indeed a satellite of Lyman $`\gamma `$ due to the blend of two close features at 992 Å and 996 Å.

The first target plotted in Fig. 2 is Wolf$`\mathrm{\hspace{0.17em}1346}`$. The comparison for the FUSE spectrum shows that the far UV is well fitted with a synthetic spectrum computed with our new profile calculations for $`\mathrm{log}g=8.0`$ and $`T_{\mathrm{eff}}=\mathrm{20\hspace{0.17em}000}`$ K. These atmospheric parameters are close to those determined from optical spectra (7.9; 20 000) and adopted by Koester et al. (1996) to fit the HUT spectrum of that white dwarf. The Lyman $`\gamma `$ satellite, already observed in the HUT spectrum, is well-reproduced by our synthetic spectrum.

The next two objects, EG$ 102$ and BPM$ 6502$, are slightly hotter white dwarfs, with smaller gravities. We used the effective temperatures and gravities determined from optical spectra (Finley et al. 1997; Bragaglia et al. 1995). The predicted spectra for EG$`\mathrm{\hspace{0.17em}102}`$, with $`T_{\mathrm{eff}}=\mathrm{20\hspace{0.17em}435}`$ K and $`\mathrm{log}g=7.87`$, and for BPM$`\mathrm{\hspace{0.17em}6502}`$, with $`T_{\mathrm{eff}}=\mathrm{21\hspace{0.17em}380}`$ K and $`\mathrm{log}g=7.86`$, fit quite well the FUSE spectra. In the case of BPM$`\mathrm{\hspace{0.17em}6502}`$, all the FUSE exposures present strong airglow emissions, so it was not possible to remove them. These two stars are not normal DA white dwarfs: EG$`\mathrm{\hspace{0.17em}102}`$ presents metallic lines (Holberg et al. 1997) and BPM$`\mathrm{\hspace{0.17em}6502}`$ is in a post common-envelope binary system (Kawka et al. 2000, 2002). These particularities seem however to have no effects on the quality of the fits.

The hottest object we present is CD$`38^{}10980`$. The FUSE observation that we used is the same as the one presented by Wolff et al. (2001), but our data reduction was performed with a more recent version of CalFUSE. Our best fit was obtained using $`T_{\mathrm{eff}}=\mathrm{23\hspace{0.17em}000}`$ K and $`\mathrm{log}g=7.9`$. The values of the atmospheric parameters are slightly smaller than the ones determined from optical spectra (Vauclair et al. 1997) to get a better agreement with the Lyman $`\gamma `$ satellite.

CD$`38^{}10980`$ is near the upper limit of the expected visibility range of the Lyman $`\beta `$ satellites. Nevertheless, Lyman $`\beta `$ and Lyman $`\gamma `$ quasi-molecular satellites have been observed in the FUSE spectrum of the ultra-massive white dwarf PG$`\mathrm{\hspace{0.17em}1658}+441`$ ($`T_{\mathrm{eff}}=\mathrm{30\hspace{0.17em}500}`$K, $`\mathrm{log}g=9.36`$) as reported by Dupuis et al. (2001). We speculate that the presence of quasi-molecular satellites in this star is related to its high surface gravity and we will investigate the question in more details in Dupuis et al. (2003). It is interesting to note that Sirius B ($`T_{\mathrm{eff}}=\mathrm{24\hspace{0.17em}790}`$K, $`\mathrm{log}g=8.57`$; Holberg et al. 1998), which also has a relatively high surface gravity, shows the quasi-molecular satellites (Holberg 2002).

Finally, the four spectra show a kink in the profile around 1037 Å that is well reproduced by our models (Fig. 2). This feature is probably the first identification of the H$`{}_{}{}^{+}{}_{2}{}^{}`$ satellite of Lyman $`\beta `$, as predicted in Fig. 7 of Allard et al. (1998a.)

## 4 Conclusion

The good agreement between the FUSE spectra and our calculations allowed the line satellites near 995 Å to be identified, although there are still some systematic discrepancies between the predicted and observed profiles. The shape of the profile in the region of the satellites is sensitive to the relative strength of these two main features. It is therefore important to get an accurate quantitative determination of the satellite amplitudes. Accurate theoretical molecular potentials have to be used to describe the interaction between the radiator and the perturber as we did in the work described here. Another important factor, not yet included, is the variation of the dipole moment during the collision. Allard et al. (1998a, 1998b, 1999) have shown that the strengths of line satellites are dependent on values of the electric-dipole moments at the internuclear separation responsible for the satellites. Large enhancements in the amplitudes of the satellites may occur whenever the dipole moment increases through the region of internuclear distance where the satellites are formed. The theoretical shape in this region of the profile is then dependent on the dipole moment of the transition. Some of the differences between the observed and the synthetic spectra may be due to our constant dipole moment approximation. That will be improved in future calculations.

As a summary we present in Fig. 3 the sum of the profiles of Lyman $`\alpha `$, Lyman $`\beta `$, and Lyman $`\gamma `$ perturbed by collisions with neutral hydrogen and protons for two different densities of neutral hydrogen. We can note that beside the quasi-molecular lines already detected, there is another predicted H2 satellite at 1150 Å (Allard et al. 2000), in the red wing of Lyman $`\beta `$, which should be observed in cool DA white dwarfs in the range of effective temperatures (11000 to 13000 K), where the variable ZZ Ceti objects are found (Hébrard et al. 2002b). This line satellite has never been observed in stars or laboratory plasmas. Future FUSE observations and ongoing study of the FUV spectrum of dense hydrogen plasmas may allow this feature to be detected.

###### Acknowledgements.

This work is based on data obtained by the NASA-CNES-CSA FUSE mission operated by the Johns Hopkins University. Financial support to U. S. participants has been provided by NASA contract NAS5-32985. French participants are supported by CNES. The molecular potential computation code was kindly provided by J. Peek; we would like to thank him.


format: text

text:
# Transition between Quantum States in a Parallel-Coupled Double-Quantum-Dot

## Abstract

Strong electron and spin correlations in a double-quantum-dot (DQD) can give rise to different quantum states. We observe a continuous transition from a Kondo state exhibiting a single-peak Kondo resonance to another exhibiting a double peak by increasing the inter-dot-coupling (t) in a parallel-coupled DQD. The transition into the double-peak state provides evidence for spin-entanglement between the excess-electron on each dot. Toward the transition, the peak splitting merges and becomes substantially smaller than t because of strong Coulomb effects. Our device tunability bodes well for future quantum computation applications.

The double-quantum-dot (DQD) is emerging as a versatile system for studying a variety of strongly correlated behaviors 1 ; 2 ; 3 ; 4 ; 5 ; 6 ; 7 . Following the experimental demonstration of the Kondo impurity-spin screening effect in single quantum dots 8 ; 9 ; 10 ; 11 ; 12 ; 13 , recent theoretical investigations of the coupled-DQD system is uncovering new correlated behaviors 1 ; 2 ; 3 ; 4 ; 5 ; 6 ; 7 . These works suggest that the DQD enables a realization of the two-impurity Kondo problem first discussed in the context of metallic systems 1 ; 14 ; 15 ; 16  in which a competition between Kondo correlations and antiferromagnetic (AF) impurity-spin correlation leads to a quantum critical phenomenon. In a different regime of parameters, a related quantum critical phenomenon can occur driven by a competition between intra-dot Kondo coupling to leads and the inter-dot-coupling 2 ; 3 . In each scenario, a transition is predicted to occur between a quantum state characterized by a single-peaked Kondo resonance, and a different state with a double-peaked resonance. Depending on model and DQD geometry–whether series or parallel coupled–both a continuous or discontinuous 6  behavior in the Kondo peak characteristics have been predicted. The quantum transition in the two-impurity Kondo problem has received wide attention in the theoretical literature in the past two decades, to a large extent because the Kondo to antiferromagnetic transition involves an unusual non-Fermi liquid fixed point. Experimental investigation of this problem thus far has not been reported.

Here we describe transport properties of an artificial molecule formed by two-path, parallel-coupled double-quantum-dots, where the inter-dot-tunnel-coupling, t, can be tuned. In the Kondo regime the differential-conductance, dI/dV, exhibits a single peak centered at zero-bias for t comparable to the lead-coupling induced level broadening. Increasing t by less than 10% resulted in a continuous evolution into a split Kondo resonance. At the same time the conductance at zero-bias exhibits a maximum in the vicinity of the transition. This peak splitting behavior in conjunction with distinct temperature dependences in the different regimes demonstrates a direct observation of an inter-dot-coupling-induced quantum transition. Moreover, on the double peak side the zero-bias conductance becomes suppressed; this suppression represents direct evidence that the localized dot spins are becoming entangled into a spin singlet.

While a double-peaked, coherent Kondo effect was observed by Jeong et al. in a series-coupled DQD geometry17 , the existence of a single-peaked Kondo effect at finite inter-dot-coupling has not previously been established. The series-geometry is unsuitable for an investigation of the transition because even a slight decrease in t can drastically reduce the tunnel current below detection. Therefore neither a single-peaked behavior nor a transition could be observed. Furthermore, direct evidence for spin-entanglement was not obtainable. In the parallel geometry the undesirable suppression of the conductance is avoided.

The Kondo effect in a single QD results from the coupling between the dot excess (unpaired) spin and the spin of the conduction electrons in the leads. The energy scale is given by the Kondo temperature $`T_K\sqrt{U\mathrm{\Gamma }}exp[\pi |\mu ϵ_0|(U+ϵ_0)/\mathrm{\Gamma }U]`$ where U is the on-site charging energy, $`ϵ_0`$ is the energy of the single-particle level, $`\mathrm{\Gamma }`$ reflects the dot level broadening from coupling to the leads, and $`\mu `$ is the chemical potential. The fully symmetric DQD system contains the additional inter-dot-coupling parameter, t. The two magnetic impurities, realized by a single excess-spin on each dot, interact through an effective AF coupling, $`J=4t^2/U`$. The new energy scales, t and J, together with U, $`ϵ_0`$, and $`\mathrm{\Gamma }`$, introduce a rich variety of correlated physics. When U is the dominant energy scale as in the case of our quantum dots, for energies below U two different scenarios are predicted for dots coupled in series 1 ; 2 ; 3 ; 4 ; 5 ; 6 . When t$`<\mathrm{\Gamma }`$, the system can be cast 1  into the two-impurity Kondo problem initially discussed by Jones et al. 14 ; 15 ; 16 . The competition between Kondo effect and antiferromagnetism appears as a continuous phase transition (or crossover) at a critical value of the coupling $`J/T_K2.5`$. If however, $`t`$ is tuned to $`t>\mathrm{\Gamma }`$ before the AF transition point can be reached, the system undergoes a continuous transition from a separate Kondo state of individual spins on each dot (atomic-like) to a coherent bonding-antibonding superposition of the many-body Kondo states of the dots (molecular-like) 2 ; 3 ; 7 . Both the AF state and coherent bonding state exhibit a double-peaked Kondo resonance in the differential conductance versus source-drain bias and involve entanglement of the dot spins into a spin-singlet. Therefore they are likely closely related to one another. The parallel-coupled case has only recently been analyzed in a model without inter-dot-tunnel-coupling, t, where the AF coupling occurs via electrostatic coupling. A discontinuous jump in the differential conductance is predicted when the AF state becomes favored 6 .

Our device was fabricated on $`GaAs/Al_xGa_{1x}As`$ heterostructure containing a two- dimensional (2D) electron gas 80nm below the surface, with electron density and mobility of $`n=3.8\times 10^{11}cm^2`$ and $`\mu =9\times 10^5cm^2/Vs`$, respectively. The lithographic size of each dot is $`170nm\times 200nm`$ (see Fig. 1(a)). The eight separate metallic gates are configured and operated with five independently tunable gate voltages. The dark regions surrounding (and underneath) gates 5 represent 120nm thick over-exposed PMMA, which serve as spacer layers to decrease the local capacitance thus preventing depletion, enabling each lead to simultaneously connect to both dots. The experiment carried out at a lattice-temperature of 30mK. Standard, separate characterization of each dot in the closed dot regime 17  yielded a charging energy U of 2.517meV (2.95meV) for the left (right) dot, with corresponding dot-environment capacitance $`C_\mathrm{\Sigma }63.6af(54.3af)`$, and level spacing $`\mathrm{\Delta }E219\mu eV(308\mu eV)`$. Modeling the dot as a metal disk embedded in a dielectric resulted in disk of radius $`r_e70nm(60nm)`$ with 58 (43) electrons. To characterize the DQD and demonstrate the full tunability of our device, in Fig. 1(b) we show the Coulomb blockade (CB) charging diagram of the conductance versus plunger gates V2 and V4 18  as t was increased, for weak coupling to the leads where Kondo correlation is unimportant. The central pincher gate V5 controlled t where a reduction of the gate voltage decreased t. At weak coupling (t small), the electrons separately tunnel through the two nearly independent dots forming grid like pattern, yielding rectangular domains in Fig. 1(b)(1). With increasing t charge quantization in individual dots is gradually lost as the domain vertices separated and the rectangles deformed into rounded hexagon. At large t when separate charge quantization is fully relaxed, the two dots merge into one and the domain boundaries become straight lines (Fig. 1(b)(6)). The evolution of the conductance pattern demonstrates the tunability of our DQD from ionic to covalent bonding states and bodes well for quantum computation applications 19 .

The parameters $`\mathrm{\Gamma }`$ and t govern the delicate DQD Kondo physics. To obtain the necessary large t, the center pincher gate V5 was set so that the zigzag pattern in the charging diagram is barely visible (Fig. 2(a)), ensuring that the Kondo valleys can be located. To accomplish the formation of Kondo states in both dots, pincher gates V1 and V3 were tuned to give a sizable $`\mathrm{\Gamma }\mathrm{\Delta }E`$ to ensure strong Kondo correlation. An estimate for t is based on the fact that the charging diagram indicates a configuration close to the limit of a merged, single large dot (Fig. 2(a), so that the level broadening $`\pi |t|^2/\mathrm{\Delta }E`$ should be comparable to the level spacing $`\mathrm{\Delta }E`$, yielding $`t150\mu eV`$, while $`\mathrm{\Gamma }`$ is deduced to be $`150\mu eV`$ from the half-width of the Coulomb blockade peaks.

Care is required when changing t. The mutual capacitive coupling between the gates and dots gives rise to a complex capacitance matrix. Operating V5 to tune t simultaneously affected the charge on the dots and other gates. Furthermore, slight residual asymmetry in the coupling to the leads caused the maximum of the Kondo zero-bias anomaly (ZBA) in the dI/dV to occasionally occur at nonzero voltages $`V_{SD}0`$ (termed the anomalous Kondo effect 20 ). By adjusting a combination of the remaining gates we tuned the ZBA to be nearly symmetrical about zero-bias. For practical purposes we relied on plunger-gates V2 or V4, which were found to follow V5 in an approximately linear manner. This procedure can be thought as experimentally diagonalizing the capacitance matrix. The desired configuration of an unpaired excess-spin on each dot could be maintained within a tuning range of $`47mV`$ in V5 without causing a sudden change in the charge configuration. Note that V5 is set typically $`70mV`$ above pinch-off. Therefore the relative tuning range in this Kondo regime is roughly 6-10% of closure and even smaller for the corresponding fractional change in t, since tunnel-coupling becomes exponentially suppressed in the small value limit.

Three distinct spin configurations may appear assuming even-odd electron filling and focusing on the topmost states in each dot (Fig. 2(b)). Our investigation was carried out in such regimes. Non even-odd behavior was also observable, but will not be discussed 21 . The Kondo valleys and spin states in each dot were identified by measuring the differential conductance, dI/dV, versus source-drain bias, $`V_{SD}`$, at an electronic temperature of $`40mK`$. For regions 1-6 in Fig. 2(a), we observed the expected appearance and disappearance of Kondo resonance peaks near zero-bias shown in Fig. 2(b).

When t was tuned two distinct regimes of behavior were evident in dI/dV (see Fig. 3(a)-(d)). The main features in the first regime (see Fig. 3(a),(d) bottom half) were the clear presence of a single peak in the ZBA, an increasing peak width and an increasing linear conductance $`g=dI/dV|_{V=0}`$ with increasing t. In the second regime where t was increased further (see Fig. 3(b), (d) upper half), the single ZBA peak developed into two peaks. In some Kondo valleys, both types of behavior were observed as shown in Fig. 3(c),(d), in which the transition is seen to take place in a continuous manner. In the transition region, the ZBA peak broadened and its zero-bias value, g, approached a maximum, became flat before dropping as the ZBA peak gradually split into two. The suppression of g on the double-peaked side can be attributed to the emergence of spin-singlet correlation between the two dot spins 1 ; 2 ; 3 ; 4 ; 5 ; 6 ; 7 . The maximum g is about $`0.1e^2/h`$, reduced below the theoretical unitary limit of $`4e^2/h`$. This reduction can occur when the energy levels in the two dots differ and the dot-lead coupling is asymmetric 3 ; 22 .

Analyzing the transition with detailed curve fitting we quantify the amplitude, splitting and width of the ZBA peak(s) at each V5 gate voltage within well-defined bounds. In the absence of precise theoretical functional forms either for the conductance background or the ZBA peak itself, we first systematically subtracted three sensible backgrounds. After subtraction of each background we fitted the ZBA using three (combined) peak shapes: one (or two) Gaussian, Lorentzian and the Breit-Wigner resonance form. The resultant parameters invariably exhibited similar trends as shown in Fig. 3(e). The double-peak feature visibly disappeared at $`V50.594V`$ as indicated by the dash line, and roughly coincides with the g maximum position. Note that the peak splitting,$`\delta `$, dramatically reduces from a maximum of $`120\mu eV`$ to $`0`$ when V5 is slightly reduced from -0.5925 V to $``$ -0.595 and correspondingly t is reduced by less than 4% from an initial value $`t150\mu eV`$.

In addition to differences in the peak shape the temperature dependence of dI/dV was also distinct in the two regimes (see Fig. 4). In the single peak regime, the zero-bias dI/dV, g(T), decreased logarithmically with T (Fig. 4(a)). In contrast, when double peaks appeared, g(T) exhibited a non-monotonic behavior (Fig. 4(b)), where with increasing T g(T) slightly increased initially, then slowly decreased before increasing again when T exceeded $`T_K`$. Together, these evidences point to qualitatively different phases in the two regimes and the existence of a quantum transition between them.

To date no theoretical work has addressed the parallel-coupled DQD with inter-dot tunnel coupling. Nevertheless, because the ZBA occurs under slightly non-equilibrium conditions it is likely that we may identify the observed transition with the quantum critical phenomenon discussed in the two-impurity Kondo problem for the series geometry based on the following evidence: (a) a continuous evolution from the single- to double peaked behavior, (b) a maximum in g, the zero-bias dI/dV, near the transition point, (c) different behaviors in the temperature dependence of g, and (d) a strong renormalization of the peak splitting, $`\delta `$, compared to the estimated t or AF coupling, J, close to the transition. These features are all in qualitative agreement with predictions for a series-coupled DQD 1 ; 2 ; 3 ; 4 ; 5 ; 6 ; 7 ; 22 , although ideally theory predicts a maximum g reaching the unitarity limit $`4e^2/h`$. Semi-quantitively, it is informative to roughly estimate key parameters and compared these to observed splitting, $`\delta 120\mu eV`$. Within the theoretical scenarios, $`\delta `$ must be compared to $`4t600\mu eV`$ or $`2J=8t^2/U180\mu eV`$. (Note that in the open dot regime U is expected to be reduced from its closed dot value by roughly 1/3 11 , yielding $`U1meV`$.) The reduction of $`\delta `$ compared to 4t and 2J are in agreement with theory and lends further credibility to our identification of the quantum transition.

Acknowledgements We thank H. Baranger, K. Matveev, S. Khlebnikov, N. Giordano, N. Windgreen, D. Cox, A. Schiller, and H. Nakanishi for discussions, and F. Altomare, L.C. Tung, V. Gusiatnikov, and H. Jeong for assistance in the experiment. Work supported by NSF grants DMR-9801760 and DMR-0135931.


text:
# From Kaluza-Klein theory to Campbell-Magaard theorem and beyond

## I Kaluza-Klein Theory

The idea that our spacetime might have five dimensions and that some observable physical effects could be attributed to the existence of a fifth dimension has perhaps no better illustration than that provided by the Kaluza-Klein theory . Indeed, in this theory a scheme is devised in which the extra dimensionality of space is combined with curvature in such a way that electromagnetic phenomena may be looked upon as a pure manifestation of geometry. Although no new physical effect was to be predicted by the theory, the very possibility of unification (even at a mathematical level) between two different subjects such as gravity and electromagnetism has had an appeal which endures until today. Of course along the history of physics there have been many attempts to formulate a unified field theory which would not resort to extra dimensions. In fact, the fifth dimension assumption was basically the cause of Einstein’s reluctance to accept the plausibility of Kaluza-Klein theory . Alternative ways to geometrize the electromagnetic field in the usual four-dimensional spacetime may be found, for instance, in the work of Weyl, who proposed a generalization of Riemannian geometry . More recently, Weyl’s approach has been revived by Novello through the WIST (Weyl Integrable Space-Time) program, where interesting applications to Cosmology have been discussed . The development of particle physics, on the other hand, led to a ressurgence of interest in higher-dimensional field theories, mainly as a possibility of unifying the long-range and short-range interactions. Indeed, inspired by the old five-dimensional Kaluza-Klein theory there appeared, around the sixties and seventies of the last century, higher-dimensional theories such as eleven-dimensional supergravity and ten-dimensional superstrings, all aiming at a unifying scheme .

## II Induced Matter theory

The original version of Kaluza-Klein theory assumes as a postulate that the fifth dimension is compact (condition of cilindricity). Recently, however, a non-compactified approach to Kaluza-Klein gravity, known as Induced-Matter theory (IMT) has been proposed by Wesson and colaborators . The basic principle of the IMT approach is that all classical physical quantities, such as matter density and pressure, should be given a geometrical interpretation . Moreover, it is asserted that only one extra dimension should be sufficient to explain all the phenomenological properties of matter. More specifically, IMT proposes that the classical energy-momentum tensor, which enters the right-hand side of the Einstein equations could be, in principle, generated by a pure geometrical means. The theory also assumes that the fundamental five-dimensional space $`M^5`$, in which our usual spacetime is embedded, should be a solution of the five-dimensional vacuum Einstein equations

$`R_{\alpha \beta }=0.`$

We shall not go into the details of the mathematical mechanism of Induced-Matter theory (the interested reader is referred to  and the references therein). For our purposes here, suffice it to say that IMT has a mathematical structure which can be better understood if it is regarded as a spacetime embedding theory . In this context our spacetime would correspond to a four-dimensional hypersurface locally and isometrically embedded in a five-dimensional Ricci-flat manifold.

## III Campbell-Magaard theorem

Embedding theories are naturally subject to embedding theorems of differential geometry. The claim that any energy-momentum can be generated by an embedding mechanism may be translated in geometrical language as saying that any semi-Riemannian four-dimensional manifold is embeddable into a five-dimensional Ricci-flat manifold. At the time IMT was proposed it was not at all apparent whether such proposition was true. Fortunately (at least for the supporters of the theory), this is essentially the content of a little known but powerful theorem due to Campbell  and Magaard , which asserts that any semi-Riemannian n-dimensional analytic manifold can be locally and isometrically embedded in a semi-Riemannian (n+1)-dimensional analytic manifold, where the Ricci tensor of the latter vanishes . Campbell-Magaard’s result has then acquired fundamental relevance for granting the mathematical consistency of five-dimensional non-compactified Kaluza-Klein gravity.

Local isometric embeddings of Riemannian manifolds have long been studied in differential geometry. Of particular interest is a well known theorem (Janet-Cartan)  which states that if the embedding space is flat, then the minimum number of extra dimensions needed to analytically embed a Riemannian manifold is $`d`$ , with $`0dn(n1)/2`$. The novelty brought by Campbell-Magaard theorem is that the number of extra dimensions falls drastically to $`d=1`$ when the embedding manifold is allowed to be Ricci-flat (instead of Riemann-flat).

## IV  New embedding theories and the need for more general theorems

The increasing attention given to the Randall-Sundrum models  in which the embedding manifold, i.e., the bulk, corresponds to an Einstein space, rather than to a Ricci-flat one, has raised the question whether Campbell-Magaard could be generalized and what sort of generalization could be done. It was conjectured that if the Ricci-flatness condition were replaced by the requirement of the embedding space being an Einstein space, then a result similar to Campbell-Magaard theorem would hold. That this conjecture is in fact a theorem was recently shown by Dahia and Romero . Specific classes of embeddings, such as those of Einstein spaces into Einstein spaces were established . This was the first extension of Campbell-Magaard theorem and more were to come.

A following question to address was whether further extensions of Campbell-Magaard theorem were possible. This led to an investigation of embeddings in spaces that are sourced by dynamical matter fields. One of the simplest forms of matter is that of a scalar field, so it was natural to consider this case first. It was then proved that any n-dimensional analytic Lorentzian or Riemannian space can be locally and isometrically embedded in a (n+1)-dimensional analytic manifold generated by any arbitrary self-interacting scalar field. As a corollary one can prove that any Lorentzian or Riemannian n-dimensional analytic manifold can be embedded in a (n+1)-dimensional space which is a solution of the vacuum Brans-Dicke field equations .

In seeking higher levels of generalization one is led to consider the more general situation of embedding spaces whose Ricci tensor is arbitrary. Pursuing this idea a bit further a third extension of the Campbell-Magaard was finally established. A theorem was proved which asserts that any n-dimensional semi-Riemannian analytic manifold can be locally embedded in a (n+1)-dimensional analytic manifold with a non-degenerate Ricci tensor which is equal, up to a local analytic diffeomorphism, to the Ricci tensor of an arbitrary specified space. As an application of this theorem, embeddings of Ricci-flat spacetimes into five-dimensional Friedmann-Robertson-Walker models were obtained in ref. .

To conclude we would like to point out that insofar as five-dimensional embedding theories are metric it appears to be of relevance to allow the embedding spaces to have different geometrical properties, which must ultimately be determined by the dynamics of the theory in question.

Therefore, generalizations of the known embedding theorems might reveal crucial in building new higher dimensional models. On the other hand, from the standpoint of epistemology, it is rather illustrative to see how modern theoretical physics can be a source of new ideas in mathematics and geometry.


text:
# Excitonic and Quasiparticle Life Time Effects on Silicon Electron Energy Loss Spectrum from First Principles

## Abstract

The quasiparticle decays due to electron-electron interaction in silicon are studied by means of first-principles all-electron GW approximation. The spectral function as well as the dominant relaxation mechanisms giving rise to the finite life time of quasiparticles are analyzed. It is then shown that these life times and quasiparticle energies can be used to compute the complex dielectric function including many-body effects without resorting to empirical broadening to mimic the decay of excited states. This method is applied for the computation of the electron energy loss spectrum of silicon. The location and line shape of the plasmon peak are discussed in detail.

One of the most challenging condensed matter problems is the prediction of excited states of materials from first principlesonida . The solution of this formidable task is doubly rewarding since (1) most of the interesting physics involves the interaction of an electromagnetic radiation with matter, and (2) most of progress in today’s nanoelectronic technologies requires the knowledge of excited states.

In the last few decades, the development of this field has been expending rapidly, and several evolutionary periods can be distinguished. In the first one, the density functional theory in the local density approximation (LDA) played a major role, and was widely used for the analysis of spectroscopic properties of materialslocalfield ; oneelectron , despite that the band gaps of semiconductors and insulators were adjusted by ad hoc methods. Later, it was realized that the GW approximation of Hedinhedin  provides a practical scheme for the determination of the quasiparticle energies, and band gaps were no longer adjustedlouie . Nevertheless, even the use of correct band gaps and the various types of exchange-correlation corrections, including local-field effects, did not improve significantly the calculated optical spectralocalfield . It was often believed that the inclusion of the electron-hole (e-h) interactions was the missing ingredient for an adequate description of the optical spectra, and this was later confirmed by model calculationshankesham . However, it is only recently that ab initio methods explicitly showed the relevance of these complex interactions to the calculation of the dielectric functiononida . At this point, calculations from first-principles were finally making direct contact with experiment and establishing the relevance of the e-h interactions.

The aim of this Letter is to add a further step to this fascinating development by incorporating, from first principles and for the first time, the life time of interband transitions into the calculation of the dielectric function. As a consequence an empirical broadening is no longer needed for a successful comparison of theoretical and experimental optical spectra. This work also improves upon existing excited state calculations based on the pseudopotential (PP) approach in conjunction with plasmon-pole (PlP) models for the screening of the Coulomb interactiononida . In those types of calculations not only the various types of matrix elements of operators are not accurate due to the use of pseudowave functions, the PlP approximation also makes it impossible to determine the imaginary part of the self-energy, and hence the spectral functions and life times of quasiparticles remain unaccessible. Indeed, it has been noticed recentlylebegue ; arnaud ; eguiluz  that GWA implementations based on PP methods led to larger and more $`𝐤`$-dependent shifts than those based on all-electron methods. To illustrate our method, we compute the plasmon resonance and the line shape of the electron energy loss spectrum (EELS) of silicon, one of the most used materials in today’s technology.

The most successful approach for describing elementary excitations probed by photoemission experiments uses many-body Green’s function theory. The energies and life times of quasiparticles are mainly determined by the pole structure of the Green’s function $`G`$ or more conveniently by solving the quasiparticle equationlayzer1 . The self-energy operator $`\mathrm{\Sigma }=\mathrm{\Sigma }^{(1)}+i\mathrm{\Sigma }^{(2)}`$ of the quasiparticles is non-local, frequency-dependant, and generally non hermitian. The non-hermitian part $`i\mathrm{\Sigma }^{(2)}`$ is related to the dominant relaxation mechanisms (scattering on crystal imperfection, electron-phonon or electron-electron interaction) which give rise to the finite life time of quasiparticles. In the present work, $`\mathrm{\Sigma }`$ is computed within Hedin’s GW frameworkhedin  which includes dynamic polarization in the random-phase approximation (RPA). In such an approach, only the electron-electron interaction contributes to the non-hermitian part of the self-energy, i.e to the instability of single-particle excitation.

The complex QP energies $`ϵ_{n𝐤}^{qp}`$ are solutions of the quasiparticle equation for a state labeled $`n𝐤`$, and are in practice obtained using a first order perturbation theorylouie . The real part of $`ϵ_{n𝐤}^{qp}`$ corresponds to the QP energy, and the life time of a single-particle excitation, given by the inverse of the full-width at half-maximum (FWHM) of the QP peak in the spectral function, is defined by $`\tau _{n𝐤}=\left[2\times |\mathrm{\Gamma }_{n𝐤}|\right]^1`$ where $`\mathrm{\Gamma }_{n𝐤}=\text{I}m[ϵ_{n𝐤}^{qp}]`$. The determination of both QP energies and their life times requires the computation of the self-energy matrix elementscomment .

The LDA eigenvalues and eigenvectors obtained by means of the all-electron Projector Augmented-Wave methodbloechl  are used as a starting point for our GWA calculations. The wave functions are expanded into plane-waves up to 20 Ry and three partial waves of $`s`$, $`p`$ and $`d`$ types are used to describe the correct nodal structure of the wave functions near the nuclei.

In Fig. 1 the spectral function of some selected high symmetry k-points ($`L`$, $`\mathrm{\Gamma }`$, $`X`$) of silicon is displayed. The FWHM of the QP peaks increases almost linearly away from the band gap because more and more decay channels of Auger type become available. Physically speaking, a QP is scattered into LDA empty states whose energies are between the QP energy and the Fermi energy.

The GW self-consistency should affect the QP life time. Indeed, the number of decay channels of Auger type is fixed by the smallest possible excitation energy contained in the screened interaction $`\stackrel{~}{W}(\omega )`$ which is defined by the absolute band gap of silicon. As the band gap increases from 0.44 eV to 1.0 eV when going from the LDA to the GW approximation (see Fig. 2.a), we expect a decrease of the imaginary part of the self-energy caused by self-consistent effects. Therefore the QP peak should be narrowed, corresponding to increased life times. Another consequence of the lack of self-consistency is that some states near the Fermi level have a spurious imaginary part because they fall between the LDA and GW Auger thresholds. Thus, the life times of the QPs in the vicinity of the gap are not correct except at the uppermost valence state and at the lowest conduction states where the non-hermitian part of the self-energy due to phonons is strictly zerolautenschlager . A proper description of the QP life time in this range of energy should include the interaction of electronic states with lattice vibrations because the dominant relaxation process is governed by electron-phonon interaction. We expect, however, that our predictions are quantitatively correct at high energy where the dominant relaxation mechanism is electron-electron interaction. Therefore, the complex QP energies obtained by means of the GW approximation represent a good starting point for the direct calculation of the complex dielectric function at high energy without using any phenomenological life time of excited states.

To include the e-h interaction in the dielectric function we solve the Bethe-Salpeter equation, which describes the propagation of a correlated quasielectron-quasihole pair, and can be turned into an effective two-particle Schrödinger Hamiltonian whose block matrix form is given by

$$H^{\text{exc}}=\left(\begin{array}{cc}R& C\\ C^{}& R^{}\end{array}\right).$$
(1)
Detailed expressions of the different block matrix elements can be found elsewhererohlfing . The diagonal blocks given by $`R`$ and $`R^{}`$ are respectively the resonant part ($`v𝐤c𝐤`$ transitions) and the anti-resonant part ($`c𝐤v𝐤`$ transitions) while the off-diagonal blocks $`C`$ and $`C^{}`$ couple positive and negative frequency transitions. Olevano and Reining have shown that these off-diagonal blocks are crucial for the correct description of the EELS of siliconolevano . Here, we follow a similar approach, but in contrast to their work, we include the life time of non-interacting e-h pairs by considering complex transition energies whose imaginary part is set by $`\text{Im}[ϵ_{c𝐤}^{qp}ϵ_{v𝐤}^{qp}]=\mathrm{\Gamma }_{c𝐤}\mathrm{\Gamma }_{v𝐤}`$. Therefore, even if one neglects the off-diagonal blocks, the effective Hamiltonian $`H^{\text{exc}}`$ is non-hermitian because the diagonal matrix elements of $`R`$ are complex.

The macroscopic dielectric function is given by

$$ϵ(𝐪=\mathrm{𝟎},\omega )=1\underset{𝐪\mathrm{𝟎}}{lim}\frac{4\pi }{\mathrm{\Omega }}\times \frac{1}{𝐪^2}\times \underset{n_1,n_2,n_3,n_4}{}n_1|e^{i𝐪.𝐫}|n_2\times [H^{\text{exc}}\omega I]_{(n_1,n_2),(n_3,n_4)}^1\times n_4|e^{i𝐪.𝐫}|n_3\times (f_{n_4}f_{n_3}),$$
(2)
with $`f_{n_i}`$ being the occupation number of the Bloch state labeled by $`n_i`$ and with $`(n_i,n_j)=(v𝐤,c𝐤)\text{or}(c𝐤,v𝐤)`$.

To build the excitonic Hamiltonian $`H^{\text{exc}}`$, we first calculate the complex QP energies for the first 16 bands and for the 216 k points sampling the Brillouin zone. Fig. 2.a shows the shifts of the real part of the QP energies with respect to the LDA energies as a function of LDA energies. Interestingly, these shifts are far from being uniform across the BZ with the exception of LDA states ranging from -5 eV to 5 eV. In this range of energy, the valence LDA states are shifted downward by about 0.45 eV while the LDA conduction states are shifted upward by about 0.15 eV. The striking feature is the decrease of the QP energy shifts for the LDA conduction states above 5 eV. This trend was also observed by Fleszar and Hankefleszar\_hanke . In addition, the QP shifts below -5 eV reflect the contraction of the occupied GW bandwidth with respect to that of the LDA. Such a contraction is not observed when using a plasmon-pole model. Fig. 2.b shows the imaginary part of the QP as a function of the QP peak positions. It’s worth mentioning that the imaginary part of the QP of Ref. fleszar\_hanke  are twice as larger as ours. Such a large difference can certainly be ascribed to the value of the broadening parameter $`\delta `$ used in Ref. fleszar\_hanke  to evaluate $`\stackrel{~}{W}(\omega )`$ along the real axis (see Ref. lebegue  for details). The broadening should be chosen as small as possible ($`\delta 0`$) to obtain correct life times. Notice that $`|\text{Im}(ϵ_{n𝐤}^{qp})|`$ is very small for the QP states ranging from -5 eV to 5 eV and increases almost linearly outside this interval. Thus, we expect a non negligible broadening $`|\mathrm{\Gamma }_{c𝐤}\mathrm{\Gamma }_{v𝐤}|`$ of about 0.5 eV for non interacting transition energies in the range of energy where the plasmon resonance occurs.

Fig. 3 shows our calculated EELS ($`\text{Im}[ϵ^1(𝐪=0,\omega )]`$) of siliconinversion . These results are compared with the experimental spectrumstiebling  characterized by a rather symmetric plasmon peak located at 16.7 eV. The dashed line represents the spectrum calculated within the RPA approximation which amounts to retaining just the complex transition energies in the excitonic Hamiltonian. The main peak of the RPA spectrum is shifted by about 1 eV towards higher energy compared to experiment. In addition, the FWHM of the plasmon peak $`\mathrm{\Delta }\text{E}_{\frac{1}{2}}`$ is underestimated by a factor of 3 while the height is overestimated by the same factor. When including both e-h attraction and local-field effects but still neglecting the coupling between positive and negative frequency transitions (open circles), the line shape of the spectrum worsens. The plasmon-peak position is now ill-defined because the spectrum exhibits two structures located at 15 and 18 eV. Finally, the inclusion of the coupling between forward and backward going quasielectron-quasihole pairs by taking into account the off-diagonal blocks of $`H^{\text{exc}}`$ has a rather drastic effect on the overall spectrum (full circles). On the one hand, the natural line shape is improved and, as expected, $`\mathrm{\Delta }\text{E}_{\frac{1}{2}}3\text{eV}`$ is underestimated with respect to experiment ($`3.7\text{eV}`$) because of the analyzer resolution included in the latter. On the other hand, the main peak of our spectrum is unfortunately shifted, by about 1 eV, towards lower energy compared to experiment, and its height is somehow overestimated. However, the inclusion of the analyzer resolution would slightly reduce the calculated height.

The position of the plasmon peak $`\omega _p`$ is determined roughly by the zero of the real part $`ϵ_1(\omega )`$ of the macroscopic dielectric function. The inset in Fig. 3 clearly shows the shift of the zero towards lower energies when the excitonic effects and coupling terms are included. It’s also of interest to notice that the FWHM of the plasmon peak is qualitatively given by $`\mathrm{\Delta }\text{E}_{\frac{1}{2}}2\times ϵ_2(\omega _p)/\frac{ϵ_1(\omega )}{\omega }|_{\omega _p}`$. Thus $`\mathrm{\Delta }\text{E}_{\frac{1}{2}}`$ is controlled by the value of $`ϵ_2`$ at the plasmon resonance energy and by the slope of $`ϵ_1`$ at the same energy as it can be seen in Fig. 3. For example, a detailed comparison of both calculations, where the exchange-correlation effects are taken into account, shows that the inclusion of the coupling yields a decrease in $`ϵ_2`$ as well as an increase in the slope of $`ϵ_1`$, i.e., it leads to a narrowing of the plasmon peak.

Obviously, the puzzling feature of our calculated spectrum is the underestimation of the plasmon-peak energy. It is therefore worth comparing our findings to the only available PP resultsolevano . The PP plasmon peak is located at about 17.3 eV while ours is located at 15.7 eV. These differences arise because we avoid most of the approximations used in Ref. olevano, . In particular, (1) we avoid using PP methods, and any kind of plasmon-pole model to calculate QP energies (such models are expected to fail at these high energies), (2) we don’t use any broadening parameter to describe the life times of excited states, and finally (3) we avoid using first order perturbation theory to include the coupling between positive and negative frequency transitions. This latter approximation is shown to have only a minor effect on the final spectrum olevano . As a consequence of the first two approximations, their plasmon-peak energy position, at the RPA level without excitonic effects, is already shifted by 0.9 eV towards higher energies with respect to ours. A self-consistent GW calculationeguiluz  should improve the agreement with experiment by enlarging the transition energies. However, self-consistency issues are beyond the scope of this work.

In conclusion, we have calculated the QP energies and life times of silicon within the all-electron GW approximation. We have also demonstrated that these QP energies can be used to evaluate the EELS of silicon without using any phenomenological parameter to describe the life times of excited states. Our results confirm that it is crucial to include excitonic effects and coupling between forward and backward going e-h pairs to obtain an improved line shape. However, the underestimation of the plasmon-peak position reflects the need for a self-consistent calculation of the QP energies or the need for an approximation beyond the GW method.

We are grateful to V. Olevano for providing us with some details about his calculations. The supercomputer time has been granted by CINES (Project No. gmg2309) on the IBM SP4.


text:
# References

Reply: In our Letter Wuetal:02  we presented a numerically exact, polynomial-time quantum computer algorithm for simulating the class of general pairing Hamiltonians (with arbitrary coupling constants $`V_{ml}^r`$)

$$H_\mathrm{p}=\underset{m=1}{\overset{N}{}}\frac{1}{2}\epsilon _m\sigma _m^z+\underset{r=\pm }{}\underset{l>m=1}{\overset{N}{}}\frac{1}{2}V_{ml}^r(\sigma _m^x\sigma _l^x+r\sigma _m^y\sigma _l^y),$$

\[see Eq. (1) in Wuetal:02  for definitions\]. The Comment by Dukelsky et al. Dukelsky:comment  reviews an interesting but numerically approximate classical algorithm \[density matrix renormalization group (DMRG)\] which they claim to be applicable for the same problem. It is this difference between numerically approximate and exact algorithms that essentially renders irrelevant the criticism of our Letter Wuetal:02  in the Comment Dukelsky:comment . There is no known classical algorithm for simulating the class of general pairing models that is both efficient (polynomial) and numerically exact. There is such a quantum algorithm: the one proposed in our Letter.

Unfortunately, there is a misleading and confusing statement in the Comment Dukelsky:comment : Dukelsky et al. write in the third paragraph that (all emphasis ours) “It was shown in this context that only by resorting to the exact numerical solution , using the DMRG …”. However, the DMRG is not an exact algorithm. It is an approximation, as in fact stated by the authors in their Comment a few lines after their claim that the DMRG is exact. Indeed, citing Dukelsky and Sierra’s own paper Dukelsky:99  (their Ref. , Ref.  in our Letter): “We show that the DMRG gives an accurate approximation to the exact ground state if the block is taken to be the set of particles while the environment is taken to be the set of holes.” And on p.174 of the same paper they write: “The DMRG is a variational method and in the region under study we expect our results to coincide with the exact ones with a relative error less than $`10^4`$.”

These statements by the authors of the Comment clearly demonstrate that the DMRG method is not a numerically exact solution of the pairing Hamiltonians. Hence the term “exact” in the the above quote from their Comment is not correct, and the basis for their criticism of our Letter is invalidated.

The authors have provided evidence that the DMRG can give a good approximation to the low-lying energy spectrum of the reduced BCS Hamiltonian (with constant pairing interaction $`V_{ml}^rV^r`$), and can do so in $`O(N)`$ steps Dukelsky:99 . We do not dispute this, and made no statement to the contrary in our Letter, since we were only concerned with a numerically exact algorithm. However, in spite of their claim to the contrary in their Comment (“…the DMRG can easily accomadate arbitrary pairing matrix elements.”), there is presently no evidence to support the hope that the quality of the approximation and the efficiency of the algorithm will remain as good for general pairing Hamiltonians (i.e., $`V_{ml}^rV^r`$ for all $`m,l`$). The burden of proof is upon the authors of the Comment. Moreover, even if they could demonstrate this, the basic point remains that theirs is a numerically approximate algorithm while ours is exact. Indeed, as we stated in our Letter, for half-filled states the dimension of the Hilbert space grows as $`N!/[(N/2)!]^2`$ (for even $`N`$), which is super-exponential in $`N`$, and there is no known way to avoid this divergence in a classical algorithm that is numerically exact. In contrast, we have given a quantum computer algorithm that is numerically exact (in the sense that it exactly diagonalizes the pairing Hamiltonian), requiring only $`O(N^4)`$ steps, for arbitrary $`N`$. Any approximate method, including DMRG and projected BCS Dukelsky:99 , must eventually be checked against numerically exact algorithms if and when we have quantum computers.

In the second part of the Comment, the authors consider a pairing model that includes a monopole interaction \[third term in their Eq. (1)\], not included in our pairing Hamiltonian, so it does not apply to our Letter. They point out that if the coupling constants $`V_{ij}^1`$ and $`V_{ij}^2`$ (in their notation) satisfy certain relations then the model is exactly solvable. This is an elegant result of which we were well aware, but that bears no relation to our work, which deals, of course, with the cases for which no analytical solution is possible. However, Dukelsky et al. write Dukelsky:comment : “We believe that most of the physical problems can be modeled with a pairing Hamiltonian within the integrable subset …”, indicating that the cases that are not analytically solvable are somehow typically not physicallyrelevant. This conclusion is invalidated by their own parameter count: the ratio of integrable to non-integrable models is $`[(6N+3)/(2N^2N)]0`$ as $`N\mathrm{}`$, so the class of integrable models occupies a negligibly small fraction of the parameter space in the interesting regime of $`N1`$.

In conclusion, the Comment Dukelsky:comment  discusses approximate classical algorithms and analytically solvable instances of the pairing model, which are certainly important, but are not relevant in the comparison to a general and numerically exact algorithm for this problem, such as ours.

L.-A. Wu, M.S. Byrd, and D.A. Lidar
Department of Chemistry, University of Toronto, Toronto, Ontario M5S 3H6, Canada

PACS numbers: 03.67.Lx,21.60-n,74.20.Fg