text:
# The Outer Atmosphere of the M-type Supergiant 𝛼 Orionis: K i 7699Å Emission

## 1 Introduction

Massive stars may evolve to become red supergiants prior to death by supernova explosion. Since mass lost by the red supergiant may dictate the time and form of the explosion and also the nucleosynthetic contribution of the massive star to the chemical evolution of the Galaxy, studies of mass loss by red supergiants are of great interest. Mass loss rates are not yet calculable from an ab initio theory. Although measurements of the mass loss rates are available from a variety of observations of the gas and dust in the circumstellar winds, uncertainties remain large. One factor that is often a contributor to the uncertainties is an assumption that the stellar surface, circumstellar shell and the wind are spherically symmetric. Observations are beginning to reveal that supergiants do not recognize fully this simplifying assumption.

Direct observational tests of the assumption are most easily executed on the nearest red supergiants. In the case of Betelgeuse ($`\alpha `$ Orionis), the nearest of these stars, departures from spherical symmetry have been extensively noted. At ultraviolet, visible, and near-infrared wavelengths, observations have shown the presence of ‘bright spots’ on or near the stellar surface (cf. Gilliland & Dupree 1996; Buscher et al. 1990; Young et al. 2000). Such spots that change with time are generally attributed to hot rising convective cells, as predicted by Schwarzschild (1975) and recently shown by Freytag (2002) in 3D hydrodynamical simulations (the hot spots appear occasionally over downdrafts), but an alternative origin has been proposed (Gray 2000). Comparison of ultraviolet and radio brightness temperatures led Lim et al. (1998) to conclude that strong temperature inhomogeneities extend to at least several stellar radii above the surface. The measured radio brightness temperature at 6cm of about 1400K corresponding to a height of about 7 stellar radii is in good agreement with the temperature of about 1200K derived by Justtanont et al. (1999) from infrared ground state fine structure lines of Si ii and Fe ii. Temperatures derived from ultraviolet spectra are much higher at these heights implying that cold and warm gas coexist at the same heights above the stellar surface (Lim et al. 1998). Recently, Harper, Brown, & Lim (2001) have constructed a semiempirical model for the inner shell extending from the photosphere to about 40 stellar radii (about 2 arcsec). It is not unlikely that the inhomogeneities detected in ultraviolet and radio observations persist to greater heights and into the stellar wind.

Dust in the shell close to the star has been resolved by infrared interferometry and direct imaging, and detected far from the star by the polarization that scattering of dust grains introduces. Infrared emission at about 11 $`\mu `$m from silicate dust grains comes from an extended region of several arcsec with an inner shell radius of 1 to 2 arcsec (cf. Sloan, Grasdalen, & LeVan 1993; Danchi et al. 1994; Bester et al. 1996; Rinehart, Hayward & Houck 1998; Sudol et al. 1999). Undoubtedly, dust - cooler grains - extends to greater distances but has yet to be resolved by far-IR imaging. Extended emission in IRAS images was noted by Stencel, Pesce, & Hagen (1988) and Young, Phillips, & Knapp (1993a, 1993b). Warm and cold dust is observable by the photospheric light that grains scatter. This scattered light is separable from light scattered by the Earth’s atmosphere and the telescope by its polarization: the polarization vector is tangential to the radius vector. McMillan and Tapia (1978) following a suggestion by Jura and Jacoby (1976) detected the polarized scattered light from Betelgeuse’s shell. Additional measurements were made by Le Borgne, Mauron, and Leroy (1986) who show that the number density of dust grains declines roughly as the inverse-square of the radial distance from the star for distances of 10 to 90 arcsec.

Circumstellar gas as detected by blue-shifted absorption lines was observed and analyzed by Weymann (1962) who provided indirect evidence that the shell extended to great distances from the star. Direct evidence in the case of the visual binary $`\alpha ^1`$ Herculis for a very extended shell off this M supergiant was presented by Deutsch (1956) who found circumstellar lines in the spectrum of the visual companion to the supergiant. Weymann’s inferences about the shell size were compatible with the direct lower limit to the outer shell radius from Deutsch’s observation. Observations of H i 21 cm emission from Betelgeuse have shown that gas extends to about 1 minute of arc from the star (Bowers & Knapp 1987).

Betelgeuse’s circumstellar shell was resolved spatially by Bernat & Lambert (1975) who detected fluorescent emission in the K I 7699Å resonance line at angular distances of 2 - 4 arcsec from the star. Subsequent observations extended the detections of the 7699Å emission to about 60 arcsec from the star and provided evidence for some departures from spherical symmetry and for structures in the shell (Bernat & Lambert 1976; Bernat et al. 1978; Honeycutt et al. 1980; Mauron et al. 1984; Mauron 1990). Emission in the Na D lines was detected by Mauron & Querci (1990), Mauron (1990), and Mauron & Guilain (1995). The most remote (from the star) manifestation of the circumstellar shell is the parsec-sized bow shock reported by Noriega-Crespo et al. (1997) from IRAS images at 60 and 100 $`\mu `$m. Emission detected from an arc with a mean radius of 6 minutes of arc is attributed to material confined by the interstellar medium’s ram pressure.

In this paper, we present and analyse a new series of observations of K i 7699Å emission obtained with a long slit at a resolving power of about 110,000 or 2.6 km s-1 with the slit placed at various impact parameters from the star and position angles around the star, and at several epochs. Our combination of spectral and spatial resolution was not achieved by earlier studies of the K i emission. We derive afresh the radial dependence of the K i emission and rediscuss the principal conclusion - the mass loss rate - that has been previously drawn from the radial dependence. A novel result of our observations is the discovery of clumps in the emission and, hence, in the neutral gas in the shell and wind. Some of these inhomogeneities are unresolved both spatially in our seeing-limited observations and also kinematically at the limit of about 2 km s-1.

## 2 Observations and Reductions

Betelgeuse was observed in 1994 February, March, and November with the W.J. McDonald Observatory’s Harlan J. Smith 2.7m reflector and its coudé spectrograph. The latter was used with an echelle grating and an interference filter to separate the order providing the K i 7699Å line. The slit length was varied with a majority of the observations taken with a length of 10 to 100 arcsec. The slit width was almost always set at 0.55 arcsec. Spectra were recorded by a Tektronix 512 X 512 CCD over a bandpass of about 17Å. Standard reduction procedures were applied to correct for the bias and pixel-to-pixel variations. The wavelength scale was set using Th-Ar lamp spectra.

During each observing session, several on-star spectra were recorded with a very short slit length. Off-star spectra were then collected with longer slit lengths. In general, the star was not included in the slit; inclusion limited the exposure time beyond which the starlight saturated a portion of the CCD chip. Since the field at the coudé focus rotates, exposure times were kept short to minimize smearing of emission inhomogeneities; typically, exposure times were 10 minutes or less. For capturing emission far from the star, a few exposures as long as 30 minutes were made. Relative spectrophotometry was attempted; observations were terminated in partly cloudy skies. The stellar flux was measured by opening the slit to its maximum width when all the visible image was contained within the slit. A guider maintained the star’s position relative to the slit to an accuracy of about 0.5 arcsec.

Our reduction techniques were adapted directly from those described earlier (Plez & Lambert 1994). In short, the spectrum at a given location off the star is presumed to be a composite of emission from the gas down the line of sight through the shell at that point and starlight scatter in the Earth’s atmosphere and the telescope. We assume that the latter contribution may be represented by a scaled version of an on-star spectrum. Fig. 1 provides an example of a raw off-star spectrum and the spectral image after subtraction of the photospheric scattered light.

For positions of the slit near the star (a few arcsec), the photospheric scattered contribution largely dominates over the light scattered by the shell. This makes the subtraction difficult, as the imperfectly removed background spectrum shows in Fig. 1. The terrestrial O2 lines appear weakly in emission in the image after subtraction. This is due to the difference in illumination across the slit when the star is on one side of the slit. For spectra taken very near the star, we could improve subtraction of the photospheric component by shifting the spectra by less than half a pixel in the dispersion direction. Slight changes in airmass between the on-star and off-star exposures may also affect the subtraction of the O2 lines.

The spectrum after subtraction of the scattered light is a velocity-position map of the emission. Rows of pixels parallel to the dispersion give the spectrum at different points in the shell. The radial distance of a point from Betelgeuse may be calculated from the impact parameter (3 arcsec in this case) and the distance in seconds of arc from the point on the slit of closest approach to the star. The orientation of the slit was derived from the hour angle of the observation. A thin uniform shell of outwardly flowing and fluorescing K i atoms appears as an ellipse in such a velocity-position map. Such an elliptical looking feature is evident in Fig. 1. The extent of the ellipse along the slit is set by the radius of the shell and the impact parameter. The maximum width of the ellipse is directly related to the expansion velocity and impact parameter; the width vanishes when the impact parameter equals the shell’s radius.

Several assumptions underlie our reductions. The recorded off-star spectrum results from a convolution of the light from points around a position on the slit with the seeing profile. A deconvolution of the observed spectra to obtain the true spectrum at a point is impossible given the varying seeing conditions. We assume that the dominant contaminant is scattered starlight with a spectrum represented by the on-star spectrum. The on-star spectrum is itself, however, not devoid of signatures of circumstellar K i absorption and emission at 7699Å. In particular, there is strong blue-shifted absorption introduced by the potassium atoms along the line of sight to the star. At some level, this absorption is weakened by emission from the immediate environs of the star. Certainly, the angular diameter of the stellar photosphere of about 0.05 arcsec is considerably smaller than the smallest aperture (0.55 x 0.55 arcsec) used for the on-star spectra. There is no obvious way to retrieve the true stellar spectrum from the data. A predicted photospheric profile could be used but this would not account for the fact that much of the shell absorption probably occurs close to the star. In view of the low inferred intensity (relative to the stellar scattered light) and the steep gradient of the emitted light from the shell, we suggest that our simple reduction technique is adequate.

Our assumption that the circumstellar shell’s contribution is solely emission in the K i line ignores a contribution from the dust grains that scatter photospheric light. This scattered light from different locations in the shell will be Doppler-shifted by different amounts off the (assumed) radially-flowing dust grains, i.e., the scattered light will be a smeared version of the photospheric spectrum. Optical depth due to dust is small and multiple scattering of photons may be neglected. Circumstellar gas is flowing out at a velocity of about 10 km s-1. If the dust and gas are well coupled, the Doppler shifts imposed on light scattered by dust grains will amount to a line broadening of less than about $`\pm `$10 km s-1. Since the intrinsic widths of the photospheric lines exceed this value, it will be difficult to distinguish the dust-scattered contribution from that scattered by our atmosphere and telescope. Such a distinction is especially hard to achieve with our spectra in which the K i line is the sole strong line. Perhaps, the isolation of the dust-scattered component can be made using a region richer in deep absorption lines. It may be necessary to exploit the polarization of the scattered light to effect a clear isolation of the dust-scattered light. We assume that the dust makes a negligible contribution to the total scattered light.

Our interpretations are based on the assumption that the K i emission results from fluorescence and not from radiative recombination or collisional excitation. Calculations using densities and temperatures expected of the stellar wind fully support this assumption.

## 3 An Asymmetric and Clumpy Shell

Inspection of the reduced spectra shows that the circumstellar envelope as traced by the K i emission about the star is not completely spherically symmetric. A particularly striking feature of the velocity-position maps is the appearance of clumps at scales extending down to the spatial and spectral resolution of our observations. This degree of clumpiness has not been revealed by previous observations. Sample velocity-position maps are shown in Fig. 2, 3 and 4 for slit positions at increasingly larger distances from the star illustrate these departures from spherical symmetry which we discuss in detail below. These Figures also show that the maximum velocity extent of a velocity-position map is not greatly different for small and large impact parameters.

### 3.1 Radial Dependence of the K i Emission and Mass Loss Rate

By azimuthal averaging to the maximum extent possible for a given sequence of observations, it is possible to determine the mean radial dependence of the emission intensity integrated over the line profile. This was done only for exposures obtained under good photometric conditions. The intensity of the scattered light at the radial distance $`r`$ is written as $`I_{scatt}(r)`$. We scale the intensity to the measured stellar flux $`I_{}`$ measured with the very wide slit. $`I_{}`$ was measured by integrating the on-star spectrum flux over the whole spectral interval (512 pixels corresponding to 17.24Å). This total flux was then divided by the wavelength range and the integration time. The resulting $`I_{}`$ is thus expressed in counts per Å per second. Note that Mauron et al. (1984) used a similar definition in averaging the flux over 20Å, while Mauron (1990), Honeycutt et al. (1980), and Bernat et al. (1978) used the photospheric flux in 1Å centered on 7699Å. There are a number of sources of uncertainty affecting the determination of the photospheric flux and in the $`I_{scatt}(r)/I_{}`$ ratio, mostly because of possible loss of photons intercepted by the slit jaws, and changes of transparency and seeing conditions between exposures. The derived normalized intensity of the K i emission, $`I_{scatt}(r)/I_{}`$, is shown in Fig. 5 for impact parameters $`R`$ out to 50 arcsec from the star. The mean profile is consistent with a power law of slope -2.36 $`\pm `$ 0.03 (least-square fit to the data).

Perhaps, one should fit the data with a two component model: a slope of -2.59 $`\pm `$ 0.2 out to 7 arcsec and of 2.55 $`\pm `$ 0.07 beyond 10 arcsec but brighter by a factor of 1.4.

Our measured slope is equal within the uncertainties to that reported earlier by Honeycutt et al. (1980) who gave the slope as -2.65 $`\pm `$ 0.2 from average measurements in the four principal quadrants at radial distances of about 10 to 50 arcsec. Estimates of the intensity of the scattered light are similarly in good agreement: we obtain $`I_{scatt}/I_{}10^{5.0}`$ at 20 arcsec from the star, and Honeycutt et al. reported a value of $`10^{4.8}`$. The agreement is less satisfactory with Mauron et al. (1984) who reported a slope of -3.5 $`\pm `$0.8 from measurements spanning distances of 20 to 55 arcsec but the intensity at 20 arcsec ($`10^{4.8}`$) is consistent with other measurements. Mauron’s (1990 - see also Mauron & Querci 1990) measurements at 5 to 8 arcsec are similarly a little higher in intensity than our mean values. Although the differences could be due to changes in the efficiency of fluorescence, the density of potassium atoms in the shell, and other physical factors, we think it more likely that they reflect observational errors encurred in the surface photometry. As noted by Mauron (1990), earlier results by Bernat et al. (1978) gave much higher intensities for the emission close to the star than are now reported: e.g., Bernat et al. found $`I_{scatt}/I_{}10^{1.3}`$ at about 3 arcsec but our value is 10-3.1. This large difference is unlikely to be due to changing conditions at the stellar surface or in the shell but reflects the difficulty of making photometric measurements of a faint light source close to a bright star with a detector (a vidicon) with a non-linear response.

A calculation of the stellar mass loss rate from the radial dependence of the 7699Å emission has been discussed previously, notably by Rodgers & Glassgold (1991). Although the observations may provide reliable estimates of the run of the density of potassium atoms, an estimate of the mass loss rate calls for knowledge the degree of ionization of potassium, the depletion of potassium onto dust grains, excitation of the atoms (fluorescence vs radiative recombination, for example) in the shell, and radiative transfer through the shell. We assume the shell to be optically thin to the resonance photons, a conclusion that is not contradicted by the ratio of emission in the two resonance lines (Mauron & Querci 1990). Here, we adopt the recipes provided by Rodgers & Glassgold and adjust their derived mass loss rate to our measured intensities.

The slope of the $`I_{scatt}(r)/I_{}`$ versus $`r`$ relation is, if spherical symmetry is assumed, directly related to the radial dependence of the number density of neutral potassium atoms. Additional assumptions include uniform illumination of the atoms by the photosphere, and a shell that is optically thin to the 7699Å photons. It is then readily shown that a relation $`N(r)r^\theta `$ provides an intensity dependence $`I_{scatt}/I_{}r^{(\theta +1)}`$ (Bernat 1976; Gustafsson et al. 1997). If neutral potassium atoms are conserved in a radial outflow, their density dependence follows $`\theta =2`$ and $`I_{scatt}/I_{}`$ follows a power law of slope 3. Conversely, the observed slope of -2.36 indicates a density dependence with a slope of -1.36. Our observed slope is in good agreement with the value of -2.33 corresponding to the model by Rodgers & Glassgold. The difference between their prediction and the slope of -3 predicted for a uniform shell with conservation of potassium atoms is due to their introduction of ionization of potassium atoms in the shell by chromospheric ultraviolet photons. Glassgold & Huggins (1986) and Rodgers & Glassgold (1991) discuss in detail the run of the neutral potassium fraction. The ionizing chromospheric radiation field is not much shielded in the envelope and decreases as $`r^2`$. The electron fraction, set constant in Glassgold & Huggins model, increases slightly outwards in Rodgers & Glassgold’s more detailed model resulting in a gradient of the electron density which is shallower than $`r^2`$. The recombination coefficient increases with distance as a result of the temperature decrease. The neutral fraction thus increases with distance from the star in the inner shell. At large distances out from the star where ionization of potassium atoms is controlled by the interstellar radiation field, the slope of $`I_{scatt}(r)/I_{}`$ must steepen. This steepening is visible beyond about 30 arcsec in our data on Fig. 5, although at first sight masked by the detached shell S2 (see Sec. 3.2.2). The slope increase (to about -4.2) is consistent with Rogers and Glassgold (1991) model.

Following Gustafsson et al. (1997) we write:

$$\frac{I_{scatt}(\beta )}{I_{}}=\frac{\pi e^2}{m_ec^2}\lambda ^210^8f\frac{d^2}{4\pi 206265^2}_{\mathrm{}}^+\mathrm{}N_{scatt}\frac{dz}{r^2}\mathrm{arcsec}^2.\mathrm{\AA },$$

where $`\beta `$ is the angular distance from the star in arcsec, $`f`$ is the oscillator strength of the line, $`N_{scatt}`$ is the number density of scattering atoms in cm-3, $`d`$ is the distance to the star in parsec, and $`z`$ is the distance along the line of sight in cm. Assuming a power law of the form: $`N_{scatt}=N_{s15}(\frac{10^{15}}{r})^\theta `$, with $`N_{s15}`$ the number density at $`r=10^{15}`$cm, the integral is readily calculated and yields: $`N_{s15}(10^{15})^\theta Y(\frac{206265}{\beta d})^{(\theta +1)}`$, with $`Y`$ a constant depending on the value of $`\theta `$. From our observations we derive $`\frac{I_{scatt}(\beta )}{I_{}}=\frac{0.0117\pm 0.0011}{\beta ^{2.36\pm 0.03}}`$. Adopting a distance of 140 pc for Betelgeuse from its Hipparcos parallax of $`7.63\pm 1.64`$ mas, we find the number density of scattering atoms to be $`N`$(K i) = 6.0 $`\times 10^4\times (10^{15}/r)^{1.36}`$ cm-3.111Gray (2000) has remarked that the accuracy of the Hipparcos parallax may be seriously underestimated because it assumes that the centroid of the stellar image remained fixed relative to the stellar disk over the period of observations; a bright spot midway between disk center and the stellar limb and contributing 25% of the light could cause the centroid to shift by half the parallax.

The column density of neutral K atoms along the line of sight to the star was derived by Bernat (1977) from the 7699 Å circumstellar absorption line: $`N`$(K i) $`4\times 10^{12}`$ cm-2. Integrating our relation for the radial dependence, Bernat’s column density is reached by integrating from infinity to an inner radius $`R_{inner}`$ = $`9\times 10^{13}cm`$ $`2R_{star}`$. This is a plausible value for the inner radius. To halve the column density requires $`R_{inner}13R_{star}`$. To double for $`R_{inner}<R_{star}`$, which makes no sense. In short, our derived radial dependence for the density of potassium atoms is consistent with the column density derived from the absorption line.

For this reason and because our observed slope for $`I_{scatt}(r)/I_{}`$ is equal to that predicted by Rodgers & Glassgold, we use the ratio of our observed to their predicted number density versus radial distance directly to correct their mass loss rate. This corresponds to a reduction of their mass loss rate from $`\dot{M}`$ = $`4\times 10^6`$ $`M_{}`$ y-1 to $`2.4\times 10^6`$ $`M_{}`$ y-1. The revision is dependent on our value of $`I_{scatt}/I_{}`$ which rests on an average of the stellar photospheric flux averaged over 17 Å around the potassium line. The stellar flux intercepted by the outflowing gas of the shell is that on the blue wing of the K i line. As already pointed out, we do not have access to the true photospheric spectrum, and the K i line profile can only be guessed at. A rather conservative estimate of a factor of 2 lower flux at line center leads to a mass-loss rate larger by a factor of two. A distance of 200 pc was assumed by Rodgers & Glassgold and Glassgold & Huggins. To correct the mass loss rate to the Hipparcos distance of $`d`$ = 140 pc, we use the scaling given by Mauron & Caux (1992), i.e., the K i brightness at a given angular distance from the star is approximately proportional to $`\dot{M}^2/d`$. Adoption of $`d`$ = 140 pc instead of 200 pc implies a reduction of $`\dot{M}`$ by 85 % or $`\dot{M}=2\times 10^6M_{}`$ y-1. Glassgold & Huggins warn that the accuracy of an $`\dot{M}`$ derived from K i emission is low (“a factor of some 3 - 5”). Although the subsequent analysis of the temperature profile (Rodgers & Glassgold 1991) may have reduced the uncertainty, a factor of two uncertainty seems likely to be optimistic. The derived rate is dependent too on the assumption of a spherically symmetric shell, an assumption probed observationally in the next section.

This uncertainty could in principle be reduced by observing other emission lines. In light of the higher abundance of Na and using the atomic concentrations predicted by Rodgers & Glassgold’s (1991), the Na D emission is expected to be stronger than that in the K i 7699 and 7665 Å lines. Mauron & Guilain (1995) find that the Na D emission is considerably fainter than expected although other mass losing M stars do show Na D emission stronger than the K i emission. These authors suggest that interstellar absorption is “the simplest explanation, though not completely convincing”. Fortunately, interstellar absorption in the K i lines is observed to be considerably weaker than in the Na D lines.

Betelgeuse has been subjected to numerous investigations of its mass loss rate using a variety of spectroscopic indicators. The K i-based rate may be compared with a couple of determinations using tracers of a major constituent of the shell. Bowers & Knapp (1987) detected 21 cm H i emission to derive $`\dot{M}1.1\times 10^6`$ $`M_{}`$ y-1 if $`d`$ = 140 pc. Huggins et al. (1994) detected the C i 609 $`\mu `$m fine structure line and inferred a mass loss rate also of $`\dot{M}1\times 10^6`$ $`M_{}`$ y-1 for $`d`$ = 140 pc. The close agreement between these results and ours is encouraging but must be largely fortuitous. More recently, Harper et al. (2001) estimated $`\dot{M}=3.1\pm 1.3\times 10^6M_{}`$ y-1 from absorption features in ultraviolet Fe ii lines.

### 3.2 Departures from Spherical Symmetry

Azimuthal averaging produces a smooth run of the emission line’s flux with radial distance but masks a fascinating structure of incomplete shells and clumps. This structure is illustrated in Fig. 2, 3, and 4 and additionally in Fig. 6, and 7. Several features of the structure are seen by inspection: the appearance of thin incomplete shells, intense clumps of small spatial extent and velocity dispersion, and lines of sight - particularly far from the star and inside the outermost detected shell - from which no 7699Å emission is detected, and a seemingly sharp boundary to the outermost shell.

The following relations are helpful in discussing the properties of the circumstellar shell. At the assumed distance of 140 pc, 1 second of arc corresponds to a projected distance of 140 AU or a light travel time of about 19 hours. The stellar angular diameter is taken to be 0.044 arcsec (Dyck et al. 1998) or a stellar radius of 3.1 AU. At an expansion velocity of 15 km s-1, gas travels 3.2 AU in 1 year or an angular distance of just 0.02 arcsec.

#### 3.2.1 Temporal Dependence?

Our observations were obtained over several months in 1994. In different observing runs, several successful attempts were made to acquire a series of spectra with the slit set relative to the star at the same positions and orientations. A goal of these observations was to search for changes in the velocity-position maps.

Expected time scales for changes to occur in the wind itself are much longer than a few months (see below). More rapid changes may result from the rise and fall of bright (or dark) spots on the stellar surface. High angular resolution images of Betelgeuse’s disk show a few (1 - 3) bright spots contributing 10 - 20 per cent of the light and changing on a timescale of weeks to months (Wilson, Dhillon, & Haniff 1997). With a light travel time of about a month to a radial distance of about 30", intensity fluctuations of (say) 20 % evolving on a time scale of a month are not going to produce dramatic changes in the fluorescent emission. Were a spectacularly bright spot to occur, a large intensity contrast would be seen in the hemisphere facing the spot. The appearance of shells and especially the presence of clumps cannot be traced to spots on the illuminating photosphere.

Changes in the density of potassium atoms depend on the rates of ionization and recombination. With rates taken from Glassgold & Huggins (1986), we find the timescale for photoionization by the stellar radiation field to be $`t_{photo}0.1\theta ^2`$ yr where $`\theta `$ is the distance from the star expressed in arcsec. The recombination timescale is necessarily much longer: $`t_{rec}30\theta ^2`$ yr (or $`\theta ^{1.3}`$ with the temperature dependence of the recombination coefficient taken into account). Except perhaps very close to the star, these estimates show that changes in K i emission occurring on a time scale of a few months are unlikely to be traceable to variations in the photoionizations by the stellar radiation field.

Fig. 3 shows pairs of velocity-position maps taken at the same slit position relative to the star but at 8 months apart. By inspection, there is considerable similarity between the March and November observations for which the impact parameter of the slit was 13 arcsec. Gross structures may be deemed identical. Smaller structures (here dubbed ‘clumps’) show no discernible change in location. Differences in shape of the smaller structures may be attributed to differences in the seeing: the seeing was about 1.2 arcsec for the March observation but about 3 arcsec in November.

Fig. 6 and 7 show similar pairs of velocity-position maps from the March and November runs but with an impact parameter of 21 arcsec.

The exposure in Fig. 7 was obtained immediately following that in Fig. 6. In these figures, there is evidence for changes in the smaller structures. Indeed, the March pair of maps clearly differ in the inner parts from the same area in the November maps. Close inspection shows a few differences for clumps in successive maps in March and in November. This is not surprising given that successive maps owing to image rotation at the coudé focus do not sample exactly the same portion of the shell. The telescope was guided on the star and the slit rotates about the star. Changes in position on the sky recorded by a particular pixel can exceed the angular size of some clumps, and, hence, these may appear clearly in one image but not in the next. This instrumental effect increases with increasing impact parameter which may explain in part why changes are more evident in Fig. 6 and 7 than in Fig. 3.

#### 3.2.2 A Distant Spherical Shell

A striking feature of the velocity-position maps with impact parameters greater than about 20 arcsec is the appearance of an elliptical shape for the outermost emission (Fig. 4. Emission is very weak, possibly absent, beyond about 1 minute of arc from the star. Ellipticity of isophotes in the velocity-position map corresponds to a spherically symmetric expanding shell. Geometry of the slit position (at 33 arcsec from the star) and the semi-major axis of the ellipse show that the shell radius is 50 $`\pm `$ 2 arcsec. The maximum velocity width, the ellipse’s semiminor axis, corresponds to an expansion (radial) velocity of 18 $`\pm `$ 2km s-1. Note that the shell is clearly seen in Fig. 5 as the plateau followed by an intensity drop at about 50 arcsec. At that distance from the star, the mean K i density is rapidly decreasing with distance, due to the interstellar UV field (see Fig. 5 and Sec. 3.1), and the shell represents a considerable density enhancement.

The thickness of the shell is 1 to 2 arcsec in position on the sky and is unresolved in velocity, i.e. less than 2.6 km s-1. Non-uniform emission over the spherical shell could result from one or more of several factors: shell ejection occurred preferentially over one hemisphere of the star; the excitation of the shell is non-uniform owing to a bright spot on the stellar disk; the ionisation of potassium atoms is greater over one hemisphere due to a chromospheric enhancement; the shell is plowing into the irregular local interstellar medium and deceleration enhances the local density non-uniformly.

Circumstellar absorption imposed upon the photospheric spectrum is composed of two components in atomic (e.g., K i 7699 Å as shown by Bernat & Lambert  and Goldberg et al. ) and the CO lines by Bernat et al. (1979) who labelled the components S1 and S2 where the former has a higher excitation temperature (200K vs 70K), substantially larger CO and K i column densities and lower expansion velocities. The expansion velocity deduced from the S2 absorption lines (20.2 km s-1 from Goldberg et al., and 13 km s-1 from Bernat et al.) are consistent with our measurement. Note that the absorption line measurements are made relative to photospheric lines but the photospheric velocity is variable. Our measurement is a true expansion velocity, relative to the center-of-mass of the star, independent of the stellar velocity, and of photospheric velocity variations.

Another notable difference between the S1 and S2 components is the very narrow line width of the latter: a Doppler width of 1 km s-1 for S2 but 4 km s-1 for S1 from visual lines (Bernat 1977). The Doppler width for S2 corresponds to a full width at half intensity of 1.7 km s-1, an estimate consistent with the upper limit obtained from our observation of the shell emission. Bernat et al. from a comparison of the collisional and radiative rates for CO excitation inferred that the former dominate and, hence, their excitation temperature was approximately equivalent to the gas kinetic temperature. Reference to a dust shell model (Tsuji 1979) showed a gas temperature of 70 K was achieved at about 55 arcsec from the star, a prediction quite consistent with our observed shell radius.

Judged by velocity of expansion, line width, and shell radius as measured or estimated from the CO lines and our 7699Å velocity-position maps, we identify the circumstellar gas contributing the S2 absorption lines with the outermost thin shell seen in Fig. 4. This identification allows interesting further deductions to be made. In particular, the volume density of CO molecules and neutral K atoms may be estimated from their column density and the shell thickness. The equivalent width of the K i 7699Å S2 line estimated from Goldberg et al. (1975, their Figure 1a) gives a column density N(K i) $`1.3\times 10^{11}`$ cm-2 on the assumption that the line is unsaturated. After correction for seeing, the shell’s thickness may be about 1 second of arc. At $`d`$ = 140 pc, a shell thickness of 1 second of arc corresponds to 2 $`\times 10^{15}`$ cm. Then, the K i density in the shell is about 6 $`\times 10^5`$ cm-3. In contrast, Rodgers & Glassgold’s (1991) model predicted a density of $`n`$(K i)$`2\times 10^6`$ cm-3 at the shell’s distance. This difference is not surprising as the shell by inspection of Fig. 4 is evidently a considerable density enhancement above the gas at smaller and greater distances In addition, the intensity variation around the shell is obvious implying that the cut through the shell to the star that gives the K i absorption line may not be representative of mean conditions in the shell.

If the shell at 55 arcsec is identified with Bernat et al.’s S2 shell, their measured column density (1.2 $`\times 10^{16}`$ cm2) translates to a volume density $`n`$(CO)$`6`$ cm-3. The ratio of the CO and K i column densities suggests a density ratio $`n`$(CO)/$`n`$(K i)$`9\times 10^4`$. The elemental abundance ratio in the photosphere is probably a little less than the solar ratio of 2700 on account of a lower C abundance resulting from the first dredge-up. If C has been reduced by about 50%, and K is not depleted onto grains in the wind, the expected elemental ratio in the wind is 1350. The higher ratio in the S2 shell implies, as expected, that K is largely ionized: the fraction of neutral K atoms is $`y_\mathrm{K}0.014`$ provided that most of the C is in CO. Huggins et al.’s (1994) prediction of ionization of C and association into CO through the shell suggests CO accounts for about 20 % of the carbon. Then, $`y_\mathrm{K}`$ 0.003. Rodgers & Glassgold (1991) ‘standard’ model predicts $`y_\mathrm{K}0.01`$ at the radial distance equivalent to 50 arcsec at $`d`$= 140 pc. This is fair agreement with observation given the various (uncertain) ingredients that enter into both estimates. The mass in the shell amounts to about 10-2 M, at a distance of 140 pc and assuming $`y_\mathrm{K}=0.01`$.

As noted above the S2 shell is not a consequence of a varying photospheric illumination; the light travel time across the shell is only about a day, and the shell is only about 50 light days from the star. In addition, the shell seen in absorption in the 1970s (Goldberg et al. 1975; Bernat et al. 1979) is plausibly identified with the K i emission shell seen in 1994. The simplest explanation for the shell may be that it is mass ejected at a substantially higher than average rate at an earlier time. If the expansion velocity has been constant, gas travelling at 18 km s-1 and now at 50 arcsec from the star was ejected 1900 years ago ($`d`$ = 140 pc is assumed). A shell thickness of 1 second of arc implies ejection lasted about 40 yr but this is a guess because the velocity of ejection presumably varied too. The shell as it appears in Fig. 4 is not uniformly bright. It is unlikely that the observed contrast around the shell is due to the presence of bright spot on the disk at the time of the observations. A possible explanation is that mass ejection occurred preferentially above a spot 1900 years ago. Stellar rotation (the rotation period of Betelgeuse is at least a few years), could have led to material being ejected into the shell at a variable rate over the lifetime of the spot. There appear to be at least two other shells of moderate completeness which suggests that significant increases in mass loss may occur every several hundred years. Multiple shells were also detected around $`\mu `$ Cep by Mauron (1997), with a characteristic time scale between two successive increases in mass loss of about one thousand years. Detached shells are also observed around AGB stars, where they are attributed to a Helium shell flash that induces both a mass loss increase and a two-wind interaction due to the increased outflow velocity (Steffen and Schönberner 2000). Clearly, an alternative mechanism must be operating in supergiant stars like Betelgeuse.

Future observations of the K i emission from the shell should be made to search for changes in the S2 shell. Measuring the S2’s expansion will be a challenge: the proper motion expected is a mere 0.03 arcsec per year.

#### 3.2.3 Clumps

A striking feature of our velocity-position maps, especially those acquired in the best seeing is the appearance of clumps or knots of emission; Fig. 3, 6 and 7 show clumps extending in size down to the velocity and spatial resolution of the map. There are locations where clumps appear connected suggesting they are condensations in a sheet or shell. In general, the collection of clumps in many velocity-position maps accounts for more than half of the emitted flux. The impression of elliptical structures is largely gained from gas at large angular distances from the star. Many of our exposures, especially those taken with small impact parameters, were too short to reveal faint distant structures.

These are the first observations of such clumps in Betelgeuse’s wind. Mauron (1997) from long-slit spectra showing K i emission from $`\mu `$ Cephei reports ‘inhomogeneities’ with an apparent size of 1′′ to 3′′ when the angular resolution was about 0.6′′. In velocity space, $`\mu `$ Cep’s clumps were unresolved at the moderate resolution of 40 km s-1 of these observations. Since $`\mu `$ Cep is about 5 times more distant than Betelgeuse, its smallest clumps would appear to be 5 to 10 larger than those we have detected.

There is a hint that the prevalence of clumps is higher close to the star than at the largest distances probed by our spectra. This encourages the speculation that the clumps’ origin may be traceable through the inhomogeneous inner shell (Lim et al. 1998) and chromosphere to the stellar surface. The greater visibility of shells at large distances may be due to the dissolution and merger of clumps at a common distance. Although clumps are seen by fluorescent emission, their shape and size are attributable to an enhanced concentration of potassium atoms in the wind and not to intensity variations over the surface of the illuminating photosphere; a bright photospheric spot will enhance fluorescence over approximately 2$`\pi `$ of solid angle not a mere few arcsec. The clumps are probably not regions of reduced ionization of potassium atoms because the characteristic timescale for photoionization is shorter that the expansion timescale (the time required to reach a given angular distance at a constant expansion velocity): it is a factor of 100 faster close to the star and about a factor of 10 at 50 arcsec from the star.

Supposing the clumps represent coherent structures related to the process of mass loss, a change of clump characteristics with radial distance provides information about the clumps’ internal structure. Consider a clump moving radially outwards at velocity $`v_{exp}`$ with an internal velocity dispersion $`v_{disp}`$. Its radial size $`\varphi `$ increases linearly with radial distance out from the star $`\mathrm{\Phi }`$, i.e., $`\mathrm{\Delta }\varphi `$/$`\mathrm{\Delta }\mathrm{\Phi }`$ = $`v_{disp}/v_{exp}`$. Identifying $`v_{disp}`$ with the sound speed, we estimate $`v_{disp}`$ = 1 km s-1. Since $`v_{exp}15`$ km s-1, a clump of about 2 arcsec at 10 arcsec from the star would expand to about 5 arcsec at a radial distance of 50 arcsec. In the inner shell, the inter-clump distances are a few arcsec. Then, if dispersion of the clump perpendicular to the radial flow occurs at about 1 km s-1, the clumps at a common distance may well merge after a few tens of arcsec.

Our observations are consistent with this result in that the maximum velocity extent of a velocity-position map is not greatly different for small and large impact parameters.

## 4 Concluding Remarks

Our study of the K i emission from Betelgeuse’s circumstellar shell reveals new aspects of the shell and its potassium atoms. Spherical symmetry is no more than a rough approximation to the distribution of potassium atoms but the radial distribution of the brightness closely follows the predictions for a model proposed by Rodgers & Glassgold (1991) implying a mass loss rate of about 2 $`\times 10^6`$ $`M_{}`$ y-1. But the striking result of our study is the presence of clumps of potassium atoms with the smallest clumps being unresolved spatially and kinematically. The clumps likely originate from cool structures close to the star. Observations suggest that they disappear at large distances from the star.

Observations at a site of excellent seeing are to be sought in order to trace the clumps back closer to the star. Other supergiants and giants are observable in both Na D and K i emission, as published reports demonstrate, but the occurrence of clumps in the winds from such stars is as yet unknown. To explore the innermost regions of these shells will require the spatial resolution of a space, or a diffraction limited ground-based telescope. Betelgeuse’s chromosphere has already been explored in the ultraviolet continuum and the Mg ii h and k lines (Uitenbroek, Dupree, & Gilliland 1998). Extension of this work to other lines, and other stars would be rewarding.

###### Acknowledgements.

We thank N. Mauron and E. Josselin for inspiring discussions of the Betelgeuse shell.


format: text

text:
# Non-Fermi-Liquid in a modified single electron transistor

## Abstract

At low temperatures, a system built from a small droplet of electrons and a larger, but still finite, droplet may display non-Fermi-liquid behavior. Stabilization of a multi-channel Kondo fixed point requires fine control of the electrochemical potential in each droplet. The desired fine control can be achieved by adjusting voltages on nearby gate electrodes. We study the conditions for obtaining this type of non-Fermi-liquid behavior and discuss the experimentally-observable consequences.

Despite the presence of strong electron-electron interactions, most “metallic” conductors are well-described by Landau‘s phenomenological Fermi-liquid (FL) theory . Among the few exceptions are various types of superconductors , Laughlin liquids in two dimensional electron systems at strong magnetic field , and Luttinger liquids in one dimension . In this Letter we argue that a simple configuration of two or more large electron droplets (see Fig. 1) attached to a small electron droplet can exhibit multi-channel Kondo (MCK) correlations , retaining non-Fermi-liquid (NFL) behavior at low temperature.

The single-channel Kondo (1CK) effect has been studied for decades in metals with magnetic impurities . The same phenomenon was observed recently in semiconductor nanostructures containing no magnetic impurities: instead, in each study an electron droplet with a degenerate ground state plays the role of a magnetic impurity, and nearby electron reservoirs play the role of the surrounding normal metal .

In contrast, due to an intrinsic channel anisotropy the MCK effect may not be observable in metals with magnetic impurities , making it an even more intriguing phenomenon to produce in artificial nanostructures . We will show below that in a certain configuration of electron droplets (Fig. 2), tuning the voltage on just one gate electrode can stabilize the two-channel Kondo (2CK) fixed point at low temperature. At this fixed point the conductance through the small droplet should have an anomalous power-law dependence on temperature, a manifestation of NFL behavior. The simplicity of the structure and the ability to tune system parameters hold out hope for detailed study of the NFL realm, including non-thermodynamic quantities such as transmission phase  and tunneling density of states.

To see how NFL behavior can be realized in a system of electron droplets, we first discuss a model in which a few large conducting droplets are attached to a single small droplet (see Fig. 1). A small central “pistil” droplet (denoted by $`d`$) has a single level of energy $`\epsilon _{ds}`$, which can be empty, or occupied by electrons of either or both spin directions $`s=,`$. Henceforth we refer to this droplet as small dot d. In the large “petal” droplets we neglect the discreteness of single-particle energy levels, while retaining a finite Coulomb energy. Thus, they behave as “interacting leads”; we refer to them as large dots.

To describe this system we use the model Hamiltonian:

$`H`$ $`=`$ $`{\displaystyle \underset{kls}{}}\epsilon _{lks}^{}l_{ks}^{}l_{ks}^{}+{\displaystyle \underset{s}{}}\epsilon _d^0d_s^{}d_s^{}+Un_d^{}n_d^{}`$ (1)
$`+`$ $`{\displaystyle \underset{l}{}}u_l^{}(n_l^{}𝒩_l^{})^2+{\displaystyle \underset{ks}{}}\left(V_{lk}^{}l_{ks}^{}d_s^{}+V_{lk}^{}d_s^{}l_{ks}^{}\right).`$ (2)
Here $`l_{ks}^{}`$ is the annihilation operator of an electron at state $`k`$, spin $`s`$ and energy $`\epsilon _{lks}^{}`$ at large dot $`l=1,\mathrm{},N`$, $`d_s^{}`$ is the annihilation operator of an electron with spin $`s`$ at small dot $`d`$, $`n_l^{}=_{ks}l_{ks}^{}l_{ks}^{}`$ is the number operator at dot $`l`$, and $`n_{d()}=d_{()}^{}d_{()}^{}`$. The parameter $`𝒩_l^{}`$ sets the equilibrium occupancy of dot $`l`$, and $`\epsilon _d^0`$ is the bare energy of level $`d`$. Each of these can be independently tuned by adjusting potentials on nearby gates.

Glazman and Raikh showed that a local site coupled to two noninteracting leads maps to a 1CK problem, in which only an even linear combination of the electron creation and annihilation operators in the two leads couples to the local spin . This result generalizes straightforwardly to arbitrary number of leads in the absence of Coulomb charging ($`u_l=0`$ for all $`l`$). The well-known 1CK system displays an interesting many-body resonance, however, at very low temperatures we can describe it simply as a FL superimposed with a resonance .

A nonzero $`u_l^{}`$ induces interactions among some of the transformed lead operators, substantially changing the low-temperature behavior of the system. In terms of the original large dots the physical role of the $`u_l^{}`$ terms is clear: at low temperatures they forbid processes in which charge is ultimately transferred from one large dot to another. Spin flip events — e.g. when an electron hops onto the small dot and then an electron with an opposite spin hops off the small dot to the same large dot — remain possible and may lead to a MCK fixed point.

To see this, we study the case $`\epsilon _d^0(\epsilon _d^0+U)<0`$, and integrate out in the renormalization group (RG) sense all energies up to a cutoff $`\stackrel{~}{D}a|\epsilon _d^{}||\epsilon _d^{}|,|\epsilon _d^{}+U|`$, with $`\epsilon _d^{}=\epsilon _d^0+\mathrm{\Gamma }/\pi \mathrm{log}(aD/|\epsilon _d^{}|)`$. Here $`\mathrm{\Gamma }`$ is the width of level $`d`$ in the absence of interactions, $`D`$ is a cutoff of order the Fermi energy, $`a`$ is a number of order 1, and we have assumed $`\epsilon _d^{}=\epsilon _d^{}\epsilon _d`$ . Now we can safely perform the Schrieffer-Wolff transformation  and find:

$`H`$ $`=`$ $`{\displaystyle \underset{kls}{}}\epsilon _{lks}^{}l_{ks}^{}l_{ks}^{}+{\displaystyle \underset{l}{}}u_l^{}(n_l^{}𝒩_l^{})^2`$ (3)
$`+`$ $`{\displaystyle \underset{lm,kq}{}}J_{lm}^{kq}\left[S_{}^+s_{lm}^{kq}+S_{}^{}s_{lm}^{+kq}+2S_{}^zs_{lm}^{zkq}\right],`$ (4)
where $`S_{}^\pm =d_{()}^{}d_{()}^{}`$, $`S_{}^z=\frac{1}{2}\left(d_{}^{}d_{}^{}d_{}^{}d_{}^{}\right)`$, $`s_{lm}^{\pm kq}=l_{k()}^{}m_{q()}^{}`$, $`s_{lm}^{zkq}=\frac{1}{2}\left(l_k^{}m_q^{}l_k^{}m_q^{}\right)`$ and

$$J_{lm}^{kq}=V_{lk}^{}V_{mq}^{}\left[\frac{1}{E_{\mathrm{elec}}^{qm}E_{\mathrm{init}}^{}}+\frac{1}{E_{\mathrm{hole}}^{kl}E_{\mathrm{init}}^{}}\right].$$
(5)
Here $`E_{\mathrm{init}}`$ is the energy of the initial state with one electron in the small dot $`d`$ and occupancy $`n_l(n_m)`$ in large dot $`l(m)`$; $`E_{\mathrm{elec}}^{qm}`$ denotes the energy of an intermediate state of a process where first an electron hops from state $`q`$ of dot $`m`$ onto dot $`d`$ and then the same electron or an electron with opposite spin hops off dot $`d`$ onto state $`k`$ of dot $`l`$; $`E_{\mathrm{hole}}^{kl}`$ is the energy of the intermediate state of a process where the temporal order of the hopping events is interchanged.

Let us assume that $`\epsilon _{l(m)qs}^{}\epsilon _{l(m)F}^{}`$, where $`\epsilon _{l(m)F}^{}`$ is the energy of the last empty (occupied) state in dot $`l(m)`$, and that $`V_{mq}=V_{lk}V`$. Then using Eq. (2) to express the energies in Eq. (5) we find:

$$J_{lm}^{kq}=J_{lm}=|V|^2\frac{U+u_m^{}+u_l^+}{\left[U+\epsilon _d+u_m^{}\right]\left[u_l^+\epsilon _d\right]},$$
(6)
where $`u_p^\pm =u_p\left[1\pm 2\left(n_p𝒩_p\right)\right]\pm \epsilon _{pF}^{}`$, and for $`1/2<𝒩_pn_p+(\mu \epsilon _{pF}^{})/(2u_p)<1/2`$ there are $`n_p`$ electrons in dot $`p`$, where $`p=m,l`$ and $`\mu `$ is the electrochemical potential of a reference reservoir. Particle hole symmetry may be absent in the large dots, so in general $`J_{lm}J_{ml}`$.

The off-diagonal terms describe transfer of an electron from dot $`m`$ to dot $`l`$; at low temperature these processes are exponentially suppressed as $`J_{lm}^{}=J_{lm}^0e^{u_{lm}^{}/(4T)}`$ with $`u_{lm}^{}=\left(u_l^++u_m^{}\right)\left(1\delta _{lm}\right)`$. Notice that $`u_{ll}^{}=0`$, since the charge on dot $`l`$ is not changed when an electron hops from dot $`l`$ onto dot $`d`$ and then back to the same dot $`l`$. At $`T<\mathrm{min}(u_{lm}^{})`$ only the diagonal terms of $`\stackrel{~}{J}_{lm}`$ do not flow to zero. Assuming that we are not at a degeneracy point where $`u_{lm}^{}=0`$, an easy condition to avoid, the RG equations are identical to the MCK RG equations :

$$\frac{d\stackrel{~}{J}_{ll}^{}}{d\mathrm{log}(\stackrel{~}{D}/T)}=\stackrel{~}{J}_{ll}^2\stackrel{~}{J}_{ll}^{}\underset{p}{}\stackrel{~}{J}_{pp}^2.$$
(7)
As in the case of classic MCK, our NFL fixed point is unstable to the introduction of channel anisotropy. If one of the coupling constants is larger than the others, the corresponding channel alone screens the local spin and forms a Kondo resonance while the other channels are decoupled from the local spin. In our model we can tune $`𝒩_l`$ to achieve $`J_{ll}=J`$ for all $`l`$. For $`N`$ large dots this requires tuning of $`N1`$ gate potentials.

These gate potentials capacitively control the energy of the last occupied level in each large dot, so excitations in each large dot will be around a different Fermi energy. This does not modify the RG equations, but will affect certain physical properties such as the small dot density of states at finite energies. A similar situation occurs in the discussion of 2CK in a dot out of equilibrium .

For large $`N`$ it would be difficult to realize MCK experimentally, requiring tuning at least $`N1`$ gate potentials to obtain the desired fixed point. We therefore propose a specific realization of the 2CK model (Fig. 2). In this “three leg” structure two noninteracting “free” leads $`L`$ and $`R`$ — effectively large dots with $`u_{L,R}=0`$ — are connected to small dot $`d`$ which is in turn connected to large dot $`m`$. Dot $`m`$ is not directly connected to the leads. This structure allows conductance measurements between leads $`L`$ and $`R`$ through dot $`d`$, and proper tuning of parameters should cause 2CK effect to proclaim itself in this conductance.

Let the potential $`𝒱_d`$ on the gate near dot $`d`$ be such that dot $`d`$ is singly-occupied in the absence of coupling to dot $`m`$. Then if $`\mathrm{\Gamma }_m`$, the bare coupling between dots $`d`$ and $`m`$, is roughly equal to $`\mathrm{\Gamma }_l`$, the coupling between dot $`d`$ and the leads, fine-tuning the gate potential $`𝒱_m`$ can stabilize the 2CK effect. To demonstrate this, we analyze the electrostatics of the circuit shown in the inset of Fig. 2.

For fixed number of electrons $`n_m(n_d)`$ on dots $`m(d)`$, the electrostatic energy of the system is:

$`E_{\mathrm{ele}}^{n_m,n_d}`$ $``$ $`E_{\mathrm{ele}}^{n_m,n_d}(𝒱_m,𝒱_d)`$ (8)
$`=`$ $`U\left(n_d𝒩_d\right)^2+u_m\left(n_d+\alpha n_m𝒩\right)^2,`$ (9)
where $`Ue^2/(2\stackrel{~}{C}_d)u_me^2/(2\stackrel{~}{C}_mc_\mu ^2/\stackrel{~}{C}_d)`$, $`|e|𝒩_dc_d𝒱_d`$, $`|e|𝒩c_m𝒱_m+c_dc_\mu /\stackrel{~}{C}_m𝒱_d`$ and $`\alpha c_dc_\mu /\stackrel{~}{C}_d`$. Here $`\stackrel{~}{C}_{m(d)}`$ is the total capacitance of dot $`m(d)`$. See Fig. 2 for definitions of the other capacitances.

Since dot $`m`$ is large we assume that $`\stackrel{~}{C}_m`$ is much larger than all other capacitances. The number of electrons on each dot can be estimated by minimizing $`E_{\mathrm{ele}}`$ with respect to $`n_m`$ and $`n_d`$ . In Fig. 3 we plot the number of electrons on each dot as a function of $`𝒱_m`$ and $`𝒱_d`$ .

To write down the full Hamiltonian $`\stackrel{~}{H}`$ of the model system we note that there is no Coulomb blockade for transfer of electrons between $`L`$ and $`R`$ leads, so we can define $`l_{ks}=\mathrm{cos}\theta L_{ks}+\mathrm{sin}\theta R_{ks},o_{ks}=\mathrm{cos}\theta R_{ks}\mathrm{sin}\theta L_{ks},\mathrm{tan}\theta =V_R/V_L,V_l=\sqrt{|V_L|^2+|V_R|^2}`$. Without loss of generality we take the coupling constants $`V_i`$, $`i=L,R,m`$, to be real. With these definitions:

$`\stackrel{~}{H}`$ $`=`$ $`{\displaystyle \underset{i=o,m,l;ks}{}}ϵ_{iks}^{}i_{ks}^{}i_{ks}^{}+{\displaystyle \underset{s}{}}ϵ_{ds}^{}d_s^{}d_s^{}+E_{\mathrm{ele}}^{n_m,n_d}(𝒱_m,𝒱_d)`$ (11)
$`+V_m{\displaystyle \underset{ks}{}}m_{ks}^{}d_s^{}+V_l{\displaystyle \underset{ks}{}}l_{ks}^{}d_s^{}+\mathrm{h}.\mathrm{c}.,`$
which is similar to Eq. (2). To obtain a 2CK fixed point we tune $`𝒱_m`$ and $`𝒱_d`$ to make $`\stackrel{~}{J}_{mm}`$, the coupling of dot $`d`$ to dot $`m`$, equal to $`\stackrel{~}{J}_{ll}`$, the coupling of dot $`d`$ to the leads:

$`\stackrel{~}{J}_{mm}(𝒱_m,𝒱_d)`$ $``$ $`\mathrm{\Gamma }_m\left[{\displaystyle \frac{1}{E_{\mathrm{ele}}^{1,0}E_{\mathrm{ele}}^{0,1}}}+{\displaystyle \frac{1}{E_{\mathrm{ele}}^{1,2}E_{\mathrm{ele}}^{0,1}}}\right]=`$ (12)
$`\stackrel{~}{J}_{ll}(𝒱_m,𝒱_d)`$ $``$ $`\mathrm{\Gamma }_l\left[{\displaystyle \frac{1}{E_{\mathrm{ele}}^{1,0}E_{\mathrm{ele}}^{0,2}}}+{\displaystyle \frac{1}{E_{\mathrm{ele}}^{0,0}E_{\mathrm{ele}}^{0,1}}}\right].`$ (13)
Here $`\mathrm{\Gamma }_{m(l)}=|V_{m(l)}|^2\nu _{m(l)}`$, where $`\nu _{m(l)}`$ is the density of states in dot $`m`$ (the leads). Eq. (13) together with the ratio $`\gamma \mathrm{\Gamma }_m/\mathrm{\Gamma }_l`$ defines a curve in the $`𝒱_m,𝒱_d`$ plane. In Fig. 3 we show several of these “2CK-lines” for different values of $`\gamma `$. On these lines 2CK physics should be realized at low $`T`$.

The 2CK fixed point can be reached experimentally by fixing $`𝒱_d`$ to give an odd number of electrons in the small dot, and fine tuning $`𝒱_m`$ so that $`\stackrel{~}{J}_{mm}(𝒱_m,𝒱_d)=\stackrel{~}{J}_{ll}(𝒱_m,𝒱_d)`$. In Fig. 3 this corresponds to tuning $`𝒱_m`$ until we hit the 2CK line for our given value of $`\gamma `$.

The current $`I`$ between the left and right leads can be measured as a function of temperature $`T`$, and as a function of $`𝒱_{LR}^{}`$, the bias applied between leads $`L`$ and $`R`$ (see Fig. 2). Drawing on the extensive literature of 2CK physics  we can predict the qualitative behavior of the $`I`$-$`𝒱_{LR}^{}`$ curve through the dot, for different values of the gate voltages $`𝒱_m`$, $`𝒱_d`$ that scan the the hexagon of Fig. 3.

The 2CK-lines divide the hexagon into three different parts. On the 2CK-lines in the unitary limit — $`T`$,$`𝒱_{RL}`$ Kondo temperature $`T_K\sqrt{U\mathrm{\Gamma }_m}e^{1/\stackrel{~}{J}_{mm}(𝒱_m,𝒱_d)}`$ — the differential conductance $`G(T,𝒱_{LR}^{})dI/d𝒱_{LR}^{}`$ should approaches its limiting value $`G(0,0)`$ as

$`2CK:G(0,0)G(T,𝒱_{LR}^{})\sqrt{\mathrm{max}(𝒱_{LR}^{},T)}.`$ (15)
In the symmetric case $`V_L=V_R`$, we get $`G(0,0)=G_Ke^2/(2\pi \mathrm{})`$, half the maximal value of $`G(0,0)`$ in the 1CK effect .

In regions of the hexagon adjoining boundaries labelled $`X`$ in Fig. 3, (i.e. for $`\stackrel{~}{J}_{mm}>\stackrel{~}{J}_{ll}`$), as temperature decreases the electrons in dot $`m`$ screen the spin of dot $`d`$, while the leads are decoupled. In the RG sense $`\stackrel{~}{J}_{ll}`$ flows to zero, so that dot $`m`$ “wins” over the leads and forms a 1CK state with dot $`d`$. In that case $`dI/d𝒱_{LR}^{}`$ is small and given by

$$\text{ Large dot wins:}G(T,𝒱_{LR}^{})\left[\mathrm{max}(𝒱_{LR}^{},T)\right]^2.$$
(16)
In contrast, in regions of the hexagon adjoining boundaries labelled $`Y`$ in Fig. 3 (i.e. for $`\stackrel{~}{J}_{mm}<\stackrel{~}{J}_{ll}`$), dot $`m`$ decouples from dot $`d`$ at low temperature, leaving the leads to form a 1CK resonance with dot $`d`$ and

$$\text{Leads win:}G(0,0)G(T,𝒱_{LR}^{})\left[\mathrm{max}(𝒱_{LR}^{},T)\right]^2,$$
(17)
where $`G(0,0)=2G_K`$ for $`V_L=V_R`$.

At sufficiently low temperature the finite level spacing $`\mathrm{\Delta }_m`$ in dot $`m`$ will cut off the RG flow of the coupling constants . We cannot make $`\mathrm{\Delta }_m`$ infinitesimal as we must retain a finite Coulomb blockade energy $`u_m>kT`$. However, the ratio between charging energy and level spacing can be made large, allowing 2CK behavior to be observed over a significant temperature range before the system finally flows to the 1CK fixed point.

In conclusion, if a small dot is coupled to multiple electron reservoirs (large dots) Coulomb blockade suppresses inter-reservoir charge transfer at low temperatures. In the simplest version of this system, electrostatic gates provide the tunability needed to stabilize a 2CK fixed point, resulting in observable NFL behavior. Softer versions of suppressing inter-reservoir tunneling could also work in place of Coulomb blockade. For example, the reservoirs could be conductors with large impedance , one-dimensional Luttinger liquids  or conductors with strongly-interacting charge carriers. Finally, while the channel asymmetry parameter is relevant in the RG sense, for realistically well-matched channel couplings we expect that the system will remain near the 2CK fixed point, and will show NFL behavior, over a wide range of temperatures.

It is our pleasure to acknowledge useful discussions with Natan Andrei, Michel Devoret, Gabi Kotliar, Leo Kouwenhoven, Yigal Meir, Arcadi Shehter, Alex Silva Ned Wingreen, and especially Avi Schiller. This work was supported by DIP grant c-7.1, by ISF grant 160/01-1, and by Stanford University.


text:
# Magnetic field induced charge and spin instabilities in cuprate superconductors

## Abstract

A $`d`$-wave superconductor, subject to strong phase fluctuations, is known to suffer an antiferromagnetic instability closely related to the chiral symmetry breaking in (2+1)-dimensional quantum electrodynamics (QED3). Based on this idea we formulate a “QED3 in a box” theory of local instabilities of a $`d`$-wave superconductor in the vicinity of a single pinned vortex undergoing quantum fluctuations. As a generic outcome we find an incommensurate 2D spin density wave forming in the neighborhood of a vortex with a concomitant “checkerboard” pattern in the local electronic density of states, in agreement with recent neutron scattering and tunneling spectroscopy measurements.

Among the open questions in modern condensed matter physics, few have inspired more theoretical effort than the emergence of a superconducting state from the doped antiferromagnetic (AF) insulator anderson1 . Recently, using an “inverted” approach to the problem balents1 ; ft1 , it has been shown that AF order arises naturally when the superconducting order in a $`d`$-wave superconductor ($`d`$SC) is destroyed by vortex-antivortex fluctuations herbut1 ; ft2 . As we shall discuss, the implications of these theories transcend the possibility of providing a route to understanding the destruction of superconductivity in strongly-underdoped cuprates; indeed, they also apply to the problem of local field-induced vortices within the superconducting state.

Recent neutron scattering lake1 ; lake2  and scanning tunneling spectroscopy (STS) davis1  experiments have revealed the presence of local AF and charge order in the vicinity of field-induced vortices. Existing theoretical treatments arovas1 ; demler1 ; ting1 ; zhang1 ; demler2  of vortex-induced AF ordering rely on the proximity of the system to a quantum critical point. Within such treatments, it is the suppression of the SC order parameter near the vortex cores that leads to the nucleation of islands of AF order. Here we present an alternative scenario in which the AF order is brought about by local quantum fluctuations of a vortex around its equilibrium position. In the present theory there is no competition between the SC amplitude and AF order: the latter arises purely from the presence of vortex fluctuations and is a genuine low-energy phenomenon taking place on lengthscales much longer than the core size.

It is a well-known fact that the low superfluid density in cuprates makes the SC order vulnerable to phase fluctuations emery1 ; corson1 ; ong1 . This observation has inspired theories in which the pseudogap state is modeled as a phase-disordered $`d`$-wave superconductor balents1 ; ft1 ; franz1 ; kwon1 , such that the demise of superconductivity is brought about by the unbinding and proliferation of the topological defects – vortices – in the phase of the SC order parameter. It has been pointed out ft1  that fluctuating vortices produce a non-trivial Berry-phase interaction between the quasiparticles of the underlying $`d`$SC. This interaction is described in terms of a massless non-compact gauge field $`a_\mu `$, minimally coupled to the Dirac fermions representing the low-energy quasiparticle excitations of the system. Within the theory of Ref. ft1  which maps the problem onto (2+1)-dimensional quantum electrodynamics (QED3), it is this interaction that destroys the Fermi liquid nature of quasiparticles in the pseudogap state and ultimately drives the AF instability herbut1 ; ft2 . Remarkably, both the ‘algebraic’ Fermi liquid describing the symmetric pseudogap phase and the antiferromagnet emerge from the same QED3 theory ft1 .

Here, we use the philosophy and formalism developed in Refs. ft1 ; herbut1 ; ft2  to model quasiparticle excitations in the superconducting state in the spatial region close to a single field-induced vortex undergoing fluctuations around its equilibrium position. We call this model “QED3 in a box”. We note that there exists direct experimental evidence that individual vortices indeed undergo significant quantum fluctuations hoogenboom1 . We find that, under generic conditions, interactions generated by such fluctuating vortex lead to local instability of the superconducting state which takes form of a 2D incommensurate AF spin density wave (SDW) with a wave vector tied to the positions of the nodes in the underlying $`d`$-wave gap.

In order to motivate our model for a single vortex we first review the treatment of the AF instability in QED3 and reformulate it in a way that will be more suitable for our present purposes. We start from the Euclidean QED3 action $`S=d^3x_D`$ with

$$_D\underset{l=1}{\overset{N}{}}\overline{\mathrm{\Psi }}_l(x)\gamma _\mu (i_\mu a_\mu )\mathrm{\Psi }_l(x)+_B[a(x)],$$
(1)
describing the low-energy fermionic excitations of a $`d`$-wave SC coupled to fluctuating vortices represented by the gauge field $`a_\mu `$ ft1 . Here, $`\mathrm{\Psi }_l(x)`$ is a four component Dirac spinor representing the fermionic excitations associated with a pair of antipodal nodes, $`x=(\tau ,𝐫)`$ denotes the space-time coordinate, and $`\gamma _\mu `$ ($`\mu =0,1,2`$) are the gamma matrices satisfying $`\{\gamma _\mu ,\gamma _\nu \}=2\delta _{\mu \nu }`$. The number $`N`$ of fermion species is equal to $`2`$ for single-layer cuprates; $`N=4,6,\mathrm{}`$ for bilayer, trilayer and multilayer materials. The Lagrangian $`_B`$ encodes the dynamics of the gauge field $`a_\mu `$ and is given by $`_B[a]=\mathrm{\Pi }_{\mu \nu }(q)a_\mu (q)a_\nu (q)`$ with

$$\mathrm{\Pi }_{\mu \nu }(q)=\left(m_a+\frac{N}{8}|q|\right)\left(\delta _{\mu \nu }\frac{q_\mu q_\nu }{q^2}\right).$$
(2)
The gauge field mass $`m_a`$ vanishes when vortices are unbound (i.e., in the pseudogap regime or, in the present situation, near a single fluctuating vortex) and is finite in the superconducting state where vortices appear only in tightly bound loops or pairs.

In the standard treatment herbut1 ; ft2  the AF order occurs via the phenomenon of chiral symmetry breaking pisarski1 ; appelquist1 ; nash1  in the QED3 Lagrangian (1). The instability is signaled by the spontaneous generation of fermion mass, $`m_D`$, which is interpreted in our context as the onset of SDW gap herbut1 ; ft2  for the original Bogoliubov quasiparticles. The most general, nonperturbative treatment of mass-generation in QED3 obtains $`m_D`$ as a solution of a self-consistent Dyson-Schwinger equation. Here we shall follow a slightly simpler route which leads to the same result and has the advantage of being more easily generalizable to the present problem. In Eq. (1) we integrate out the gauge field to obtain the following fermionic effective action:

$`S_{\mathrm{eff}}`$ $`=`$ $`{\displaystyle d^3x\overline{\mathrm{\Psi }}(x)\gamma _\mu i_\mu \mathrm{\Psi }(x)}`$
$``$ $`{\displaystyle d^3xd^3yJ_\mu (x)D_{\mu \nu }(xy)J_\nu (y)},`$
where $`J_\mu (x)\overline{\mathrm{\Psi }}(x)\gamma _\mu \mathrm{\Psi }(x)`$ is the fermion 3-current and $`D_{\mu \nu }(x)`$ is the Fourier transform of the gauge boson propagator $`D_{\mu \nu }(q)=\mathrm{\Pi }_{\mu \nu }^1(q)`$. Henceforth we shall focus on a single pair of nodes and thus drop the nodal index $`l`$. The integrand of the interaction term may be written as $`D_{\mu \nu }(xy)\mathrm{Tr}[\mathrm{\Psi }(y)\overline{\mathrm{\Psi }}(x)\gamma _\mu \mathrm{\Psi }(x)\overline{\mathrm{\Psi }}(y)\gamma _\nu ]`$ where the trace is taken over the spinor indices. This form suggests a Hartree-Fock (HF) approach in which we decouple the 4-fermion interaction to obtain $`D_{\mu \nu }(xy)\mathrm{Tr}[\mathrm{\Psi }(y)\overline{\mathrm{\Psi }}(x)\gamma _\mu G_0(x,y)\gamma _\nu ]`$ with $`G_0(x,y)=\mathrm{\Psi }(x)\overline{\mathrm{\Psi }}(y)`$ and the average is taken with respect to the HF effective action to be specified shortly. To make the structure of the interaction term more transparent we utilize the relative and center of mass coordinates $`r=xy`$ and $`R=(x+y)/2`$ to write it as

$$d^3Rd^3r\mathrm{Tr}\left[\mathrm{\Psi }(R_+)\overline{\mathrm{\Psi }}(R_{})\gamma _\mu G_0(R,r)\gamma _\nu \right]D_{\mu \nu }(r),$$
(4)
where $`R_\pm R\pm \frac{r}{2}`$. In the uniform system the Green’s function is independent of $`R`$, $`G_0(R,r)=G_0(r)`$. Furthermore, both $`G_0(r)`$ and $`D_{\mu \nu }(r)`$ are strongly peaked at $`r0`$. The dominant contribution to the interaction therefore comes from this region and we may write (4) as remark1

$$d^3R\overline{\mathrm{\Psi }}(R)\mathrm{\Psi }(R)d^3r\gamma _\mu G_0(r)\gamma _\nu D_{\mu \nu }(r).$$
(5)
We have dropped the trace since the interaction is proportional to the unit matrix in the spinor space.

Inspection of Eq. (5) suggests the following HF effective action and self-consistency condition:

$`S_{\mathrm{HF}}`$ $`=`$ $`{\displaystyle d^3x\overline{\mathrm{\Psi }}(x)(\gamma _\mu i_\mu im_D)\mathrm{\Psi }(x)},`$ (6a)
$`im_D`$ $`=`$ $`{\displaystyle \frac{1}{4}}\mathrm{Tr}{\displaystyle d^3r\gamma _\mu G_0(r)\gamma _\nu D_{\mu \nu }(r)}.`$ (6b)
The last integral is easily evaluated by going to momentum space and Eq. (6b) becomes $`m_D=(8m_D/N\pi ^2)\mathrm{ln}(\mathrm{\Lambda }/m_D)`$, where $`\mathrm{\Lambda }`$ is the high-momentum cutoff. This yields a nontrivial solution

$$m_D=\mathrm{\Lambda }e^{N\pi ^2/8},$$
(7)
in agreement with the classic result of Pisarski pisarski1 . More sophisticated treatments appelquist1 ; nash1  based on the Schwinger-Dyson equation give a finite critical value of $`N_c`$ above which no mass is generated; however, for our purposes the level of approximation embodied by Eq. (7) will be sufficient.

We have thus seen that, in a uniform system, fluctuating vortices lead to the formation of SDW order. The challenge we now face is twofold: (i) we must adapt the above treatment to the case of a single fluctuating vortex, and (ii) since we seek to study the commensuration effects present in real materials, we must formulate the corresponding theory on the lattice. To address (i) let us denote by $`\mathrm{}_v`$ the characteristic length scale over which the vortex fluctuates around its classical equilibrium position. Within this length scale, the Berry-phase interaction between quasiparticles (and hence tendency towards SDW ordering) will be strong. We model this by taking in this region the gauge field to be massless. On the other hand at distances well beyond $`\mathrm{}_v`$, quasiparticles feel no interaction and we model this by gauge field having a large mass $`m_a`$. In particular we take,

$$m_a(𝐑)=\mathrm{\Delta }_0\left(\frac{|𝐑|}{\mathrm{}_v}\right)^n,$$
(8)
where $`\mathrm{\Delta }_0`$ is an energy scale which we take to be the maximum superconducting gap, $`|𝐑|`$ is the distance from the vortex equilibrium position and $`n`$ is a positive exponent. (We use $`n=2`$ but our numerical calculations below are largely insensitive to the exact value of $`n`$.)

To address (ii), (i.e., to put the theory on the lattice) we recall that the effective action (1) and its HF version Eq. (6a) descend from a model of a lattice $`d`$SC linearized near the nodes of the gap. We therefore consider the corresponding lattice Hamiltonian enriched by the “mass” term present in Eq. (6a), to represent the HF decoupled Berry phase interaction. Thus we have,

$$H_{\mathrm{HF}}=\underset{\sigma }{}\underset{ij}{}\mathrm{\Phi }_{i\sigma }^{}_{ij}^\sigma \mathrm{\Phi }_{j\sigma }.$$
(9)
Here $`\mathrm{\Phi }_{i\sigma }^{}=(c_{i\sigma }^{},c_{i\overline{\sigma }})`$, $`c_{i\sigma }^{}`$ represents the electron creation operator at lattice site $`i`$, spin index $`\overline{\sigma }=\sigma `$, and

$$_{ij}^\sigma =\left(\begin{array}{cc}t_{ij}+\delta _{ij}(m_{i\sigma }\mu )& \mathrm{\Delta }_{ij}\\ \mathrm{\Delta }_{ij}^{}& t_{ij}\delta _{ij}(m_{i\overline{\sigma }}\mu )\end{array}\right),$$

with $`t_{ij}`$ the tight binding hopping amplitude, $`\mathrm{\Delta }_{ij}`$ the SC gap, $`\mu `$ chemical potential, and $`m_{i\sigma }`$ the local spin magnetization representing the mass gap $`m_D`$ in Eq. (6a).

We diagonalize $`H_{\mathrm{HF}}`$ by the generalized Bogoliubov transformation $`c_{i\sigma }=_n[u_{n\sigma }(𝐫_i)\gamma _{n\sigma }+\sigma v_{n\overline{\sigma }}^{}(𝐫_i)\gamma _{n\overline{\sigma }}^{}]`$, where $`\chi _{n\sigma }(𝐫_i)[u_{n\sigma }(𝐫_i),\sigma v_{n\sigma }(𝐫_i)]^T`$ satisfy

$$\underset{j}{}_{ij}^\sigma \chi _{n\sigma }(𝐫_j)=ϵ_{n\sigma }\chi _{n\sigma }(𝐫_i).$$
(10)
In terms of the $`\chi _{n\sigma }`$, the self-consistency condition (6b) can be written as

$`m_{i\sigma }`$ $`=`$ $`{\displaystyle \underset{n\sigma j}{}}\sigma f(ϵ_{n\sigma })V_i(𝐫_j)u_{n\sigma }^{}(𝐑_i+𝐫_j)u_{n\sigma }(𝐑_i𝐫_j),`$
$`V_i(𝐫)`$ $``$ $`{\displaystyle \frac{1}{4}}{\displaystyle _0^{\mathrm{}}}𝑑\tau e^{\tau ϵ_{n\sigma }}\mathrm{Tr}\left[\gamma _\mu D_{\mu \nu }(\tau ,𝐫)\gamma _\nu \right].`$ (11)
For $`m_a0`$ the above integral cannot be evaluated in closed form. However, we find that it can be accurately approximated by a simple interpolation formula

$$V_i(𝐫)𝒱_0\frac{c_1}{r(rm_a(𝐑_i)+c_1)(rϵ_{n\sigma }+c_1)}$$
(12)
where $`r=|𝐫|`$, $`c_1=2/\pi `$ and $`𝒱_0=16/N\pi ^2`$ for the case of an isotropic Dirac cone ($`t=\mathrm{\Delta }`$). In the physical case $`t>\mathrm{\Delta }`$ the constant $`𝒱_0`$ will be modified somewhat and in what follows we treat it as an adjustable parameter of the model measuring the strength of the interaction. It is interesting to note that, as seen from Eq. (12), in (2+1)D a gauge field mass does not lead to exponentially decaying interactions on long length scales.

To capture the effect of vortex fluctuations on the local superconducting order we solve the eigenproblem (10) numerically on a lattice of $`M\times M`$ sites and iterate to self-consistency using Eq. (11). For simplicity we consider only nearest-neighbor hopping amplitudes $`t_{ij}=t`$ and a uniform $`d`$-wave gap $`\mathrm{\Delta }_{ij}=\pm \mathrm{\Delta }_0`$, with $`+`$ and $``$ signs referring to vertical and horizontal bonds respectively. We emphasize that the vortex, fluctuating around the equilibrium position at the center of our lattice, enters through the position dependent gauge-field mass $`m_a(𝐑)`$ given by Eq. (8), which in turn enters the potential $`V_i(𝐫)`$ given by Eq. (12). In the spirit of our working philosophy that the SDW order arises from the vortex fluctuations (and therefore from the gauge field), we neglect at this stage any effects of the superflow around the vortex or suppression of $`\mathrm{\Delta }_{ij}`$ in the core. Such effects are well understood and will be included in a future publication. We also neglect the effects of changes in the fermionic spectrum due to the onset of SDW on the interaction mediated by Berry gauge field Eq. (12).

The diagonalizations are performed using standard LAPACK routines, which allow us to handle systems up to $`40\times 40`$ sites. Typically, 10-15 iterations are needed to ensure self-consistency in $`m_{i\sigma }`$. We use both periodic and free boundary conditions and find that they have negligible effect on the results reported below. Our typical results are summarized in Figs. 1 and 2, showing the spatial distributions of the spin magnetization $`M_i=_\sigma \sigma c_{i\sigma }^{}c_{i\sigma }`$, staggered spin magnetization $`M_i^S=M_i(1)^{x_i+y_i}`$, local electron charge density $`n_i=_\sigma c_{i\sigma }^{}c_{i\sigma }`$, and energy integrated local density of states (LDOS) $`S_{E_1}^{E_2}(i)=_{E_1}^{E_2}\rho _i(E)𝑑E`$ where $`\rho _i(E)`$ is the LDOS at site $`i`$, as well as their respective Fourier transforms (FTs).

Panel (a) in Fig. 1 illustrates the “2D” incommensurate SDW pattern emerging in the vicinity of a fluctuating vortex with an $`8\times 8`$ unit cell containing islands of AF order separated by anti-phase domain walls \[apparent in panel (c)\]. The FT displayed in panel (b) reveals that this pattern can be thought of as a superposition of four 1D SDWs with wave vectors $`𝐐_{\mathrm{SDW}}=\pi (1\pm \delta _{\mathrm{SDW}},1\pm \delta _{\mathrm{SDW}})`$, $`\delta _{\mathrm{SDW}}=\frac{1}{4}`$. The size of $`\delta _{\mathrm{SDW}}`$ is doping dependent: it shrinks with increasing $`\mu `$ and vanishes at half filling $`(\mu =1)`$, giving rise to perfectly commensurate AF SDW. We also find that for $`\mu <1.7`$ the SDW becomes very weak for reasonable values of coupling $`𝒱_0`$: overdoped samples are less susceptible to AF instability.

According to general symmetry arguments, a spatial modulation in the spin density generates a modulation in the charge density, $`\delta n_iM_i^2`$. For our 2D SDW pattern this implies that the corresponding CDW will have a unit cell with half the area, rotated by $`45^{}`$ relative to the unit cell of $`M_i^S`$. Indeed, panels (a) and (b) of Fig. 2 confirm this general expectation, showing a “checkerboard” CDW at principal wavevectors $`𝐐_{\mathrm{CDW}}=\pi (\pm \delta _{\mathrm{CDW}},0),\pi (0,\pm \delta _{\mathrm{CDW}})`$ with $`\delta _{\mathrm{CDW}}=\frac{1}{2}`$. A similar checkerboard pattern arises in the integrated LDOS and is displayed in panels (c) and (d).

Our findings of a checkerboard pattern in LDOS are consistent with the recent STS experiments performed on Bi2Sr2CaCu2O8+δ (BSCCO) crystals davis1 . Our prediction is that the period of the pattern should increase with underdoping and that the effect should vanish in the overdoped samples. Also, if the observed LDOS pattern is associated with electron density modulation in a single CuO layer, we predict that the corresponding neutron scattering peaks should be found at wavectors $`𝐐_{\mathrm{SDW}}=\pi (1\pm \frac{1}{4},1\pm \frac{1}{4})`$. We note that neutron experiments lake1 ; lake2  on La1.84Sr0.16CuO4 (LSCO) show peaks at different $`k`$-space positions, $`\pi (1\pm \frac{1}{4},1)`$ and $`\pi (1,1\pm \frac{1}{4})`$. Although the findings of STS and neutron experiments are generally cited as being mutually consistent our analysis above indicates that this is not necessarily so: for a genuine 2D SDW illustrated in Fig. 1 determination of the corresponding CDW must consider the interference terms which cause the apparent $`45^{}`$ rotation of the latter. In the absence of neutron measurements on BSCCO we see two possible resolutions of this difficulty. First, it may be that the CDW pattern is truly 2D and neutron scattering on BSCCO would find a pattern illustrated in Fig. 1(b). Second, it could be that STS sees an incoherent superposition of two orthogonal 1D CDWs originating in two CuO layers comprising the BSCCO bilayer. This would explain the $`xy`$ anisotropy reported in Ref. davis1  and the neutron pattern would be consistent with that observed in LSCO. Within our simple model such 1D solutions have slightly higher energy than the 2D solutions reported above but it is possible that in a more complete model (using e.g. a more realistic band structure) the situation will be reversed.

To summarize, we have presented a “QED3 in a box” theory for field induced spin and charge instabilities in cuprates, driven by local fluctuations of a pinned vortex. Without any need to fine-tune parameters we find LDOS patterns in detailed agreement with the tunneling data on BSCCO davis1  and we relate them in a plausible way to existing neutron experiments lake1 ; lake2 . Other models arovas1 ; demler1 ; ting1 ; zhang1 ; demler2  rely on the local suppression of the superconducting order parameter and would therefore predict similar instabilities in the vicinity of impurities, grain boundaries, and sample edges where thus far no such effects have been observed.

This research was supported in part by NSERC (MF,DES) and NSF Grant DMR00-94981 (ZT).