text:
# Gravitational Microlensing of Gamma Ray Bursts at Medium Optical Depth

## 1 Introduction

Press & Gunn (1973) were the first to propose that a cosmological abundance of dark compact objects could be detected by gravitational lensing of more distant sources. If a point mass lies along the observer-source line of sight, the relative motion between the lens, source and observer produces a change in the magnifications of two lensed images. The presence of foreground compact objects is therefore detected through a change in the observed flux of the background source (termed microlensing). This effect has been used successfully in the search for compact objects in the halo of the Milky Way galaxy (e.g. Alcock et al. 2000). Also, microlensing due to stars in a galaxy at moderate redshift has been observed in the gravitationally lensed quasar Q2237+0305 (Irwin et al. 1989; Corrigan et al. 1991). The short durations of gamma ray bursts (GRBs) offer an alternative way to study microlensing and hence search for a cosmological population of compact objects, through observation of repeating bursts.

In this paper we assume GRBs to be at cosmological distances (Paczynski 1995), and therefore probable sources for gravitational lensing (Turner, Ostriker & Gott 1984). Paczynski (1986, 1987) noted that while multiple images of a single GRB cannot be angularly resolved by present day detectors, their relative delay may be longer than the burst duration so that the lensed images could be resolved temporally (a pair of lensed GRB images that are produced by a single mass $`(M_{CO})`$ have a relative time delay that is $`\mathrm{\Delta }t50sec\times (M_{CO}/10^6M_{})`$ (Mao 1992)). Microlensing of a GRB by a compact object is therefore observed as a GRB that repeats. The utility of GRBs to explore cosmological dark matter, as well as the possible inference of properties of the GRB population itself was discussed in detail by Blaes & Webster (1992). The short event duration, as well as the transparency of the universe to gamma rays make GRBs ideal probes of dark matter in the form of compact objects over a wide range of masses.

Microlensing of existing and potential catalogues of GRBs have been used to discuss the cosmological abundance of compact objects. Marani et al. (1999) use non-detections of lensed images from the BATSE and Ulysses catalogues to set conservative limits on dark compact objects with masses between $`10^{16}`$ and $`10^7M_{}`$. Also a universe proposed by Gnedin & Ostriker (1992) with $`M_{CO}10^{6.5}M_{}`$ and $`\mathrm{\Omega }_M=\mathrm{\Omega }_{CO}=0.15`$ was ruled out at a confidence level of 90%. This scenario had been previously investigated in detail, using a 3-D lensing code by Mao (1993). He found results that did not depart significantly from those obtained by a single-screen approximation. Mao (1992) and Grossman & Nowak (1994) estimate that the waiting time for one lensed pair due to an intervening galaxy to be observed in the BATSE catalogue is between one and 10 years, and conclude that it is not certain that such a lensed pair will be found by BATSE.

As in quasar lensing, the detection of a lensed pair is separated from spatially coincident events by the comparison of light-curves and spectra. However the presence of noise and the faintness of images will alter the light curves so that macro images appear dissimilar (Wambsganss 1993; Nowak & Grossman 1994). In addition, if a foreground galaxy is responsible for the lensed pair then Williams & Wijers (1997) find that microlensing by individual stars in that galaxy can smear out image light curves, further increasing the chance of their being mis-classified as two spatially coincident events. Light-curve similarity of a lensed pair requires the source to be much smaller than the Einstein ring radius of the lens (projected into the source plane), a condition likely to be filled in the case of GRBs due to the short event duration for all but the smallest lenses. However an anisotropic source (eg. a beamed source) is effectively viewed at two different angles in the two images. This raises the possibility that the time variation of the two images may be different (eg. Babul, Paczynski & Spergel 1987), although Blaes & Webster (1992) find that a source should be isotropic to lensing provided that the beaming is not too strong. Paczynski (1987) noted that the potential difficulties in identifying lensed pairs will make confirmation of the cosmological origin of the GRBs difficult through a lensing argument. On the other hand, if we assume that GRBs have a cosmological origin, the mis-classification of events will lead to an underestimate of a compact object population.

The microlensing effect of a cluster of objects on a GRB was first discussed by Paczynski (1987). He found that a single instantaneous burst lensed by a screen of objects such as that produced by a galactic halo or a cluster of galaxies results in the observation of many repeats of the original burst separated by different delays. Examples of the light-curves produced can be found in Paczynski (1987) and Williams & Wijers (1997). At moderate optical depths the first images have magnifications that are not correlated with their delay, however the faint images arriving later have observed fluxes that decrease monotonically with time. Williams & Wijers (1997) investigated the flux weighted $`rms`$ of the delay of clusters of microimages and found this to be a sensitive function of optical depth and shear. The temporal spread of the cluster of microimages was shown to be related to the longest side of the area of the image plane containing microimages that account for most of the macroimage flux. This scale length is considerably larger than that of the separation of a lensed pair produced by a single microlens, and so the spread in the flux weighted arrival time is also considerably larger than the relative delay between images produced by a single mass. This phenomenon allows smaller masses to be probed modulo a detector resolution and intrinsic GRB duration.

A non-negligible fraction of the baryon content of the universe may exist in stellar remnant form (eg. Kerins & Carr 1994). In addition, Galactic microlensing searches (e.g. Alcock et al. 2000) provide evidence for a significant fraction of the Galactic halo mass being in the form of stellar mass compact objects, an interpretation which is supported by observations of high proper motion, cold white dwarfs (Ibata et al. 2000; Hodgkin et al. 2000). While intervening galaxies produce relatively large delays, GRBs with micro-second durations or variability would resolve microlensing by cosmological stellar mass objects (Nemiroff, Norris, Bonnell & Marani 1998). Nemiroff et al. (1998) note that there is no known fundamental reason for the non-existence of micro-second GRBs, and find that a significant rate of these might exist and go undetected by current telescopes. They find that at a given flux level there may be more GRBs with durations between 1 and 2 milliseconds than between 8.192 and 16.384 seconds (BATSE duration bin), and an abundance of microsecond spikes that is an order of magnitude greater. However, Nemiroff et al. (1998) also note that the converse possibility of short duration spikes not existing at all is equally consistent with the analysis of current data.

We propose that if short duration bursts were to exist then they could be used to probe the cosmological stellar mass compact object population without resolving and identifying individual micro-images. Also, in the absence of short spikes, the light-curves of GRBs can place upper limits on typical compact object mass and optical depth. Rather than considering the probability of lensing by a single mass on individual GRBs, we consider the microlensing effect on an ensemble of GRB light curves of a collection of point masses at low to moderate optical depth. For a sample of microlensed spike GRBs we look at the average value of the cumulative flux as a function of time.

This paper presents calculations corresponding to the microlensing effect of populations of stellar mass objects on the observed light-curves of hypothetical, very-short duration GRBs. However, the results are also applicable to the presently known population of GRBs, with durations of 10s of seconds as probes of compact objects having masses of $`10^6M_{}`$. The relevant scaling is pointed out where appropriate.

The remainder of the paper is presented in 3 parts. Section 2 describes the microlensing models and section 3 proposes a new method for analysing microlensing in GRBs. A brief summary is presented in section 4.

## 2 Microlensing models

As a crude approximation to the lensing effect of compact objects distributed along the line of sight, we assume a screen of point masses distributed randomly in a disc. The model does not include a continuous matter distribution. We use standard notation for gravitational lensing. The Einstein radius of a microlens in the image plane is denoted by $`\xi _o`$ and when projected into the source plane by $`\eta _0`$. The normalised convergence or optical depth is denoted by $`\kappa `$. The normalised lens equation for a field of point masses is

$$\stackrel{}{y}=\left(\begin{array}{cc}1\gamma & 0\\ 0& 1+\gamma \end{array}\right)\stackrel{}{x}+\underset{j=0}{\overset{N_{}}{}}m^j\frac{(\stackrel{}{x}^j\stackrel{}{x})}{|\stackrel{}{x}^j\stackrel{}{x}|^2}$$
(1)
Here $`\stackrel{}{x}`$ and $`\stackrel{}{y}`$ are the normalised image and source positions respectively, $`\gamma `$ is the applied shear, and the $`\stackrel{}{x_i}^j`$ and $`m^j`$ are the normalised positions and masses of the individual microlenses. To construct a microlensed light-curve for an instantaneous microlensed GRB, Eqn. 1 is solved for the micro-image positions at many points along a predefined source line through the inversion technique of Lewis et al. (1993) and Witt (1993). The time delay and magnification are then determined for each micro-image $`i`$, given by

$`\mathrm{\Delta }T_i`$ $`=`$ (2)
$`{\displaystyle \frac{\xi _o^2}{c}}{\displaystyle \frac{D_s}{D_dD_{ds}}}(1+z_d)\left({\displaystyle \frac{(\stackrel{}{x}_i\stackrel{}{y})^2}{2}}{\displaystyle \underset{j=0}{\overset{N_{}}{}}}m^jln(|\stackrel{}{x}_i\stackrel{}{x}^j|)\right)`$
and

$`\mu _i`$ $`=`$ $`{\displaystyle \frac{1}{|detA(\stackrel{}{x}_i)|}}`$
$`\mathrm{w}heredetA(\stackrel{}{x}_i)={\displaystyle \frac{y_1}{x_1}}|_{\stackrel{}{x}_i}{\displaystyle \frac{y_2}{x_2}}|_{\stackrel{}{x}_i}\left({\displaystyle \frac{y_2}{x_1}}|_{\stackrel{}{x}_i}\right)^2.`$
Here $`D_s,D_d`$ and $`D_{ds}`$ are the angular diameter distance between the observer and the source, the observer and the lens, and the lens and source. $`c`$ is the speed of light and $`z_d`$ the lens redshift. We note that the delay is proportional to ($`1+z_d`$), but is not explicitly dependent on the source redshift. Also, due to the dependence on $`\xi _o^2`$, $`\mathrm{\Delta }T`$ is proportional to the mean mass.

At a sufficiently large angle from the point source, there is a low magnification image located very close to each point mass. All solutions of the lens equation must therefore be found in a region that contains a sufficient percentage of the total macroimage flux. The region of the lens plane in which image solutions need to be found to ensure that $`99.9\%`$ of the total macro-image flux is recovered from all points on the source line is known as the shooting region. The number of stars in the region about any point which collects 99.9% of the macro-image flux was calculated by Katz, Balbus & Paczynski (1986) and is given by:

$$N_{}=300\frac{m^2}{m^2}\frac{\kappa ^2}{|(1\kappa )^2\gamma ^2|}.$$
(4)
In the absence of shear ($`\gamma =0)`$ these stars are distributed in a disc with a radius $`R_{sd}=\sqrt{\kappa /(N_{}\times m)}`$. The shooting region is defined by the union of these discs centred on the points $`(x_1=y_1/(1\kappa ),x_2=y_2/(1\kappa ))`$ corresponding to all parts of the source line. The minimum number of stars ($`N_{min}`$) required for the model are contained in a disc having a radius $`R_{}`$ which covers the shooting region. We have used 500 stars in each of our models which is larger than $`N_{min}`$ in all cases. We assume that the source is stationary with respect to the microlenses for the duration of the microlensed light-curve. We find the light curves for 10 bursts per microlens field distributed along a source track of length 10$`\eta _o`$. At optical depths between $`\kappa =0.025`$ and 0.25 (at intervals of 0.025) we compute the light curves for $`10^4`$ bursts.

The most probable scenario for lensing has the lens lying at a distance that is a reasonable fraction of that of the source (Turner, Ostriker & Gott 1984). In addition, the distribution of optical depth with redshift is reasonably sharply peaked. With this in mind we stress that our calculation uses a single screen to approximate the lensing effect of all masses along the line of sight. As an example we place the population of model GRBs at a redshift of $`z_s=1`$ and a screen of compact objects at a redshift of $`z_d=\frac{2}{3}`$. For simplicity we assume that all compact objects have a unique mass $`m`$.

## 3 Method of analysis

In this section we consider the analysis of a hypothetical sample of spike GRBs that have an intrinsic duration smaller than 0.1 milliseconds which we take as the detector binsize. From our sample of model microlensed light curves we find the cumulative flux (summed over all located microimages) that has arrived by the end of each bin. One minus the cumulative flux divided by the total flux gives the fraction of total flux that is yet too arrive $`(F_i(\mathrm{\Delta }T))`$ for each burst $`i`$ as a function of time. These curves are averaged over many bursts:

$$F(\mathrm{\Delta }T)\underset{i=0}{\overset{N_{bursts}}{}}F_i(\mathrm{\Delta }T)/N_{bursts}$$
(5)
Eqn. 5 describes the average behaviour of the arrival time of lensed flux for a set of microlensing parameters $`\kappa `$ and $`m`$. Note the assumption that all lines-of-sight have the same average optical depth $`\kappa `$. From the work on clusters of GRB microimages by Paczynski (1987) and Williams & Wijers (1997) we expect $`F(\mathrm{\Delta }T)`$ to have the following characteristics. Since the microlensing spread of the GRB is related to the size of the shooting region, which is proportional to the Einstein radius of the microlenses and therefore to $`\sqrt{m}`$, $`\mathrm{\Delta }T(m)`$ such that $`F(\mathrm{\Delta }T)=const`$ is proportional to $`m`$. Secondly, at larger optical depths there are more microimages with large magnification. A given fraction of microlensed flux therefore arrives later at higher optical depths and so $`F(\mathrm{\Delta }T)`$ is larger for increased $`\kappa `$. Note that when in a region where $`(1\kappa )^2\gamma ^20`$, the magnification becomes very large, and flux from the burst will be observed at a comparable level for a time similar to the geometric delay for a trajectory at the edge of the region.

Fig. 1 shows the fraction of flux that is yet to arrive as a function of time. Cases are shown for mean masses of 1$`M_{}`$ and 10$`M_{}`$, at optical depths of $`\kappa =`$0.025 (thin lines) and $`\kappa =`$0.250 (thick lines). At each optical depth, the larger mass produces longer microimage delays, causing a given fraction of flux to arrive later on average. In addition, at a fixed mass a given fraction of flux arrives later at larger optical depths producing a curve that has a smaller initial drop in the first bin followed by a decline that is more rapid than that at smaller optical depth. However at higher $`\kappa `$, $`F(\mathrm{\Delta }T)`$ is larger for all $`\mathrm{\Delta }T`$.

We define $`T_{0.99}`$ as the time $`\mathrm{\Delta }T`$ at which 99% of the total microlensed flux from the GRB has arrived. For binsizes that are small with respect to the microlensed spread of the event, $`T_{0.99}`$ scales linearly with microlens mass. $`T_{0.99}`$ is also a function of $`\kappa `$ and provides a natural scaling unit for time. The left hand panel in Fig. 2 re-displays the curves from Fig. 1 with time normalised by $`T_{0.99}`$ in each case ($`F(\mathrm{\Delta }T_{norm})`$ where $`\mathrm{\Delta }T_{norm}=\mathrm{\Delta }T/T_{0.99}`$). At the end of each bin the normalised curves are independent of mass. The curvature therefore provides a probe of the optical depth. The right hand panel of Fig. 2 displays $`F(\mathrm{\Delta }T_{norm})`$ for optical depths between $`\kappa =`$0.025 and 0.250 at a microlens mass of $`10M_{}`$. We suggest the value of $`F(\mathrm{\Delta }T_{norm}=0.1)`$ as an indicator of optical depth.

For values of $`\kappa `$ and $`m`$ in the ranges $`0.025\kappa 0.25`$ and $`0.1m100`$, we made 1000 simulations of $`F(\mathrm{\Delta }T)`$ (each $`F(\mathrm{\Delta }T)`$ was calculated using 100 model microlensed GRB light-curves). Fig. 3 shows average values of the quantities $`T_{0.99}`$ (dark contours) and $`F(\mathrm{\Delta }T_{norm}=0.1)`$ (light contours) over a range of optical depths and microlens masses. In regions combining small values of optical depth with microlens masses $`m<1M_{}`$, the contours are nearly parallel, however for most of the parameter space, values of $`T_{0.99}`$ and $`F(\mathrm{\Delta }T_{norm}=0.1)`$ describe a unique set of $`\kappa `$ and $`m`$. For samples containing 100 GRBs the variance in values of $`T_{0.99}`$ and $`F(\mathrm{\Delta }T_{norm}=0.1)`$ are an order of magnitude lower than the corresponding mean.

We have demonstrated that a sample of GRBs could be used to measure the quantities $`T_{0.99}`$ and $`F(\mathrm{\Delta }T_{norm}=0.1)`$, and that this combination corresponds to measurements of $`\kappa `$ and $`m`$. However it will not be known whether the observed spread in the arrival time of GRB flux is due solely to microlensing effects or whether it is intrinsic to the source (unless the individual micro-images are resolved). $`T_{0.99}`$ and $`F(\mathrm{\Delta }T_{norm}=0.1)`$ will therefore be upper bounds on the values for instantaneous bursts, and so will exclude a region of mass - optical depth parameter space rather than measure probable values. In addition, we note that Fig. 3 presents results covering only 3 orders of magnitude in mass. However the binsize, mass and $`T_{0.99}`$ are all linearly related, allowing scaling of the mass and $`T_{0.99}`$ by the binsize divided by 0.1 milliseconds. Estimates of upper limits for $`T_{0.99}`$ and $`F(\mathrm{\Delta }T_{norm}=0.1)`$ which will probe $`10^6M_{}`$ compact objects are therefore possible using current data with GRB durations of $`10`$ seconds.

The approach described requires that the source be integrated for long enough to be sure that the entire event flux (or some appropriately large fraction) has been accumulated. We now describe an alternative approach that does not require measurement of the total flux. Fig. 2 demonstrates that the shape of $`F(\mathrm{\Delta }T)`$ is unique for each combination of $`m`$ and $`\kappa `$. For example, a curve produced by a sample of GRBs lensed by a low optical depth of high masses may have the same value $`F(\mathrm{\Delta }T=T_{0.99})`$ as a sample lensed by smaller masses having a higher optical depth. However at all $`\mathrm{\Delta }T<T_{0.99}`$ the former will have a smaller value of $`F(\mathrm{\Delta }T)`$. We have produced the average function $`F_{av}(\mathrm{\Delta }T)`$ for masses and optical depths covering the parameter space shown in Fig 3. At each combination of mass and optical depth ($`m_{true}`$ and $`\kappa _{true}`$), we make 1000 mock observations of $`F_i(\mathrm{\Delta }T)`$ each calculated from 100 GRBs. We look for the values of $`m`$ and $`\kappa `$ such that $`F_{av}(\mathrm{\Delta }T)`$ best fits each mock observation $`F_i(\mathrm{\Delta }T)`$. Our criteria for best fit is to minimise the value

$$D=max(|F_iF_{av}|).$$
(6)
This procedure provides a likelihood for measuring $`m`$ and $`\kappa `$ given the true values $`m_{true}`$ and $`\kappa _{true}`$:

$$p_{lh}(m,\kappa |m_{true},\kappa _{true})$$
(7)
The construction of $`F_{obs}(\mathrm{\Delta }T)`$ for an observed set of GRBs then provides estimates of $`m`$ and $`\kappa `$ through application of Bayes’ theorem:

$$p(m,\kappa )=Np_{lh}(m,\kappa |m^{},\kappa ^{})p_{prior}(m^{},\kappa ^{})𝑑m^{}𝑑\kappa ^{},$$
(8)
where $`p_{prior}(m,\kappa )`$ is the assumed prior probability for $`m`$ and $`\kappa `$, and $`N`$ is a normalising constant.

To demonstrate the statistical uncertainty we have plotted in Fig. 4 the mean and variance of retrieved values of $`m`$ and $`\kappa `$ verses their assumed values ($`\kappa _{true},m_{true}`$) at fixed mass and optical depth respectively (ie. sections through $`p_{lh}`$). As in the previous calculation, the figure demonstrates that a sample of 100 GRBs is sufficient to obtain a consistent result. The systematic bias (which is smaller than the statistical uncertainty) in the recovered values is due to the finite grid of $`m`$ and $`\kappa `$ over which $`F_{av}`$ and $`p_{lh}`$ are computed. This calculation would obtain upper limits on the true values of $`m`$ and $`\kappa `$ due to possible intrinsic spread of the GRB. Also, the binsize and mass are linearly related allowing scaling of the case presented here.

## 4 Discussion

If gamma ray bursts have a cosmological origin then they provide a unique opportunity to study the cosmological abundance of compact objects through the identification of lensed pairs. However there are several mechanisms that may lead to the mis-classification of lensed pairs as spatially coincident but independent bursts. In addition, the identification of a lensed pair may be ambiguous if the relative delay is smaller than the event duration. Upper limits on the contribution to the mass density of compact objects determined from a lack of lensed pairs may therefore be underestimated.

We have demonstrated that simultaneous upper limits on the average optical depth and mass of compact objects can be estimated by measuring the average cumulative flux as a function of time for a collection of bursts and without resolving the individual lensed images. By selecting sub-samples of bursts, specific populations of compact objects could be examined. For example, a collection of bursts that are spatially coincident with fore ground galaxies (like GRB 971214 (Kulkarni et al. 1998; Diercks et al. 1999)) will probe the compact object population in galactic halos. In particular, if a population of short duration (micro-second) cosmological gamma ray bursts are identified, limits will be placed on the abundance of stellar mass compact objects.

## Acknowledgements

This work was supported by NSF grant AST98-02802. JSBW acknowledges the support of an Australian Postgraduate award and a Melbourne University Overseas Research Experience Award. We would like to thank Rachel Webster for numerous helpful discussions.


format: text

text:
# Distance Constraints for High Velocity Clouds from Optical Emission Lines

## 1. Introduction

High velocity clouds are detected primarily in H i; no resolved optical counterparts such as stars have been detected in HVCs, and the only distance upper limits are for two HVCs seen in absorption against background halo stars. The lack of distance constraints makes the nature of and models for HVCs extremely uncertain; see the review of Wakker & van Woerden (1997). Currently favored models include recycling of disk gas through a fountain (e.g. Bregman 1980); stripping from Galactic satellites; and infall of possibly primordial gas (e.g. the Local Group model of Blitz et al. 1999). These models place HVCs at from $`<10`$ kpc to $`1`$ Mpc respectively, a range of 100 in distance and $`10^4`$ in gas mass.

Faint, diffuse optical recombination emission lines are observed from some HVCs. H$`\alpha `$ emission must be caused by ionization on the cloud, either from photo-ionization by Lyman continuum radiation, or another process such as collisional ionization. Measurements of H$`\alpha `$ flux can constrain HVC distances, given models for the ionization processes.

## 2. Observations

Observing diffuse emission lines from HVCs is difficult because the emission is faint, night sky emission lines are strong and fluctuating, and HVCs are degrees across, larger than the field of view of most optical telescopes. Fabry-Perot observations of diffuse emission provide large collecting solid angle with moderate to high spectral resolution, needed to separate HVC emission from sky lines. FPs have previously been used to detect H$`\alpha `$ emission from four HVCs and the Magellanic Stream (Weiner & Williams 1996, Tufte et al. 1998, Bland-Hawthorn et al. 1998). Chopping by several degrees between object and sky fields is necessary to obtain decent sky subtraction.

Figure 1 shows example spectra from the Wide Field Camera F-P at the Las Campanas 2.5-m, with a field of view of 25 and etalon FWHM 1.2 Å ($`R=5500`$). The top panel is a extracted spectrum of a 900 sec exposure of HVC 343+32–140, one of the HVCs brightest in H$`\alpha `$, before subtraction of a sky-field spectrum. The second panel is after a two-step sky subtraction: sky lines in the object-field and sky-field spectra are fit and subtracted, then the sky-field continuum is subtracted from the object-field continuum. The H$`\alpha `$ flux is strong, 1060 milli-Rayleighs (mR), as is \[N ii\] 6583, with \[N ii\]/H$`\alpha `$ = 1.1. (1 Rayleigh = $`10^6`$ photons cm-2 s-1 into $`4\pi `$.) The lower two panels show two HVCs very faint in H$`\alpha `$, HVC 165–43–280 (41 mR) and HVC 230+61+165 (48 mR). The H$`\alpha `$ detections agree well with the H i velocities. No \[N ii\] is detected. Note the tremendous difference in strength of HVC H$`\alpha `$ and night-sky lines.

Figure 2 compiles our results from the LCO FP and the Rutgers FP at the CTIO 1.5-m, and HVC H$`\alpha `$ detections from Tufte et al. (1998) and Bland-Hawthorn et al. (1998). There is a wide range of HVC H$`\alpha `$ intensity, but clouds in the same complex tend to have similar intensities, which suggests that the variations between complexes are due to HVC properties (e.g. distances) rather than extrinsic variations (e.g. spatial variations in the ionizing field escaping from the Galaxy). On the other hand, the Magellanic Stream points vary widely: some points have strong emission while others have weak or no H$`\alpha `$ despite a high H i column density. Strong emission in the MS is often located near cloud edges.

A fundamental result of our H$`\alpha `$ survey is that we have not found any cloud near or below our photon-limited noise (generally $`1020`$ mR at 2$`\sigma `$). The faintest H$`\alpha `$ detections, at 41–50 mR, are well above our noise limit. We always either detect H$`\alpha `$ (15 of 20 HVCs), or are clobbered by residuals from sky-subtraction when HVCs are not well separated from sky lines, raising the detection limit.

## 3. Models for the source of the emission

If H$`\alpha `$ emission from HVCs is due to photo-ionization by flux from the Galaxy, we can infer distances to the HVCs. We use the HVCs A and M with known distance brackets (4–10 kpc and $`<5`$ kpc; van Woerden et al. 1999, Danly et al. 1993) and H$`\alpha `$ fluxes (Tufte et al. 1998) to normalize a model for the ionizing flux escaping from the Galaxy. Figure 3 shows contours of the ionizing photon flux $`F_{LC}`$ in the model; it has a total ionizing luminosity of $`L_{LC}=2.7\times 10^{53}`$ photons s-1, distributed in an exponential disk, and models the Galactic absorbing layer with a one-sided face-on optical depth to ionizing photons of $`\tau =2.5`$, yielding an overall, angle-averaged escape fraction $`<f_{esc}>=2\%`$. (See also Bland-Hawthorn & Maloney 1999 and Bland-Hawthorn, these proceedings.)

The inferred distances of HVCs are indicated; the error bars are for systematic variations by a factor 1.5 up or down in H$`\alpha /LC`$ ratio (statistical errors on the fluxes are much smaller). The brightest clouds are within a few kpc of the Galactic plane but fainter clouds are well away from the plane. These clouds at $`|z|>10`$ kpc are inconsistent with a Galactic fountain origin, especially given their high velocities.

However, the observed H$`\alpha `$ in the Magellanic Stream cannot be explained by this model; MS II and MS IV are much too bright compared to HVCs A and M. In fact, at a distance of $`D50`$ kpc, the MS H$`\alpha `$ emission cannot be explained by a reasonable model of photoionization from the Galaxy. The H$`\alpha `$ photon flux is $`F(H\alpha )=0.46F_{LC}0.46L_{LC}/(4\pi D^2)e^{\tau \mathrm{csc}b}`$, so that H$`\alpha `$ of 300–600 mR at MS II requires $`\tau =0.30`$ and $`<f_{esc}>=60100\%`$. The required escape fraction is unrealistically large (even if $`L_{LC}`$ were increased by $`\times 2`$) since most ionizing radiation must be absorbed in the Galaxy to power H ii regions. (In agreement with the authors, we find that Bland-Hawthorn & Maloney 1999 overestimated the ionizing flux incident on the MS.) Furthermore, there are strong spatial variations in the MS H$`\alpha `$, shown by the large dispersion of points in Figure 2. Photoionization from the Galaxy should produce roughly the same H$`\alpha `$ flux on any H i surface optically thick to ionizing radiation, since the ionizing photons travel nearly along our line of sight to the MS. Factors of 2–3 variation between different points are tolerable either from variations in the incident field or geometry of the absorbing surface, but factors of 30 are not.

A second source of ionization in the Galactic halo is needed to explain the MS H$`\alpha `$. Simply put, why are some points in the Magellanic Stream brighter than Complexes A and M, despite being much further away? The only likely source is collisional ionization, probably from interaction with hot halo gas through ram pressure and turbulent mixing/thermal conduction (Weiner & Williams 1996). The MS radial velocities are large and the space velocities with respect to halo gas are probably above 200 km s-1; energy input from ram pressure goes as $`v^3`$, so it may explain the brightness of the MS with respect to HVCs A and M, but it is barely adequate for reasonable halo gas densities (Weiner & Williams 1996); perhaps the Magellanic Stream clouds are colliding with each other. Collisional ionization is required, but the brightest MS points remain a puzzle.

Recent FUSE detections of O vi absorption in the Magellanic Stream and in a few HVCs provide further evidence for collisional ionization in HVCs (Sembach et al. 2000). O vi can be produced collisionally in hot gas at $`T10^{5.5}`$ K, which could arise from the interaction of an H i cloud with hot halo gas, but O vi cannot realistically be produced in HVCs by photoionization from O stars.

## 4. Conclusions

Given the existence of a second source of ionization in the Galactic halo, simple photoionization models may not yield reliable distances. However, robust predictions can still be made for a relation of H$`\alpha `$ flux to distance.

HVCs which are bright in H$`\alpha `$ (as HVCs go) are likely only a few kpc from the Galactic plane. The brightest HVCs plotted in Figure 2 (L, GCP, CR) have inferred distances which give them low “deviation velocities” (Wakker & van Woerden 1991), so they are not too far from Galactic differential rotation, despite having high $`V_{LSR}`$; these HVCs could well be recycled disk gas. High \[N ii\]/H$`\alpha `$ ratios suggest high-metallicity gas photoionized by a hardened radiation field, as in the extraplanar diffuse ionized gas of NGC 891 (Rand 1998).

Conversely, if an HVC is very far from the Galaxy, the only source of ionization is the metagalactic UV background. Then the HVC should be very faint in H$`\alpha `$, approaching the upper limits on H$`\alpha `$ emission from isolated extragalactic H I clouds: $`<20`$ mR at $`2\sigma `$ from the Giovanelli-Haynes cloud (Vogel et al. 1995) and $`<30`$ mR at $`2\sigma `$ from the Leo Ring (Donahue et al. 1995).

We do not find any HVCs that are very faint, even though our photon-limited noise is generally 10–20 mR (2$`\sigma `$): we either detect H$`\alpha `$, or are limited by sky subtraction residuals. The faintest HVCs we find are 165–43–280, 47–52–129, 227–34+114, and 230+61+165, with H$`\alpha `$ emission at 41–48 mR. In the photoionization model of Figure 3 their distances are $`D=15,40,11`$, and 26 kpc from the Sun. We can push the distances to the limit by considering solely collisional ionization: the HVC H$`\alpha `$ fluxes are some 5-15 times smaller than the fluxes typical of the Magellanic Stream. If we attribute the H$`\alpha `$ to interaction with hot halo gas, and assume that the hot gas density falls off as $`r^2`$ to $`r^3`$, the HVCs can be at most 2 to 4 times farther than the MS, or 100-200 kpc.

At distances of somewhere between 20 and at most 200 kpc, the H$`\alpha `$-faint HVCs are not Local Group objects as defined by the Blitz et al. (1999) model. However, their low \[N ii\]/H$`\alpha `$ ratios suggest low metallicity, and they are not easily produced by a Galactic fountain: for example, HVC 165–43–280, with $`V_{GSR}=240`$ km s-1, has faint H$`\alpha `$, is at least 10 kpc away (20 kpc Galactocentric radius), and has a large mass and kinetic energy; perhaps it is infalling tidal debris. The best explanation for these H$`\alpha `$-faint HVCs is still likely to be gas, of uncertain origin, accreting onto the Galaxy.

## References

Bland-Hawthorn, J. & Maloney, P.R. 1999, ApJ, 510, L33

Bland-Hawthorn, J., Veilleux, S., Cecil, G.N., Putman, M.E., Gibson, B.K. & Maloney, P.R. 1998, MNRAS, 299, 611

Blitz, L., Spergel, D.N., Teuben, P.J., Hartmann, D. & Burton, W. B. 1999, ApJ, 514, 818

Bregman, J.N. 1980, ApJ, 236, 577

Danly, L., Albert, C.E. and Kuntz, K.D. 1993, ApJ, 416, L29

Donahue, M., Aldering, G., & Stocke, J.T. 1995, ApJ, 450, L45

Rand, R.J. 1998, ApJ, 501, 137

Sembach, K.R., et al. 2000, preprint astro-ph/0005012

Tufte, S.L., Reynolds, R.J. & Haffner, L.M. 1998, ApJ, 504, 773

van Woerden, H., Schwarz, U.J., Peletier, R.F., Wakker, B.P. & Kalberla, P.M.W. 1999, Nature, 400, 138

Vogel, S.N., Weymann, R. , Rauch, M. & Hamilton, T. 1995, ApJ, 441, 162

Wakker, B.P. & van Woerden, H. 1991, A&A, 250, 509

Wakker, B.P. & van Woerden, H. 1997, ARA&A, 35, 217

Weiner, B.J., & Williams, T.B. 1996, AJ, 111, 1156


text:
# Phase diagram of the quarter–filled extended Hubbard model on a two–leg ladder

## I Introduction

Quasi–one–dimensional systems present a unique opportunity to study the interplay between strong quantum fluctuations on the one hand and a tendency to charge or spin ordering on the other. Examples include Peierls or spin–Peierls behavior as well as charge ordering due to electron–electron interaction. One material in which such behavior has been found is the two–leg ladder material NaV2O5. NaV2O5 undergoes a phase transition at $`T_c=34`$ K that is characterized by the opening of a spin gap and a doubling of the unit cell. Although this transition was originally thought to be spin–Peierls, recent studies have found evidence for charge order. Above the transition, the material seems to be best described as a quarter–filled ladder. Such two–leg–ladder structures have also been found in a number of other materials, including the Vanadates MgV2O5 and CaV2O5, and the cuprates SrCu2O3 and Sr14Cu24O41. For a more detailed description of ladder materials and models as well as a discussion of the extensive theoretical work on ladder models, we direct the reader to a recent review  and the references contained therein.

In NaV2O5, the character of the charge ordering and the nature of the transition are currently under debate. While evidence is growing  that the charge–ordered state has a zigzag–like charge arrangement, the question of whether this charge ordering is continuous or discontinuous has not yet been settled experimentally. Theoretically, the magnetic properties of the compound for $`T<T_c`$ have been calculated for different charge–ordering patterns  and compared to experiments. Again, zigzag charge ordering on ladders (possibly with intervening disordered ladders ) leads to magnon dispersions in good agreement with neutron scattering data. Additionally, the possibility of two separate transitions at nearby temperatures has been raised; in this scenario the charge order would set in at one temperature while the spin gap would open at a second, lower critical $`T`$. Very recently, a scenario of singlet cluster formation instead of zigzag charge ordering has been suggested.  However, this model seems not to be able to reproduce the experimentally found spin gap and magnon dispersion data.  It is clear that more experimental work is necessary to clarify the nature of the low–temperature phase of NaV2O5.

One of the simplest models of interacting electrons that allows for charge ordering is the extended Hubbard model, i.e. the Hubbard model supplemented by an additional nearest–neighbor (NN) Coulomb repulsion, $`V`$. This model has been studied in one dimension (1D) in the strong–coupling limit, at quarter filling  and at half filling  and in between, in the 2D system at half filling,  and within the Dynamical Mean Field Theory (the limit of infinite dimensions) at quarter  and half filling.  A variety of techniques, such as mean–field approximations, perturbation theory, as well as numerical methods as quantum Monte Carlo and the Density Matrix Renormalization Group (DMRG) have been employed.

The result of investigations of the charge–order transition can be summarized as follows: At the mean–field level, the transition between a homogeneous state and a charge–density wave (CDW) state at half filling in a hypercubic lattice occurs at $`V_c=U/z_0`$, where $`z_0`$ denotes the number of nearest neighbors ($`z_0=2d`$) and $`U`$ is the on–site interaction. Numerical studies  indicate a slightly higher value of $`V_c`$, at least in 1D. Interestingly, the transition at half filling in 1D has been found to be second order at small $`U/t`$ and first order at large $`U/t`$ with the tricritical point located at $`U_c/t46`$.  Here we use the term “first order” to denote discontinuous behavior of the charge order parameter as a function of microscopic parameters such as $`V`$ or band filling, and “second order” to denote continuous behavior. For fillings below half–filling in 1D, the situation is more complicated because a number of phases compete at large $`V`$.  For dimension larger than one, indications are that the charge–order transition is generally first order.  However, conclusive studies that can reliably distinguish between first– and second–order transitions are lacking. At small $`U`$ and large $`V`$, the extended Hubbard model undergoes phase separation (PS) rather than a transition to a CDW state. For the 1D model between quarter– and half–filling, it has been established that PS occurs for $`|U|/t<4`$ in the $`V=\mathrm{}`$ limit, whereas for $`U/t>4`$ the system undergoes a transition to a $`q=\pi `$ CDW state for sufficiently large $`V`$. Phase separation in higher dimensions has also been discussed.

For the 1D system in particular there has been recent interest in the possibility of dominant superconducting correlations in the uniform ground–state away from half–filling when $`VUt`$,  i.e., in the proximity of the phase–separated region. We note here that the uniform phase in 1D off half–filling is metallic and can in general be described within the Luttinger–liquid picture. Although dominant superconducting correlations have not been established in the ground state of the 1D extended Hubbard model to date, a number of non–Luttinger–liquid effects have been observed.

Some of the present authors have previously studied  the charge–order transition in the extended Hubbard model on the two–leg ladder at various band fillings for $`U/t=4`$ and 8. A transition to a checkerboard charge–ordered state was found for all fillings between quarter– and half–filling. The transition is second–order near quarter filling and first–order near half–filling for sufficiently large $`U`$.

The focus of the present paper is on this model at quarter filling, with a two–fold purpose: First, we present a comprehensive study of the phase diagram as a function of $`U/t`$ and $`V/t`$ for repulsive $`U`$ and $`V`$ and discuss the properties of the ground state phases as well as the nature of the phase transitions. Second, we derive an effective Hamiltonian for the spin degrees of freedom in the charge ordered state at strong coupling, and compare the predictions of this effective low–energy theory with our numerical results.

### A Phase diagram

Our main results, the phase diagrams deduced from the numerical calculations, are summarized in Figs. 1 (isotropic hopping) and 2 (anisotropic hopping). The phases are distinguished by the presence or absence of a gap for spin and/or charge excitations. To denote this, we employ the following labeling: HI$`_{\text{sg}}`$ denotes a homogeneous insulator (nonzero charge gap) with a spin gap, HI a homogeneous insulator without a spingap, HM$`_{\text{sg}}`$ (HM) a homogeneous metallic phase having zero charge gap with (without) spin gap, and CDW$`_{\text{sg}}`$ (CDW) is a charge–ordered state with (without) spin gap. The CDW states are always insulating in the present quarter–filled model. The phase diagrams can be roughly divided into four regions: (i) Weak coupling: for small $`U`$ and $`V`$ we find homogeneous phases similar to the ones in the “bare” Hubbard model (see discussion in Sec. II and results in Sec. III B). (ii) Large $`U`$, small $`V`$: These homogeneous strong coupling phases have characteristics similar to the weak–coupling region. (iii) Small $`U`$, large $`V`$: phase separation, this is discussed further in Sec. III D. (iv) Strong coupling: large $`U`$ and $`V`$ lead to an insulating checkerboard charge–ordered with either gapless or gapped spin excitations depending on the ratio of $`V/U`$.

For isotropic hopping we have numerically determined the phase boundaries as shown in Fig. 1. The precise location of the CDW–CDW$`_{\text{sg}}`$ boundary at which the spin gap closes at intermediate coupling could not be obtained by the methods used here; we have indicated this uncertainty by a question mark in the phase diagram.

In the case of anisotropic hopping we have not mapped out the full phase diagram, but the numerical results (discussed in Sec. III) provide the schematic pictures shown in Fig. 2. Varying the ratio of the rung to leg hopping strengths, $`t_{}/t_{}`$, has two effects: (a) For small $`t_{}/t_{}`$ there appears a metallic phase (HM$`_{\text{sg}}`$) with spin gap and dominating d–wave–like singlet pair correlations (as in the “bare” $`V=0`$ Hubbard ladder). (b) The existence of a spin gap depends strongly on the hopping ratio, i.e., there is a transition as function of $`t_{}/t_{}`$ where a spin gap opens (in both the homogeneous and CDW phases). The critical hopping ratio may depend on the interaction strength, but is near unity in the homogeneous phases for reasonable values of the interactions.

The rest of the paper is organized as follows: In Sec. II, we introduce the extended Hubbard model and discuss some results known for the $`V=0`$ case, i.e., the “bare” Hubbard model on a two–leg ladder. In Sec. III, we present our numerical results, discuss the properties of the phases shown in Figs. 1 and 2, and examine the transition to the charge–ordered state. At large $`V`$, where charge order is well established, it is possible to derive an effective Hamiltonian for the residual spin degrees of freedom; this is done in Sec. IV. A summary and a discussion of the relevance of our results to experimental systems (especially NaV2O5) terminates the paper.

## II Extended Hubbard model

The single–band extended Hubbard model has the Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \underset{ij\sigma }{}}t_{ij}(c_{i\sigma }^{}c_{j\sigma }+h.c.)`$ (2)
$`+U{\displaystyle \underset{i}{}}n_in_i+{\displaystyle \underset{ij}{}}V_{ij}n_in_j.`$
Here we consider a lattice consisting of two chains of length $`L`$, i.e., a ladder, and restrict ourselves to the band filling $`n=N/(2L)=1/2`$, where $`N`$ is the number of electrons. The summation $`ij`$ runs over all pairs of nearest–neighbor sites on the ladder, taking open boundary conditions between the chains. The hopping matrix elements along the legs and rungs of the ladder are denoted $`t_{}`$ and $`t_{}`$, respectively, and the nearest–neighbor Coulomb interactions are similarly denoted $`V_{}`$ and $`V_{}`$. Unless otherwise noted, we will use $`t_{}`$ as unit of energy. In this work, we will treat primarily the “isotropic” case, $`t_{}=t_{}=t`$ and $`V_{}=V_{}=V`$.

The non–interacting Hamiltonian ($`U=V=0`$) can be diagonalized by a Fourier transform (for periodic boundary conditions in the chain direction), leading to the single–particle energies

$$ϵ_𝐪=2t_{}\mathrm{cos}q_x+t_{}\mathrm{cos}q_y$$
(3)
where $`𝐪=(q_x,q_y)`$, $`q_x`$ is the momentum along the chains and the momenta $`q_y=0,\pi `$ correspond to bonding and antibonding symmetry, respectively. Either one or both of the bands can be occupied in the noninteracting system, depending on the total particle density and the ratio of $`t_{}`$ and $`t_{}`$. At quarter filling, the transition occurs at isotropic hopping: for $`t_{}<t_{}`$ both bands are less than half–filled, whereas for $`t_{}>t_{}`$ the bonding band is half–filled and the antibonding band is unoccupied.

The effect of the Hubbard interaction $`U`$ on this system has been extensively studied. In the weak–coupling limit, $`Ut`$, the phase diagram has been investigated using the perturbative renormalization group (RG).  A variety of phases have been shown to exist as a function of band filling and hopping anisotropy. In general, the two–band system can have four possible modes (symmetric and antisymmetric charge and spin modes), each of which can be either massive or massless. The phases can therefore be classified using the notation C$`n`$S$`m`$ where $`n`$ and $`m`$ designate the number of gapless charge and spin modes, respectively ($`0n,m2`$). At quarter filling, the weak–coupling RG  for the “bare” Hubbard model ($`V=0`$) yields the following results: For $`t_{}>t_{}`$ the system behaves as a half–filled Luttinger liquid; Umklapp scattering in the bonding channel is a relevant perturbation which leads to a charge–gapped C0S1 phase at small $`U`$. In contrast, deep in the two–band region, $`t_{}t_{}`$, one finds a metallic C1S0 phase in the weak–coupling limit. Near isotropic hopping, the bottom of the antibonding band just “touches” the Fermi surface, and the curvature of the dispersion becomes important, leading to additional narrow regions of C2S2 and C2S1 phases. Several of these weak–coupling predictions have been verified by numerical DMRG calculations in the intermediate and strong coupling regimes for a wide range of filling.  Systematic studies of the phases of the extended Hubbard model on a two–leg ladder away from half–filling have to our knowledge not yet been carried out.

## III Numerical results

In this section we present the results of our numerical investigations and discuss the characteristics of the phases shown in Figs. 1 and 2.

### A Technique and observables

The numerical results have been calculated with the DMRG technique  on lattices of up to $`2\times 80`$ sites with open boundary conditions at the ends of the chains as well as between the two chains. Most data shown are obtained by keeping 600 states per block, resulting in the sum of the discarded density matrix eigenvalues being typically $`10^8`$ or less. For small system sizes we have checked the convergence by using up to 1000 states per block. Unless otherwise noted, we estimate the errors in the gap energies and correlation functions obtained using the DMRG procedure to be less than a few percent.

Important ground–state properties are the static charge and spin correlation functions: we have calculated the static charge structure factor

$`C(𝐪)={\displaystyle \frac{1}{2L}}{\displaystyle \underset{i}{}}e^{i𝐪𝐑_i}\overline{C}(𝐑_i)`$ (4)
where

$$\overline{C}(𝐑_i)=\frac{1}{N_{\mathrm{av}}}\underset{\{j\}}{}\delta n_{j+i}\delta n_j,$$
(5)
$`\mathrm{}`$ denotes the ground–state expectation value, $`\delta n_j=n_jn_j`$, and we average over typically $`N_{\mathrm{av}}=6`$ sites to remove oscillations due to the open boundaries. The spin structure factor $`S(𝐪)`$ is defined similarly in terms of the spin–spin correlation function $`S_{j+i}^zS_j^z`$.

The nature of the low–lying excitations can be determined by calculating the energy gaps of the system. In particular, we will consider the charge and spin gaps, defined as

$`\mathrm{\Delta }_c`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[E_0(L,N+2)+E_0(L,N2)2E_0(L,N)\right],`$ (6)
$`\mathrm{\Delta }_s`$ $`=`$ $`E_0(L,N,S_z=1)E_0(L,N,S_z=0)`$ (7)
where $`E_0(L,N)`$ is the ground state energy of ladder system with $`2L`$ sites and $`N`$ electrons. Since we calculate the gaps $`\mathrm{\Delta }(L)`$ on finite systems, the gaps must be extrapolated to $`L\mathrm{}`$; we do this by performing a polynomial fit in $`1/L`$ through the data points from the larger system sizes ($`L24`$). Although we include both $`1/L`$ and $`1/L^2`$ terms, the coefficient of the quadratic term is quite small in most cases. The uncertainty of the extrapolated value depends strongly on finite–size effects which become large when the correlation length becomes large; the results for the gaps are most accurate in the strong coupling region, $`U,Vt`$, and for not too small $`t_{}`$.

### B Homogeneous phases

First, we concentrate on the states without charge order, i.e., the region of small $`V`$ as shown in Fig. 1. The calculated charge and spin correlation functions (see e.g. Fig. 1 of Ref. ) indicate antiferromagnetic correlations peaked at ordering vector $`(\pi ,\pi )`$. The nature of the phases is best probed by calculating spin and charge gaps. The extrapolation to the thermodynamic limit is illustrated in Fig. 3, in which we show results for charge and spin gaps at $`U/t_{}=8`$, $`V/t_{}=2`$ and different values of the hopping ratio $`t_{}/t_{}`$. As noted above, finite–size effects increase with decreasing interchain coupling $`t_{}`$. However, we have verified (by keeping more DMRG states per block and/or using a fit with a $`1/L`$ term only) that, within the numerical accuracy available, $`\mathrm{\Delta }_s`$ in Fig. 3 vanishes for $`t_{}/t_{}=1`$, but is nonzero for $`t_{}/t_{}=0.7`$.

The extrapolated DMRG results for charge and spin gaps in the homogeneous phase are displayed in Fig. 4. First we discuss the $`V=0`$ case, i.e., the “bare” Hubbard ladder. According to the weak coupling RG,  both spin and charge gaps vanish in the case of isotropic hopping. Varying the hopping ratio $`t_{}/t_{}`$ tunes the system through the one–band to two–band transition: For $`t_{}>t_{}`$ Umklapp scattering opens a charge gap (C0S1 phase); with decreasing $`t_{}/t_{}`$ one finds narrow regions of C2S2 and C2S1, followed by a C1S0 phase. Our results for $`U/t=8`$ agree with these predictions: we find a C0S1 (HI) phase for $`t_{}>t_{}`$ and a C1S0 (HM$`_{\text{sg}}`$) phase for $`t_{}<0.9t_{}`$. The data also indicate a narrow region where both $`\mathrm{\Delta }_c=0`$ and $`\mathrm{\Delta }_s=0`$ (HM), possibly corresponding to the C2S2 and C2S1 phases. Turning to $`V>0`$, the described behavior continues to small nonzero values of $`V`$, but the transition points shift to smaller $`t_{}/t_{}`$. A further increase of $`V`$ suppresses the metallic phase, and only the spin gap transition remains. Data for $`V/t_{}=2`$ is shown in Fig. 4: the behavior of the spin gap is similar to the $`V=0`$ case, i.e., it is finite for small $`t_{}`$ and vanishes for $`t_{}/t_{}`$ larger than some critical value. However, the charge gap is found to be nonzero for all hopping ratios examined here (see also Figs. 6 and 7 below).

These data can be understood from the RG analysis of Ref. : the additional nearest–neighbor repulsion does not introduce a new relevant operator, but only changes the scaling dimension of the perturbation introduced by $`U`$. This implies that small $`V`$ does not modify the phases found at $`V=0`$. However, our data indicate that relatively small values of $`V`$ are enough to drive the system to an insulating state even for $`t_{}<t_{}`$.

### C CDW phases and charge ordering transition

As the nearest–neighbor repulsion $`V`$ is increased, we expect a transition to a checkerboard charge–ordered state. As this has been examined in our earlier paper,  here we summarize the main findings: At large $`V`$, an insulating CDW state with ordering wavevector $`𝐐=(\pi ,\pi )`$ occurs for all fillings between quarter and half–filling. At quarter–filling, the transition is second order, i.e., the order parameter $`\eta =lim_L\mathrm{}C(𝐐)/n^2`$ vanishes continuously upon approaching a critical $`V_c(U)`$ from above. Interestingly, the transition has been found to change from second–order to first–order at higher band filling as a function of $`U/t`$;  such tricritical behavior has also been observed in the 1D case at half–filling. In the quarter filled CDW state, the spin correlations indicate zigzag antiferromagnetic ordering of the spins $`\frac{1}{2}`$ on the occupied sites; at larger filling the spin correlations become incommensurate and are gradually suppressed.

We now turn to our new results and a more detailed discussion of the quarter–filled system. In Fig. 5, the finite–size scaling for $`\eta `$ as well as the extrapolated values for $`\sqrt{\eta }`$, which corresponds to the relative difference of the sublattice occupancies in the broken–symmetry charge–ordered state, are shown for a typical second–order transition. (Data for first–order transitions at larger filling are shown in Ref. , Figs. 2 and 3.) The inclusion of the quadratic fit term turns out to be important near the transition. The $`V_c(U)`$ values obtained from the numerics are displayed as phase boundary in Fig. 1. Note that $`V_c`$ decreases with increasing $`U`$ in the quarter–filled case treated here, similar to the behavior found in 1D. Our results for large $`U`$ suggest that $`V_c`$ is nonzero in the $`U\mathrm{}`$ limit as in 1D; an extrapolation based on data up to $`U/t=64`$ yields $`V_c(U=\mathrm{})2t`$. It is interesting to contrast this with half filling, for which weak– and strong–coupling approximations  as well as numerical studies  yield $`V_cU/z_0`$ for a hypercubic lattice with $`z_0`$ being the number of nearest neighbors, i.e., in the half–filled case $`V_c`$ increases with $`U`$.

The behavior of the low–energy electronic excitations in the vicinity of $`V_c`$ provides further information on the character of the charge–order transition. Since the energy gaps show different behavior for different values of the hopping ratio $`t_{}/t_{}`$, as seen in Fig. 4, we focus on two representative values of the hopping anisotropy.

The energy gaps as function of $`V/t_{}`$ together with the order parameter at $`U/t_{}=8`$ are shown in Figs. 6 and 7 for $`t_{}/t_{}=0.7`$ and 1.4. The main observation is that no gap opens or closes at the transition to the charge–ordered state. Small $`t_{}/t_{}`$ leads to fully gapped C0S0 phases on both sides of the transition (HI$`_{\text{sg}}`$-CDW$`_{\text{sg}}`$ transition), whereas at large $`t_{}/t_{}`$ both the homogeneous and the charge–ordered phase have zero spin gap and nonzero charge gap (HI-CDWtransition). For the former case (Fig. 6), the spin gap decreases appreciably as the charge–ordered state is entered. This indicates a change in the spin dynamics from a “weak–coupling” regime dominated by band effects to a “strong–coupling” regime determined by the physics of a frustrated spin chain. This behavior will be discussed in detail in Sec. IV.

In any case, it appears that low–lying fermionic excitations do not play a role in the critical dynamics at the charge–order transition. This implies that this zero–temperature transition must be in the universality class of the 1+1–dimensional Ising model. Note, however, that inter–ladder couplings will increase the effective dimensionality of this transition in experimental ladder systems at low enough energies or temperatures.

### D Phase separation

For small $`U`$ and large $`V`$, phase separation is expected: In the $`V\mathrm{}`$ limit, existing double occupancies are immobile and cannot be broken up, whereas single fermions can move in the “unoccupied” space. For $`U<U_{\mathrm{PS}}(V)`$, the system then separates into a region with double occupancies at every second site (i.e., checkerboard order with two electrons per occupied site) and a region in which the other electrons can gain kinetic energy by hopping. For the one–dimensional system it is possible to solve the $`V=\mathrm{}`$ problem exactly because it maps onto non–interacting spinless fermions moving on open chain segments.  While the mapping to spinless fermions is similar for the ladder system, the geometry leads to interactions among the fermions which preclude an exact solution. Nevertheless, the qualitative behavior should be similar to the 1D case, i.e., for $`V=\mathrm{}`$ there should be a critical $`U_{\mathrm{PS}}(\mathrm{})`$ ($`U_{\mathrm{PS}}(\mathrm{})=4t`$ in 1D) below which the system phase separates. For small $`V`$, the phase separation should disappear.

The numerical results for charge and spin correlation functions at small $`U`$ are shown in Fig. 8. Incommensurate peaks appear in both the $`q_y=\pi `$ channel of $`C(𝐪)`$ and the $`q_y=0`$ channel of $`S(𝐪)`$ as $`V`$ is increased. At the largest value shown, $`V/t=8`$, there are strong oscillations and side peaks in $`C(𝐪)`$, an indication of PS.

To examine the thermodynamic stability of the system, we have numerically computed the compressibility of the system which is defined as

$$\kappa =\frac{4L}{N^2}\left[E_0(L,N+2)+E_0(L,N2)2E_0(L,N)\right]^1.$$
(8)
Our results clearly show the occurence of phase separation in the large $`V`$, small $`U`$ region indicated by (i) a diverging compressibility, (ii) oscillating incommensurate spin and charge correlations with wave vectors strongly dependent on the system size, and (iii) the occurence of site charge densities greater than unity. Note that at quarter filling no double occupancies occur even in the perfectly charge–ordered state. The appearance of doubly occupied sites in the phase–separated state is clearly consistent with the phase–separation mechanism explained above. The criteria (i)–(iii) give consistent results and allow for a reasonably accurate determination of the PS boundary (see Fig. 1), even though finite–size effects in the calculation of the compressibility are large.

In contrast to that found in the 1D chain, the phase–separation boundary has non–monotonic behavior, i.e. $`U_{\mathrm{PS}}(V)`$ shows a maximum at around $`V/t8`$, $`U/t2.4`$, as can be seen in Fig. 1. The described “re–entrant” (non–monotonic) behavior of the PS boundary is illustrated in Fig. 9 in which we display the charge gap as a function of $`V`$ for $`U/t=2.2`$. Here we find a homogeneous phase at small $`V`$, a charge–ordered phase at large $`V`$, and a region of phase separation in between, for $`4.8<V/t<13.5`$.

Another difference with the behavior of the single–chain extended Hubbard model is that no homogeneous phase appears for large $`V`$ and small or intermediate values of $`U`$: increasing $`U`$ in the PS region drives the system directly into the charge–ordered state (Fig. 1). In contrast, in the 1D system, a homogeneous phase is present at any $`V`$ and the boundaries to the charge–ordered and to the PS phases merge (at $`U/t=4`$) only in the $`V\mathrm{}`$ limit.

For the present ladder system, the behavior at the boundary between the charge–ordered state and the PS region is quite interesting: The charge gap appears to vanish continuously at this boundary (Fig. 9). However, the charge–density wave order parameter $`\sqrt{\eta }`$ does not tend to zero when approaching the PS boundary from the charge–order phase. Moreover, the numerical results for small $`U`$ (in the PS region) indicate charge–density oscillations in the spatial regions without double occupancies. This suggests that strong CDW correlations exist on both sides of the PS boundary, and the transition can be interpreted as “continuous”.

## IV Spin dynamics in the strong–coupling limit

This section focuses on the low–energy spin dynamics at large $`U`$ and $`V`$ in the checkerboard charge–ordered state. The charge gap in this state can be estimated to be $`\mathrm{\Delta }_c=\mathrm{min}(U,3V)`$ when $`tU,V`$, by neglecting the kinetic energy. In this limit, each occupied site in the charge–ordered state carries charge $`e`$ and spin $`\frac{1}{2}`$, and the spin states are degenerate for $`V=\mathrm{}`$. We would like to discuss the spin ordering arising from effective exchange interactions which occur for small but finite $`t/V`$. We can do this by treating $`t/V`$ as a perturbation, in a manner similar to the derivation of the effective spin exchange in the large–$`U`$ Hubbard model at half–filling which leads to the mapping to an antiferromagnetic Heisenberg model. However, the present problem is slightly more complicated because the degeneracy is lifted in fourth order in the hopping rather than in second order as in the half–filled Hubbard model.

The aim is to find an effective Hamiltonian for the residual spin degrees of freedom. It is easy to see that this model will be a frustrated antiferromagnetic Heisenberg $`J_1`$$`J_2`$ chain where $`J_1`$ and $`J_2`$ are a diagonal (1,1) and a horizontal (2,0) coupling between the spins in the checkerboard ordered state. It is well–known  that this model has a zero–temperature phase transition as a function of $`\alpha =J_2/J_1`$. For $`\alpha <\alpha _c`$ the ground state is gapless with power–law correlations. For $`\alpha >\alpha _c`$ a spontaneous dimerization occurs which leads to a spin gap and a doubly degenerate ground state. The numerical estimate for $`\alpha _c`$ is 0.2411. At $`\alpha =1/2`$ (the Majumdar–Ghosh point), the ground state has been shown to be an exact product of nearest–neighbor singlets. Therefore, a corresponding spin–gap transition is also possible in the charge–ordered state of the $`tUV`$ ladder (i.e., a CDW – CDW$`_{\text{sg}}`$ transition) provided that the effective $`\alpha `$ can be tuned through the critical value by changing the system parameters.

We use a recently developed method  based on cumulants to derive an effective Hamiltonian for the spin degrees of freedom in the charge–ordered state of the quarter–filled model. We give the derivation in the appendix, and here state only the final result for isotropic nearest–neighbor repulsion, $`V_{}=V_{}=V`$:

$`J_1`$ $`=`$ $`{\displaystyle \frac{2t_{}^2t_{}^2}{V^2}}\left\{{\displaystyle \frac{2}{U}}+{\displaystyle \frac{2}{U+2V}}+{\displaystyle \frac{1}{V}}\right\}`$ (9)
$`J_2`$ $`=`$ $`{\displaystyle \frac{t_{}^4}{V^2}}\left\{{\displaystyle \frac{1}{U}}+{\displaystyle \frac{2}{U+2V}}\right\}.`$ (10)
The lowest–order nonzero contributions to $`J_1`$ and $`J_2`$ are indeed of order $`t^4/V^3`$ and terms of order $`t^4/(V^2U)`$ also appear. For anisotropic $`V`$, the general expressions become more complicated and are given in the appendix. It turns out that the ratio $`\alpha =J_2/J_1<1/4`$ for the isotropic case $`t_{}=t_{}`$ and $`V_{}=V_{}`$; it approaches 1/4 for $`VUt`$, as shown in Fig. 10. The plot shows that the dimerization transition will take place at $`V/U2`$. However, this transition is hard to observe numerically since the induced gap is very small as we discuss below. To access larger values of $`\alpha `$ a hopping anisotropy $`t_{}/t_{}<1`$ is necessary. The parameter $`\alpha `$ can be tuned to any value by varying the hopping ratio.

To verify the expressions for $`J_1`$ and $`J_2`$ given above, we have studied the behavior of the charge–ordered state in the strong–coupling limit for different hopping anisotropies. In order to interpret results for the spin gap, it is important to note that $`\mathrm{\Delta }_s`$ in the $`J_1`$$`J_2`$ chain vanishes exponentially near the critical point $`\alpha _c`$, leading to nonzero but very small values for $`\alpha <0.3`$. Therefore, the spin gap calculations in the charge–ordered state require an anisotropy in $`t`$ or $`V`$ in order to reach $`\alpha `$ values significantly larger than 0.3. Furthermore, they are feasible only in a window of intermediate values of $`U/t,V/t`$: overly small values do not lead to a charge–ordered state, whereas overly large values of $`V`$ lead to an unobservably small spin gap (of order $`Jt^4/V^3`$).

The nature of the magnetic ground state can also be probed using static spin correlation functions. For the $`J_1`$$`J_2`$ chain it is known  that the static structure factor $`S(q)`$ is peaked at $`q=\pi `$ for $`\alpha <1/2`$. For $`\alpha >1/2`$, the peak position shifts to smaller $`q`$ as $`\alpha `$ is increased, approaching $`\pi /2`$, the value for two uncoupled chains with a doubled lattice constant, as $`\alpha `$ becomes large.

The spin correlations at $`U=V=8t_{}`$ and various hopping ratios obtained from the DMRG calculations are shown in Fig. 11. (Note that we use $`t_{}`$ rather than $`t_{}`$ as the energy reference in this section.) A shift of the maximum from $`q_x=\pi `$ to $`q_x=\pi /2`$ (at $`q_y=0`$) with increasing $`t_{}/t_{}`$ is clearly visible. To compare quantitatively with the strong–coupling picture of the frustrated spin chain, we show in Fig. 12 the results for the spin gap $`\mathrm{\Delta }_s`$ and for the peak position $`q^{}`$ in the spin structure factor (see Fig. 11) for different parameter sets in the charge–ordered phase. Data for the frustrated spin chain from Ref.  are also shown for comparison. Note that the data are plotted as a function of the ratio $`J_2/J_1`$ with the values of these effective couplings taken from the strong–coupling expressions (10). The spin gap value follows the strong coupling prediction closely even for intermediate values of $`U/t`$ and $`V/t`$. The peak position also shows the expected behavior, i.e., it deviates from $`\pi `$ when the effective $`J_2/J_1`$ exceeds a certain value. For large $`U`$ and $`V`$, the agreement with the results from the frustrated spin chain is nearly perfect, clearly indicating that the spin dynamics in the strong–coupling charge–ordered state is correctly described by the $`J_1`$$`J_2`$ spin chain. For smaller values of $`U`$ and $`V`$, there are slight deviations in the peak position from the spin chain data: the region of incommensurate spin order becomes narrower with decreasing interaction. This might be expected because there is no incommensurability at half–filling in the non–interacting limit. A similar behavior for $`S(q)`$ has been found for the half–filled Hubbard chain with next–nearest–neighbor hopping  which can also be mapped onto an effective frustrated spin chain in the large–$`U`$ limit.

Now we turn to a discussion of the special case of isotropic hopping for which the $`V=0`$ weak coupling system is near the one–band to two–band transition. The numerical results obtained by DMRG indicate a zero spin gap and nonzero charge gap for any finite $`V`$ (outside the phase–separation region). This is in disagreement with the strong–coupling analysis presented above which predicts a CDW$`_{\text{sg}}`$ phase at large $`V/U`$ due to spontaneous dimerization. Since $`\alpha =J_2/J_1`$ is close to $`\alpha _c`$, however, the spin gap would be very small and therefore hard to observe using numerical methods. By assuming that the strong–coupling picture is also valid in the intermediate–coupling regime, we can locate the CDW – CDW$`_{\text{sg}}`$ boundary as shown in Fig. 1. Since it is not possible to deduce the behavior of the spin gap close to the charge–order transition from the current numerical results, we cannot decide whether the charge–order transition line and the spin–gap transition line meet for the case of isotropic hopping. Additional numerical approaches (e.g., based on level–crossing methods) could be used to check the spin–gap scenario and to determine the precise location of the boundary of the spin–gap phase.

It is worth pointing out that although a spin gap is present in both the homogeneous and the charge–ordered phases at small $`t_{}/t_{}`$, the mechanisms for the spin gap opening appear to be quite different: The strong–coupling CDW – CDW$`_{\text{sg}}`$ transition involves spontaneous dimerization in a spin model and is described by a sine–Gordon theory, whereas the weak–coupling case at $`V=0`$ is more complicated due to the presence of low–lying charge modes (see Ref.  for a discussion on the RG for $`V=0`$), however, not much is known about the $`V>0`$ case.

## V Conclusions

In summary, we have studied the phase diagram of the extended Hubbard model on a quarter–filled two–leg ladder. At very small $`V`$, the system behaves as in the $`V=0`$ case , while slightly larger values of $`V`$ lead to an insulating state with either zero (HI) or nonzero (HI$`_{\text{sg}}`$) spin gap depending on the hopping anisotropy. For $`U`$ and $`V`$ both strong, the ground state shows zigzag charge order. In this phase each occupied site carries spin $`\frac{1}{2}`$ and the residual kinetic energy leads to effective antiferromagnetic exchange interactions between the spins. We have rigorously established the mapping of the spin degrees of freedom to a frustrated spin–$`\frac{1}{2}`$ chain in the strong coupling limit, $`U,Vt`$. This effective spin chain can be either in the gapless regime (CDW) with algebraic spin correlations or in the spontaneously dimerized regime (CDW$`_{\text{sg}}`$) with gapped spin excitations. For $`t_{}<t_{}`$, the spin gap in the charge–ordered state could be numerically observed. Its magnitude is in good agreement with the results for a corresponding $`J_1`$$`J_2`$ spin chain down to $`U/t=4`$. The dimerization of the effective spin chain can be interpreted as bond–order wave  in the original Hubbard model, so the system has an insulating CDW$`_{\text{sg}}`$ ground state with coexisting bond–order and charge–density waves. Finally, at small values of $`U`$ and moderate to large values of $`V`$, the system phase separates into a phase of immobile double occupancies on every second site and a phase of mobile single electrons.

We have identified a purely electronic mechanism for the opening of a spin gap in a quarter–filled CDW system on a ladder based on the physics of a frustrated spin chain. However, we note here that the spin–gap physics discussed in Sec. IV probably cannot be realized in NaV2O5 since this material has $`t_{}2t_{}`$ (see Refs. ) leading to $`J_1J_2`$ for the effective spin chain. It is likely that the spin gap opening in NaV2O5 is driven by the interplay of charge ordering and phonons, as suggested in Ref. . Other effects which are important for the spin dynamics in NaV2O5 are hopping processes between the ladders which may lead to quite large exchange terms across the ladders.

We have found no evidence for “exotic” phases in the examined ladder system like the ones present in the single–chain model. There appears to be no metallic phase in the quarter–filled model except for the one at very small $`V`$ and $`t_{}<t_{}`$ (Fig. 2). This phase has been discussed  in the context of the “bare” Hubbard model on the ladder.

Any effects of interladder couplings have been neglected in the present treatment, as well as the interplay of electron and lattice effects which is known to lead to further interesting ordering effects;  these should be investigated in the future. Also, a more detailed study of the spin dynamics in a partially charge–ordered state ($`V<\mathrm{}`$) could be performed.

###### Acknowledgements.

The authors thank D. Baeriswyl, R. Bulla, P. G. J. van Dongen, R. E. Hetzel, and A. P. Kampf for useful conversations. M.V. acknowledges support by the DFG (VO 794/1–1) and by US NSF Grant No DMR 96–23181, R.M.N. has been supported by the Swiss National Foundation under Grant No 20–53800.98. The calculations were performed on the Origin 2000 at the Technical University Dresden.

## A Effective Hamiltonian for the charge–ordered state

This appendix provides the derivation of the effective exchange Hamiltonian for the quarter–filled ladder in the strongly charge–ordered regime, $`(V,U)t`$. The charge degrees of freedom are projected out, i.e., the effective Hamiltonian $`_{\mathrm{eff}}`$ acts in a Hilbert space where every second site is singly occupied. This is analogous to the derivation of the Heisenberg model as large–$`U`$ limit of the half–filled one–band Hubbard model. The present problem maps onto a frustrated $`J_1`$$`J_2`$ spin chain. The effective exchange arises from fourth–order hopping processes which makes the problem more complicated than the half–filled Hubbard model for which the lowest non–trivial contributions arise at second order in $`t`$.

We apply a recently developed cumulant method  to derive $`_{\mathrm{eff}}`$. It is useful to split $`=_0+_1`$ where $`_0`$ contains the dominating interaction terms and $`_1`$ the perturbation caused by hopping. We start with broken translational symmetry from the outset and define a projection operator $`𝒫`$ which projects onto the low–energy space where the charges show perfect checkerboard charge order, i.e., $`n_i=[1+\mathrm{exp}(i\mathrm{𝐐𝐑}_i)]/2`$ with $`𝐐=(\pi ,\pi )`$. This order defines two sublattices which we will denote as $`A`$ and $`B`$ for occupied and unoccupied, respectively. Transitions between states within the $`𝒫`$ space are only possible with four or more hopping processes; the fourth order processes only involve intermediate states outside the $`𝒫`$ space. The fourth–order Hamiltonian can be obtained by fourth–order perturbation theory and is given by

$`_{\mathrm{eff}}`$ $`=`$ $`𝒫_1𝒬{\displaystyle \frac{1}{_0}}𝒬_1𝒬{\displaystyle \frac{1}{_0}}𝒬_1𝒬{\displaystyle \frac{1}{_0}}𝒬_1𝒫,`$ (A1)
where $`𝒬1𝒫`$.

This expression can be easily used to identify the processes contributing to the effective exchange. We first discuss the fourth–order processes leading to $`J_2`$, i.e., the processes coupling two spins located on the same leg of the ladder. They involve exactly one site in between the two originally occupied sites. Therefore, any fourth order hopping process must involve a temporary double occupancy. A suitable classification of the possible processes is shown in Fig. 13.

Any process involves three transition states giving rise to the energy denominator of the final expression for $`J_2`$. It is easy to see that (a) and (b) have transition states with energies $`V_{}+V_{}`$, $`U`$, and $`V_{}+V_{}`$ (the $`V_{}`$ arises from the nearby occupied site on the second leg) whereas (c)–(f) have transition state energies $`V_{}+V_{}`$, $`U+2V_{}`$, and $`V_{}+V_{}`$. The sum of all processes has the form $`J_2(n_{i,\sigma }n_{j,\sigma }+c_{i,\sigma }^{}c_{i,\sigma }c_{j,\sigma }c_{j,\sigma }^{})`$. Examination of the signs shows that the resulting exchange is antiferromagnetic, $`J_2>0`$. For the $`J_1`$ processes there are two intermediate empty sites in between the two occupied sites Therefore, processes involving one or both of these two sites are possible. Particularly interesting are the circular hopping processes contributing to the diagonal exchange $`J_1`$ in which both intermediate sites are involved and no temporary double occupancy occurs. These processes include four sites and the two electrons under consideration are never on the same site. Nevertheless, these processes lead to an effective spin–spin interaction $`𝐒_1𝐒_2`$ (and not only a constant energy shift): For parallel spins the process reproduces the original state, whereas antiparallel spins are always exchanged. This leads to the exchange form $`S_1^{}S_2^++S_1^+S_2^{}+2(S_1^zS_2^z+1/2)`$.

For isotropic hopping and interaction, these considerations can be summarized to $`J_1=2J_2+J_{\mathrm{circ}}`$, where $`J_{\mathrm{circ}}`$ arises from the circular hopping processes; in the general anisotropic case the $`J_1`$ processes will have energy denominators different from the ones quoted above. Collecting all terms leads to the following (general) result for the exchange couplings:

$`J_1`$ $`=`$ $`4t_{}^2t_{}^2\{({\displaystyle \frac{1}{U}}+{\displaystyle \frac{1}{U+2V_{}}})({\displaystyle \frac{1}{(2V_{})^2}}+{\displaystyle \frac{1}{V_{}(V_{}+V_{})}}+`$ (A3)
$`+{\displaystyle \frac{1}{(V_{}+V_{})^2}})+{\displaystyle \frac{1}{V_{}^2(V_{}+V_{})}}\}`$
$`J_2`$ $`=`$ $`{\displaystyle \frac{4t_{}^4}{(V_{}+V_{})^2}}\left\{{\displaystyle \frac{1}{U}}+{\displaystyle \frac{2}{U+2V_{}}}\right\}.`$ (A4)
For isotropic repulsion, $`V_{}=V_{}`$, these expressions reduce to Eq. (10) quoted in the body of the paper.