text:
# Towards the automated reduction and calibration of SCUBA data from the James Clerk Maxwell Telescope

## 1 INTRODUCTION

The Submillimetre Common User Bolometer Array (SCUBA) (Holland et al., 1999) consists of two arrays of bolometers (or pixels); the Long Wave (LW) array has 37 pixels operating in the 750- and 850-$`\mu `$m atmospheric transmission windows, while the Short Wave (SW) array has 91 pixels for observations at 350- and 450-$`\mu `$m. Each of the pixels has diffraction-limited resolution on the telescope (approximately 14.5 and 8 arcsec FWHM respectively), and are arranged in a close-packed hexagon. Both arrays have approximately the same field-of-view on the sky (diameter of 2.3 arcmin), and can be used simultaneously by means of a dichroic.

SCUBA was commissioned in late 1996 and there now exists a large searchable archive of observations111The public archive is hosted by the Canadian Astronomy Data Center. As of February 2001 the archive contained approximately 80,000 separate observations. These data provide an excellent historical record of the telescope and instrument behaviour as well as being an archive of astronomical data. The major challenge of an archive of this size is the data processing. In order to extract useable historical data from the archive it is not feasible to reduce the data manually and an automated approach must be used.

There are three reasons why automated use of the archive is required. Firstly, it provides a unique resource for the investigation of instrument and telescope parameters, leading ultimately to an optimisation of observing practices, secondly the proposed Virtual Observatories will require reduced and calibrated data (e.g. Jenness et al., 2002), and thirdly it allows monitoring of variable sources, especially extragalactic radio sources that are observed frequently for pointing purposes. In this latter respect the present paper forms a companion to our recent work on this topic (Robson, Stevens & Jenness, 2001).

In order to determine its accuracy, any automated technique under consideration must be used to process a standard calibration source of known properties. We present our choice of such a reduction method below, using Uranus as the primary calibrator and then apply the technique to Mars and several compact secondary calibration sources as a consistency check.

## 2 OBSERVATIONS

In order to test any automated reduction system it is important that a suitable dataset is extracted from the archive. This dataset should contain frequent observations of a bright, compact source of known flux taken with a standard observing mode. Although Mars is the primary calibrator at submillimetre wavelengths, Uranus is a better choice for our approach because of its compactness relative to the JCMT beam (semi-diameter approximately 1.7 arcsec over the period in question). The best sampled, homogeneous dataset comes from jiggle-map mode (e.g., Holland et al., 1999) observations of Uranus that are performed at 850/450-$`\mu `$m for calibration purposes, and at 850-$`\mu `$m for pointing checks. We have searched the archive for all jiggle-map observations of Uranus; between April 1997 (the date when the SCUBA data headers were finalised) and May 2001 there were approximately 1500 850-$`\mu `$m observations that were performed in conditions suitable for accurate extinction correction, with 300 of these including data at 450-$`\mu `$m. These data use chop throws of between 60 (all the pointing observations) and 120 arcsec (most of the maps). In October 1999 both the 450- and 850-$`\mu `$m filters were updated to new wideband versions (denoted with a ‘W’ suffix cf. ‘N’ for the previous filters). We include analysis for both filter combinations.

## 3 DATA PROCESSING TECHNIQUES

All SCUBA data are reduced using the SURF (SCUBA User Reduction Facility) package (Jenness & Lightfoot, 1998, 2000) in conjunction with the Starlink software environment (Bly, 2000). To allow us to process the many thousands of observations described in this paper we automated the reduction using the ORAC data reduction pipeline (ORAC-DR) (Economou et al., 1999; Jenness & Economou, 1999). The general reduction method for SCUBA data is documented elsewhere (e.g., Sandell et al., 2001). In this section we comment briefly on the steps that we have tailored to our automated procedure, or which require some discussion in this context.

### 3.1 Extinction correction

Once the data have been demodulated, corrected for beam switching and flatfielded a critical calibration step is to correct for the atmostpheric attenuation. The difficulties of calculating the atmospheric opacity at submillimetre wavelengths are well documented (e.g., Stevens & Robson, 1994) but recent progress has been presented by Archibald et al. (2002) who use skydips taken at both the JCMT and the Caltech Submillimeter Observatory (CSO) to provide a reliable record of the opacity spanning a period of several years. We refer the reader to Archibald et al. (2002) for a detailed discussion but give a brief outline below.

Automated reduction relies critically on the quality of the opacity measurement – simply applying the nearest skydip or CSO tau measurement to the data is dangerous for a number of reasons. Specifically, SCUBA skydips are usually not taken frequently enough to be useful in variable conditions. In addition, although the CSO tau measurements are taken much more frequently (every 10 mins), the instrumental noise on the CSO tau meter is quite high, it produces erroneously large readings occasionally and operates at a fixed azimuth.

With these points in mind we proceeded as follows: a polynomial was fitted to the re-scaled (Archibald et al., 2002) CSO tau data in regions where the SCUBA skydips were in reasonable agreement; times when the CSO tau meter either did not agree with the SCUBA value or when the signal was varying rapidly were not included in these fits. This method ensured that the atmosphere seen by the tau meter was the same as that seen by SCUBA. In practice, the pipeline ignores data taken in periods when no fitted extinction value is available. Whilst this means that no data are reduced when, for example, the CSO tau meter is broken it is a necessary requirement of the automation process. The error of the 225 GHz opacity from the fit is usually of order 0.005; if this error is assumed to be the dominant error in the extinction then for a source at airmass 1.5 this contributes to the final result an error of approximately 3 per cent at 850-$`\mu `$m and 20 per cent at 450-$`\mu `$m.

### 3.2 Despiking

Occasionally, the time-series data contains spikes (for example from cosmic rays) and these must be removed. In the simplest case this involves analysing each bolometer in turn, removing data that are significantly different to the mean. The complication for SCUBA mapping is that each bolometer sees different parts of the target during the observation (“jiggling” on and off a source) such that a simple mean no longer provides valid statistics. For long observations it is possible to compare data points that are from similar points on the sky (effectively equivalent to subtracting source signal from the data prior to despiking) but for short observations, such as those from bright calibrators described below, there are simply too few data points for a reliable detection of spikes when near the source in question. For this work, we have compromised by performing the despiking for each bolometer independently but with a large enough clipping level (5 sigma) that we can detect the largest spikes without fear of clipping real source structure. This approach does mean that we are still susceptible to small spikes in the data but the overall effect on the integrated sum will be minimal when averaged over all the observations in this sample.

### 3.3 Sky-noise removal

Nodding and chopping do not remove all the atmospheric emission noise from the data so an additional step is required that takes into account the fact that the source signal is constant in time and the sky-noise is varying in time but coherent across the array. Removal of sky-noise from jiggle maps requires knowledge of regions of the map that are uncontaminated by source flux. Our criterion of selecting only compact sources for analysis simplifies this procedure enormously. For pointing observations it is safe to use the median of all the bolometers to determine the sky-noise signal; the pointings use an azimuthal 60 arcsec chop-throw so that each map contains the source at the centre of the map and two negative images with half the flux. The latter images are present because the chop throw is smaller than the field-of-view with the result that the negative beams of some of the ‘sky’ bolometers fall onto the central source. For standard mapping observations the chop throw is generally chosen to be larger than the field-of-view in order to avoid this effect. It was, however, found that median sky removal was equally well applicable to these observations. See Jenness, Lightfoot & Holland (1998) and Archibald et al., (2002) for examples of the sky-noise removal algorithm.

### 3.4 Residual sky removal

The sky-noise removal algorithm currently assumes that the sky signal adds a d.c. offset across the whole array. In some extreme cases this is not true and a linear gradient is visible across the array. To compensate for this a plane is fitted to the regridded image and subtracted before the data are processed. For the data considered here this procedure has only a minor effect because the source is in the centre of the image but is performed for completeness (when the sky-noise is a simple offset this step will have no effect).

### 3.5 Signal determination

The penultimate step in the reduction process is to determine the signal from the source. In the present case there is always a single point-like source at the centre of the map and this simplifies the analysis.

Having first made sure that we have accurately located the centroid position there are two ways to determine the signal from a source: (1) measure the peak flux, and (2) measure the flux in an aperture. Both of these approaches are valid, we calculated both, and compare them in Section 4.1.

The first step for case (2) is to select an aperture of suitable diameter. For the special case of pointing observations the aperture diameter can be no larger than 60 arcsec otherwise contamination from the negative images of the central source becomes a problem. In general, there are two competing considerations, the effects of which have to to traded-off to obtain the optimum aperture size for the particular problem. On the one hand, it is desirable to use as large an aperture as possible so that most of the time dependent error beam pattern is enclosed for both source and calibrator. This is especially important at 450-$`\mu `$m where the error lobe contribution typically exceeds 100 per cent for chop throws larger than about 2 arcmin. On the other hand, increasing the aperture also increases the measured noise. The final trade-off is thus between accuracy of calibration and signal-to-noise. Experimenting with various apertures we found that 40 arcsec is a good choice for point-like sources. Increasing the aperture to 60 arcsec results in an increase in signal of only 5–6 per cent at 850-$`\mu `$m under most conditions considered, and correspondingly 13 per cent at 450-$`\mu `$m whereas the signal-to-noise degrades by a factor of 2/3 (the error on the integrated sum is proportional to the radius of the aperture). Apertures of significantly less than 40 arcsec produced unacceptable scatter in the signals measured for Uranus with the standard deviation increasing by almost a factor of 2 (see Table 1).

The signal-to-noise is determined by looking at small regions towards the edge of the map and determining the standard deviation of the pixels in the aperture. Bolometers with uncharacteristically high noise are often apparent when maps are made of faint sources. In order to avoid these and other effects like spikes, which are not always easy to remove from small datasets, a number of apertures are used and the one with the smallest standard deviation is selected to determine the overall error in the integrated sum (the aperture measures the noise-per-pixel). We verified that this procedure is valid by reducing a number of maps by hand.

Once the data have been reduced automatically there are inevitably a small number of points that must be removed by hand because of problems with the data acquisition, or with the instrument or telescope. The main problems are occasional pointing observations where the source is not centred on the array, where the dish has been heated by the Sun and is grossly out of shape, when SCUBA is not operating correctly, or when the polarimeter is attached. For these cases, the observing logs and final data products are inspected and the point is removed if warranted by the circumstances.

## 4 Results

### 4.1 Flux conversion factor

The final step in the reduction process is the application of a suitable ‘flux conversion factor’ (FCF) to the reduced data so that the measured signal in volts is converted into physically useful units such as Janskys. In order to automate the reduction of arbitrary data we need to know the FCF history of the telescope and instrument. We do this by measuring the response of our primary calibrators (Mars and Uranus) throughout the period of interest and comparing this with the expected flux derived from a standard model for planetary fluxes (Privett, Jenness & Matthews, 1998; Griffin & Orton, 1993; Wright, 1976). The FCF can take two forms: as mentioned in Section 3.5 we can either calibrate on the ‘peak’ of a point source giving the traditional radio astronomy units of Jy beam-1 ($`FCF_{peak}=S_{peak}/V_{peak}`$ where $`S_{peak}`$ is the flux in the beam and $`V_{peak}`$ is the measured peak voltage), or we can ‘integrate’ in an aperture giving units of Jy arcsec-2 ($`FCF_{int}=S_{tot}/[V_{int}A]`$ where $`S_{tot}`$ is the total flux of the calibrator, $`V_{int}`$ is the integrated sum in volts and $`A`$ is the pixel area in arcsec2). In this section we discuss the merits of the two schemes by applying the automated pipeline reduction to Uranus.

The results of our Uranus reduction for the 850N and 450N filters are shown in Figs. 1 and 2 respectively. Similar trends are seen for the 450W/850W filter pair, and the FCF values are summarized in Table 1. The top left panel of Fig. 1 shows the FCF calculated by summing the signal in a 40 arcsec aperture plotted against time for a period of $`23`$ yrs. In general the FCF is remarkably stable given the variable weather/telescope/instrument parameters over this relatively long time period. There is evidence for the FCF to be slightly higher during 1997 and early 1998 than in subsequent times. It is presently unclear as to why this should be the case (if real) but we note that SCUBA was warmed up in early 1998 because of varying base temperature problems and a similar, small, change in FCF was seen following a warm up in the middle of 2000. Also apparent from this panel is that a small fraction of the FCF values are discrepant by a large amount ($`>50`$ per cent). The effect is shown clearly in the top right panel of Fig. 1 where we show a histogram of the integrated FCF values. The dotted line shows a Gaussian fit to the data giving a standard deviation of around 10 per cent as used by Robson et al. (2001). The solid line shows the result of fitting a double Gaussian to the data (Table 1 lists the average for the filter over the entire period) which gives an improved fit. Note the tail of high FCF values.

The lower left panel of Fig. 1 shows the ‘integrated’ FCF versus ‘peak’ FCF. Variations in the sensitivity of SCUBA connected with warm-ups must contribute significantly to this general trend shown in this plot. However, on any given night we believe that the linear trend, and hence the dominant source of uncertainty in the calibration of SCUBA data is produced by the thermal state of the antenna; it is well known that the focus of the antenna changes rapidly around sunset and sunrise reflecting its changing shape as it thermally relaxes. During these periods the observed beam pattern is also variable because power is removed from the main beam and spread out into the error lobes. The fact that we see a trend demonstrates that the error beam spreads out beyond the $`40`$ arcsec aperture, but the key point to note is that whereas the ‘peak’ FCF varies by about a factor of two, the ‘integrated’ FCF varies by $`<50`$ per cent. Furthermore, the size of the chop throw can vary the ‘peak’ FCF by as much as $`510`$ per cent (the beam smears out a small amount in the direction of the chop throw and the sensitivity is reduced) but affects the ‘integrated’ FCF much less. The best reduction method for our automated approach is thus to use the ‘integrated’ FCF for flux determination. A corollary is that photometry mode observations, which are made with a single pixel and a $`2`$ arcsec, $`3\times `$3 square jiggle pattern, and hence are calibrated in an analogous manner to the FCF ‘peak’ method described above, can at times be susceptible to large calibration uncertainties. A quantitative assessment would require detailed analysis of calibration data from individual nights which is beyond the scope of the automated approach described in this paper, but if high values of the FCF are measured then it may well be good observing practice to use SCUBA in ‘jiggle-map’ mode rather than photometry mode. Although mapping observations are intrinsically less efficient than photometry observations, the ability to measure the signal in a suitably sized aperture will allow observations to be calibrated accurately during the period when the antenna is changing shape.

The extreme values plotted in the lower left panel of Fig. 2 probably occur when the antenna is grossly out of shape which may reflect the nature of the observing, for example if observing was carried out in the daytime or if the dish was pointed towards the Sun in early evening. These data points are calculated from maps in which the beam shape is so obviously poor that they are often removed by hand in the post-pipeline phase of the data reduction. Evidence that this is indeed the case is provided in the bottom right panel of Fig. 1 where we plot ‘integrated’ FCF versus time of day Here we see that the majority of high FCF values ($`>1.3`$) occur during the first half of the night or after sunrise. Note also that the FCF can remain high until the middle of the night in some cases but is almost always lower than about 1.2 during stable night-time conditions after midnight. This may indicate that once the antenna has been heated it can take several hours to recover its optimum shape. Similar trends are seen at 450-$`\mu `$m (Fig. 2) where a single Gaussian provides a good fit the ‘integrated’ FCF values giving a standard deviation of about 20 per cent. Experience shows that the uncertainties derived from our automated approach are similar to those applied by observers when processing their data ‘by hand’, demonstrating the validity of the method for point-like sources.

The ratio of ‘peak’ to ‘integrated’ FCF provides some idea of the magnitude of the error beam within the aperture at these wavelengths. At 850-$`\mu `$m the ratio for a $`40`$ arcsec aperture is approximately 250 arcsec2 beam-1, an area equivalent to a $`15`$ arcsec Gaussian beam (the measured value is approximately 14.5 arcsec). At 450-$`\mu `$m the ratio is approximately $`120`$ arcsec2 beam-1, equivalent to a $`10`$ arcsec beam (the main beam is approximately $`7.5`$ arcsec). The average error lobe contribution within a $`40`$ arcsec aperture is thus about $`4`$ per cent at 850-$`\mu `$m and $`40`$ per cent at 450-$`\mu `$m.

### 4.2 Mars

In order to test the output from the automated reduction, including a test of the validity of our calibration, we tested the complete system with observations of Mars. An aperture of 60 arcsec was used to maximize signal (signal-to-noise is not an issue) and to guarantee the inclusion of the majority of the flux (Mars had a maximum diameter of 15 arcsec during the period). The results are shown in Fig. 3 along with the expected flux provided by the fluxes software (Privett et al., 1998) which uses the model for mean surface temperature derived by Wright (1976) and knowledge of the SCUBA filters. The excellent agreement between the observed and predicted fluxes provides convincing support for our automated reduction procedure and the Mars model.

### 4.3 Secondary calibrators

Now that the system has been shown to generate reasonable results we can apply these techniques to archive data. Here we present the results from processing of the most common JCMT secondary calibrators (see Table 2); the results from extragalactic radio pointing sources are presented in Robson et al. (2001). Table 3 shows the results of the automated reduction of the selected calibrator sources using a 40 arcsec aperture. These data are selected in the same way as for Uranus with only pointing and mapping observations retrieved from the data archive. This list does not include potential candidates, such as VY CMa and TW Hya, mentioned in Sandell (1994) simply because the archive does not currently contain enough observations to justify reliable automated reduction.222By the end of the year 2000 there were only 12 suitable observations of VY CMa and 20 suitable observations of TW Hya in the archive. Of those only 2 of the VY Cma and 6 of the TW Hya observations had 450-$`\mu `$m data. At 850-$`\mu `$m the results from automated reduction were 2.1$`\pm `$0.1 Jy for VY Cma and 1.46$`\pm `$0.2 Jy for TW Hya Since all these sources are used for pointing there are significantly more observations available at 850-$`\mu `$m than at 450-$`\mu `$m. IRAS 16293$``$2422 has not been included in this analysis because of the extended nature of the source (e.g., Sandell, 2001).

Table 3 indicates that the automated reduction shows excellent agreement with Sandell (2001), especially at 850-$`\mu `$m. The fundamental difference is that the numbers presented here are integrated over an aperture and will therefore also include flux from extended emission whereas Sandell (2001) calibrated his data in Jy/beam.

#### 4.3.1 HL Tau

The flux density shows excellent agreement with Sandell (2001) at both wavelengths indicating that there is no evidence for extended structure out to $`20`$ arcsec radius when performing chopped observations. Fig. 4 shows the light-curve. Whilst the plot raises the possibility of source variation (possibly with a maximum at 1998.5, amplitude of 100 mJy and period of $`600`$ d) a chi-square fit indicates that the data are consistent with zero variability.

#### 4.3.2 CRL 618

CRL 618 is the most popular secondary calibrator in use at the JCMT. The flux density shows excellent agreeement with Sandell (2001) and Sandell (1994) at both wavelengths (the latter taken using UKT14 (Duncan et al., 1990)) indicating that there is no evidence for extended structure out to $`20`$ arcsec radius when performing chopped observations. Knapp, Sandell & Robson (1993) indicate a possibility of source variability in the millimetre region but Fig. 5 shows that no obvious variability can be seen at 850-$`\mu `$m between mid-1997 and mid-2001. Apparent short term trends (of the order of $`0.2`$ yr) in the data stream are most likely connected with the observing/reduction methodology or with FCF variations due to dish heating rather than to source variability.

#### 4.3.3 OH 231.8+4.2

We get excellent agreement with the results of Sandell (2001). OH231.8 is a long period variable star and the time series show clear evidence for variability at 850-$`\mu `$m. A least-squares fit of a sine-wave to the time series indicates a period of $`630\pm 20`$ d. This is somewhat lower, but not necessarily inconsistent with, that seen in the near-infrared: $`650`$ d (Feast et al., 1983), $`684\pm 40`$ (Bowers & Morris, 1984) and $`708\pm 6`$ (Kastner et al., 1992), but our data are relatively noisy. Indeed, Kastner et al. (1992) suggest that the period may have changed slightly over the years and it is therefore possible that the period has shortened recently.

Fig. 6 includes a sine-wave with a period of $`630`$ d of the form:

$$S(Jy)=2.75+0.4\mathrm{sin}\left(2\pi \frac{Y1999.08}{1.75}\right)$$
(1)
where $`Y`$ is the year, or alternatively:

$$S(Jy)=2.75+0.4\mathrm{sin}\left(2\pi \frac{JD2451207}{630}\right)$$
(2)
where $`JD`$ is the Julian day. If the sine wave is subtracted from the time series data the error in the residual is reduced to $`0.2`$ Jy. The uncertainty in the period makes it difficult to extrapolate the light curve from earlier studies with any accuracy. Therefore it is not possible to reliably determine the lag time between the light maximum in the near-infrared and that of the submillimetre. The 450-$`\mu `$m data are too noisy to show source variability.

#### 4.3.4 IRC+10216 (CW Leo)

This is a variable star with a period of $`635`$ d as seen in the optical (Alksnis, 1989) and at 1.1-mm (Sandell, 1994). When compared to Sandell (2001) and Sandell (1994) this source is clearly extended in our aperture. Using a source size of $`8.9`$ arcsec at 850-$`\mu `$m and $`4`$ arcsec at 450-$`\mu `$(Sandell, 2000) the results from Sandell (2001) would be equivalent to 8.4- and 16.4-Jy in excellent agreement with our results.

Fig. 7 shows the measured light-curve at 850-$`\mu `$m along with a sine-wave of period $`635`$ d of the following functional form:

$$S(Jy)=8.8+0.95\mathrm{sin}\left(2\pi \frac{Y1999.125}{1.74}\right)$$
(3)
where $`Y`$ is the year, or alternatively:

$$S(Jy)=8.8+0.95\mathrm{sin}\left(2\pi \frac{JD2451224}{635}\right)$$
(4)
where $`JD`$ is the Julian day. If the sine-wave is subtracted from the time-series data the error in the residual is reduced to $`0.56`$ Jy. These data indicate that the source is at a maximum approximately $`120`$ d after the maximum expected in the optical from Alksnis (1989) and $`72`$ d ($`0.2`$ yr) after the maximum expected from the 1.1-mm observations of Sandell (1994). It is likely that the discrepancy in phase between the 850-$`\mu `$m and 1.1-mm data is not significant given the difficulty in obtaining the 1.1-mm data with a single pixel instrument. The 450-$`\mu `$m data are too noisy and do not show reliable evidence of variability.

#### 4.3.5 CRL 2688

CRL 2688 has a $`25`$ per cent higher 450-$`\mu `$m flux here than presented in Sandell (2001) suggesting that the source is extended at this wavelength. A map was generated combining the data from a number of stable nights from the archive to test this, taking particular care that any pointing shifts are accounted for during the reduction. The measured ‘source size’ of $`9.5`$ arcsec is larger than that measured from CRL 618 (known to be unresolved), using the same techniques, of $`8.0`$ arcsec (see also, Coulson, 2001). This corresponds to a deconvolved source size of approximately $`5`$ arcsec and would account completely for the lower flux measured by Sandell (2001). With a source size of 5 arcsec we would expect a much larger discrepancy at 450-$`\mu `$m than at 850-$`\mu `$m which is indeed the case. The time-series in Fig. 8 indicates that the flux density has been constant since the start of 1998.

## 5 CONCLUSIONS

This paper presents a means of automating the reduction and calibration of 850- and 450-$`\mu `$m SCUBA data. The limitations of the technique are investigated by analysing archival data for the standard submillimetre calibration source, Uranus. We show that our method gives excellent agreement between observed and predicted submillimetre fluxes for Mars. It provides a useful tool for investigating trends in instrument and telescope parameters – these studies are essentially made practical for the first time with our automated procedure, and will ultimately lead to a better understanding of the uncertainties involved in the calibration of data from the JCMT and other submillimetre observatories.

We use our automated procedure to calculate light-curves and integrated fluxes for selected, compact JCMT secondary calibrators. The periodic variability of IRC+10216 is confirmed and there is strong evidence to support similar behaviour for OH231.8 at 850-$`\mu `$m; the first time this variability has been confirmed in the submillimetre. The remaining calibrators have been constant during the period 1997.5 to 2001.5.

## ACKNOWLEDGMENTS

The James Clerk Maxwell Telescope is operated by the Joint Astronomy Centre in Hilo, Hawaii on behalf of the parent organizations PPARC in the United Kingdom, the National Research Council of Canada and The Netherlands Organization for Scientific Research. J.A.S. acknowledges the support of a PPARC PDRF. We acknowledge the support software provided by the Starlink Project which is run by CCLRC on behalf of PPARC. Thanks to Wayne Holland, Remo Tilanus, Iain Coulson and Göran Sandell for useful discussions on this topic.


format: text

text:
# Superconducting charge-ordered states in cuprates

## I Introduction

A series of recent experiments have highlighted the importance of spin and charge ordering tendencies in the cuprate superconductors. Static stripe order has been established in Nd-doped La2-xSrxCuO4 which appears to co-exist with superconductivity at very low temperatures. Recently, superconducting samples of underdoped YBa2Cu3O6.35 and nearly optimally doped Bi2Sr2CaCu2O8+δ with $`T_c`$ up to 90 K have been found to show signatures of charge ordering. In addition, a variety of compounds display strong dynamic spin and charge fluctuations, which are enhanced by the application of moderate magnetic fields . This emphasizes that even materials, where no static order can be observed in the absence of an external field, are in close proximity to a critical point where spin and/or charge order is established in the superconducting state. Theories based on this assumption  successfully explain a number of NMR  and neutron scattering  experiments.

In this paper, we will focus on superconducting states with static charge order, but dynamic spin fluctuations – this appears to be realized in the experiments of Ref. . Such states are found upon doping paramagnetic Mott insulators on the square lattice . The undoped quantum paramagnet has broken translational symmetry associated with spontaneous bond charge (or spin-Peierls) order; at small carrier concentration, $`\delta `$, this order persists and co-exists with anisotropic superconductivity; a $`d`$-wave superconductor with full square lattice symmetry appears above a critical $`\delta `$. \[“Charge order” is defined very generally as spatial modulation in any SU(2)-invariant observables, such as local density of states (LDOS) per site, or kinetic exchange energy per lattice bond; the modulation in the total site charge density can be small due to long-range Coulomb interactions.\] A particular feature of the superconducting charge-ordered states found in Refs.  is that the modulation is bond-centered and its real-space period $`p`$ always takes even integer values (in units of the Cu lattice spacing); this is in contrast to the continuous doping evolution of the ordering wavevector, $`1/p`$, usually assumed in the so-called “Yamada plot” . Interestingly, the experimental results of Refs.  appear to be remarkably well described by the states proposed in Refs. , as has also been discussed in recent work  which appeared while this paper was being completed: Neutron scattering on underdoped YBa2Cu3O6.35 shows charge order with a real-space period of 8 lattice sites, whereas Bi2Sr2CaCu2O8+δ close to optimal doping displays $`p=4`$ modulation ; furthermore the STM data indicate an even-period modulation not only in the site charge density, but also in the bond kinetic energy and perhaps the bond pairing amplitude .

We note that numerical studies of the $`tJ`$ model  have observed bulk charge order with period $`p=4`$ at doping level 1/8; and paired hole states with different types of charge order in both insulators and superconductors have been discussed elsewhere .

In the past, most theoretical work has been focussed on states with one-dimensional (1d) charge modulation, often referred to as stripes. However, recent STM experiments  indicate LDOS modulations in both $`x`$ and $`y`$ directions in a single CuO2 plane (although a small anisotropy is observed). Theoretically, charge density wave (CDW) fluctuations are expected in both directions  on the disordered side of the charge ordering transition (if the system has no intrinsic lattice anisotropy). Moving to the ordered side, it depends on microscopics whether CDW order in one or in two directions condenses, leading to stripe-like or two-dimensional (2d) modulations, respectively.

The purpose of this paper is twofold: In Sec. II, we re-examine the mean-field theory of Refs. , to investigate the possible existence of and the mechanism leading to superconducting charge-ordered states with 2d charge modulation. In Sec. III, we turn to a detailed discussion of the charge modulation pattern, by calculating the energy dependence of the LDOS Fourier component at the ordering wavevector and comparing it with the measurements of Ref. .

## II Mean-field theory

We start with re-analyzing the large-$`N`$ theory of Refs. , which provides a microscopic description of doping mobile charge carriers into a paramagnetic Mott insulator . We consider an extended $`tJ`$ Hamiltonian for fermions, $`c_{i\alpha }`$, on the sites, $`i`$, of a square lattice with spin $`\alpha =1\mathrm{}2N`$ ($`N=1`$ is the physical value):

$`_{tJV}={\displaystyle \underset{i>j}{}}`$ $`[{\displaystyle \frac{t_{ij}}{N}}c_{i\alpha }^{}c_{j\alpha }+\mathrm{H}.\mathrm{c}.+{\displaystyle \frac{V_{ij}}{N}}n_in_j`$ (2)
$`+{\displaystyle \frac{J_{ij}}{N}}(𝐒_i𝐒_j{\displaystyle \frac{n_in_j}{4N}})].`$
Here $`n_i=c_{i\alpha }^{}c_{i\alpha }`$ is the on-site charge density, and the spin operators $`𝐒_i`$ are fermion bilinears times the traceless generators of $`\mathrm{Sp}(2N)`$. For most of the following, the fermion hopping, $`t_{ij}`$, and exchange, $`J_{ij}`$, will be restricted to nearest-neighbor terms, $`t`$ and $`J`$; for the detailed comparison with the Ref.  we will introduce 2nd neighbor hopping, $`t^{}`$. The electronic Coulomb interaction is represented by the on-site constraint $`n_iN`$ and the off-site repulsive interactions $`V_{ij}=V/|r_{ij}|`$. The average doping $`\delta `$ is fixed by $`_in_i=NN_s(1\delta )`$, where $`N_s`$ is the number of lattice sites. To proceed, we represent the spins by auxiliary fermions $`f_{i\alpha }`$ and the holes by spinless bosons $`b_i`$, such that the physical electrons $`c_{i\alpha }=b_i^{}f_{i\alpha }`$, and the necessary Hilbert space constraint is implemented by Lagrange multipliers $`\lambda _i`$. Via a Hubbard-Stratonovich decoupling of the antiferromagnetic interaction we introduce link fields $`Q_{ij}`$, defined on the bonds of the square lattice. After taking the limit $`N\mathrm{}`$, the slave bosons $`b_i`$ condense, $`b_i=\sqrt{N}b_i`$, the $`Q_{ij}`$ and $`\lambda _i`$ take static saddle-point values, and we are left with a bilinear Hamiltonian which can be diagonalized by a Bogoliubov transformation. At the saddle point, the slave boson amplitudes fulfill $`_ib_i^2=N_s\delta `$, and the link fields are given by $`NQ_{ij}=𝒥^{\alpha \beta }f_{i\alpha }^{}f_{j\beta }^{}`$, where $`𝒥^{\alpha \beta }`$ is the antisymmetric $`\mathrm{Sp}(2N)`$ tensor; for further details see Ref. .

Various ground states obtained from the numerical solution of the above mean-field equations have been discussed in Refs. . At $`\delta =0`$ the ground state is a fully dimerized, insulating spin-Peierls state, i.e., it has period-2 bond charge order. At low doping $`\delta `$, the bare large-$`N`$ $`tJ`$ model tends to phase separation, and the inclusion of moderate Coulomb interaction, $`V`$, leads to the formation of bond-centered, superconducting stripes with a 1d charge modulation. Large doping destroys charge order and leads to a pure $`d`$-wave superconducting ground state, in this regime the large-$`N`$ approach reduces to the usual BCS mean-field theory with renormalized hopping matrix elements.

Motivated by the STM results of Refs.  we have searched for additional saddle-point solutions with 2d charge modulation; we have restricted the attention to states in which the charge distribution respects the 90 degree rotation symmetry of the lattice. Interestingly, there are several such saddle points which were overlooked in Ref. . In most of the low and intermediate doping region the states with 1d and 2d modulation are close in energy; at small doping the 1d stripe-like states are preferred, whereas the 2d checkerboard-likes states are lower in energy in a certain interval of intermediate doping and Coulomb repulsion.

We have therefore concentrated on states with a real-space periodicity $`p=4`$, leading to a unit cell of $`4\times 4`$ sites . (Recall that bond-centered stripes with the spatial period pinned to four sites were found over a rather large range of doping values in Refs. .) The most favorable mean-field states with 2d charge modulation are characterized by holes arranged in intersecting bond-centered stripes, i.e., the hole concentration is large on 12 sites and zero on 4 sites, see Fig. 1. The actual filling in the hole-rich regions is smaller than 1 and depends on microscopic parameters; due to the strong pairing correlations the 2d CDW states are good superconductors (with a $`d`$-wave-like pairing symmetry). Such $`4\times 4`$ plaquette states occur as large-$`N`$ ground states of the Hamiltonian (2) for doping levels between approximately 13% and 25%. A sample phase diagram for doping $`\delta =20\%`$ is shown in Fig. 1 (see also Fig. 3 of Ref. ). For small $`t/J`$ phase separation tendencies do not occur, therefore the ground state is a doped spin-Peierls state with homogeneous site charge distribution. At intermediate $`t/J`$, frustrated phase separation leads either to stripes or to plaquette states; at larger $`t/J`$ the kinetic energy becomes dominant and weakens the phase separation tendencies, consequently a homogeneous $`d`$-wave state is reached. At very large Coulomb repulsion, the holes arrange into an insulating Wigner crystal.

The occurence of plaquette-modulated CDW states in favor of stripes can be understood as an interplay of exchange, kinetic, and Coulomb energy as follows: The exchange term prefers dimerization between spins on neighboring sites and tends to expel holes – this produces the overall dimer structure of the $`Q_{ij}`$, and leads to a fraction of sites being undoped. The kinetic energy is lowered by possible hopping processes, i.e., by neighboring sites with non-zero hole density. Clearly, the exchange term prefers stripes, whereas the kinetic term prefers plaquettes, where hopping in two directions is possible. Now, the Coulomb energy of the 2d modulated state is significantly lower than that of a stripe state, because the charge inhomogeneity is smaller in the plaquette state (which is also closer to a crystalline arrangement of charges). From this discussion it is clear that at low doping, where the physics is dominated by the exchange term, stripe states are preferred. With increasing doping the kinetic energy becomes more and more important, which leads to 2d modulated CDW states at moderate values of $`V`$, before a homogeneous $`d`$-wave superconductor becomes the ground state.

In all CDW states superconductivity competes with charge order. In the stripe states superconductivity is very weak due to the strong anisotropy: bulk superconductivity is established only by Cooper pair tunneling between the stripes. In contrast, the plaquette CDW states are much better superconductors, due to the full 2d character of the charge distribution. This trend is consistent with the low $`T_c`$ in Nd-doped La2-xSrxCuO4 compared to apparently charge-ordered YBa2Cu3O6.35 and Bi2Sr2CaCu2O8+δsamples .

## III Modulation in the local density of states

After having established the possible occurence of plaquette CDW states in the large-$`N`$ theory for the $`tJ`$ model, we turn to a detailed analysis of the corresponding STM signal. Hoffman et al.  have introduced a STM technique of atomically resolved spectroscopic mapping, which allowed to detect LDOS modulations around vortex cores with real-space period 4, i.e., at wavevectors $`𝐊_x=(\pi /2,0)`$ and $`𝐊_y=(0,\pi /2)`$. Howald et al.  used this technique to map the energy dependence, $`\rho _𝐊(\omega )`$, of the spatial Fourier component of the LDOS at the ordering wavevectors $`𝐊_{x,y}`$. This energy dependence has been recently discussed within a model for the pinning of SDW/CDW fluctuations by inhomegeneities , and in an analysis of different patterns of translational symmetry breaking in $`d`$-wave superconductors .

For a comparison with the experimental situation  we restrict ourselves again to states with period-4 modulation; we will employ hopping parameters $`t`$, $`t^{}`$ that yield a realistic band structure. Furthermore we have to keep in mind the shortcomings of the large-$`N`$ theory: the precise location of the phase boundaries is not reliable, and the theory underestimates fluctuations. Therefore we will work with superconducting gap values close to the experimentally observed ones, and discuss both self-consistent mean-field solutions as well as states where translational symmetry breaking is imposed by hand in the Hamitonian.

To obtain initial information about the possible forms of $`\rho _𝐊(\omega )`$ we start by considering $`d`$-wave superconductors with additional modulation in one of the following quantities: site charge density, bond charge density (kinetic energy), pairing amplitude. Such an analysis has also been independently performed by Podolsky et al. , but here we are interested in 2d modulations and furthermore diagonalize the mean-field Hamiltonian exactly for the $`4\times 4`$ unit cell. The Hamiltonian thus has the form $`_{\mathrm{BCS}}+_{\mathrm{mod}}`$, where

$$_{\mathrm{BCS}}=\underset{𝐤}{}ϵ_𝐤c_{𝐤,\sigma }^{}c_{𝐤,\sigma }+\underset{𝐤}{}\mathrm{\Delta }_𝐤(c_{𝐤,}^{}c_{𝐤,}+h.c.)$$
(3)
in standard notation, with $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }_0(\mathrm{cos}k_x\mathrm{cos}k_y)/2`$; $`_{\mathrm{BCS}}`$ is equivalent to the $`\mathrm{Sp}(2N)`$ mean-field theory presented above in the region where the large-$`N`$ ground state is a pure $`d`$-wave superconductor (with the correspondence $`\mathrm{\Delta }_{ij}=J_{ij}Q_{ij}`$ where $`\mathrm{\Delta }_{ij}`$ is the real-space Fourier transform of the energy gap $`\mathrm{\Delta }_𝐤`$).

The modulation is introduced via $`_{\mathrm{mod}}`$: for site CDW we add $`_{\mathrm{mod}}^{\mathrm{site}}=V_0^{\mathrm{site}}_if(𝐑_i)c_{i\sigma }^{}c_{i\sigma }`$, for bond CDW we have $`_{\mathrm{mod}}^{\mathrm{kin}}=V_0^{\mathrm{kin}}_{ij}f[(𝐑_i+𝐑_j)/2]c_{i\sigma }^{}c_{j\sigma }`$, and a pairing modulation is given by $`_{\mathrm{mod}}^{\mathrm{pair}}=V_0^{\mathrm{pair}}_{ij}f[(𝐑_i+𝐑_j)/2](c_i^{}c_j^{}+h.c.)`$. The function $`f(𝐑)`$ describes modulation strength and pattern, and we will concentrate on 2d bond-centered period-4 modulations with $`f(𝐑)=[\mathrm{cos}(\pi R_x/2+\pi /4)+\mathrm{cos}(\pi R_y/2+\pi /4)]/2`$.

Diagonalization of $`_{\mathrm{BCS}}+_{\mathrm{mod}}`$ yields the local density of states for each site of the unit cell, from which we find $`\rho _𝐊(\omega )`$ by Fourier transformation (the real-space origin is chosen such that $`\rho _𝐊(\omega )`$ is real); note that $`𝐊=(\pi /2,0)`$ and $`(0,\pi /2)`$ are equivalent with the above choice of $`f(𝐑)`$. Results for $`\rho _𝐊(\omega )`$ are shown in Fig. 2 for the three modulation cases listed above. If we compare the curves in Fig. 2 with the STM result in Fig. 3 of Ref. , which shows a peak in the magnitude of $`\rho _𝐊(\omega )`$ at subgap energies $`|\omega |/\mathrm{\Delta }_02/3`$, it is clear that the experiments are not well described by a site charge modulation alone. In contrast, our result for a modulation in the pairing amplitude comes closest to the curves of Ref. . We note that our results in Fig. 2 are somewhat different from the ones of Ref. , this may be due to the 2d character of the modulation considered here and due to the approximations employed in Ref. .

It is clear that the experimentally realized CDW state will have modulations in all quantities invariant under spin rotations and time reversal. This can – at least in part – be captured by a self-consistent solution of the mean-field equations. A natural candidate is given by the plaquette state found as large-$`N`$ ground state above (Fig. 1). Thus, we employ the Hamiltonian $`_{\mathrm{BCS}}+_{\mathrm{mod}}^{\mathrm{site}}`$, where the pair-field $`\mathrm{\Delta }`$ in $`_{\mathrm{BCS}}`$ (3) is determined self-consistently from $`\mathrm{\Delta }_{ij}=J_{ij}𝒥^{\alpha \beta }c_{i\alpha }^{}c_{j\beta }^{}`$, and $`_{\mathrm{mod}}^{\mathrm{site}}`$ imposes a weak modulation of the site charge density as shown in Fig. 3 (the strong modulation found in the large-$`N`$ limit will certainly be weakened by fluctuation corrections beyond the large-$`N`$ theory). Fig. 3 displays a corresponding LDOS modulation $`\rho _𝐊(\omega )`$. The agreement with the available experimental data  is not satisfying, in particular $`\rho _𝐊(\omega )`$ does not show a large peak at energies below the bulk superconducting gap.

This fact and the results in Fig. 2 led us to consider an additional effect not captured in the mean-field calculations: On general symmetry grounds, a static charge modulation will lead to a real-space modulation in the effective pairing interaction, because the CDW influences the local spin fluctuation spectrum. On the mean-field level, this can be phenomenologically accounted for by a modulation in the exchange interaction $`J`$. Therefore, we have studied self-consistent solutions of the mean-field theory, $`_{\mathrm{BCS}}+_{\mathrm{mod}}^{\mathrm{site}}`$, as above, but in addition to a weak site charge modulation in $`_{\mathrm{mod}}^{\mathrm{site}}`$ we imposed a modulation of the $`J_{ij}`$ exchange interaction, $`J_{ij}=J_0+V_0^Jf[(𝐑_i+𝐑_j)/2]`$ for nearest neighbor sites $`i`$ and $`j`$, which leads to corresponding modulations in both the pair fields $`\mathrm{\Delta }_{ij}`$ and the bond charge density (kinetic energy). Results for the LDOS Fourier component $`\rho _𝐊(\omega )`$ are shown in Fig. 4, with a rather good agreement with the experiments of Ref. .

## IV Conclusions

Summarizing, we have studied superconducting charge-ordered states of doped Mott insulators. Within a large-$`N`$ theory we have established that ground states with 2d charge modulation can occur at intermediate doping where they are preferred over stripes. By analyzing the energy dependence of the LDOS modulation as observed in STM, we have found the data of Ref.  to be well described by combined modulations in charge density as well as exchange and pairing energy, caused by a modulation of the pairing interaction.

###### Acknowledgements.

The author thanks S. Davis, E. Demler, A. Polkovnikov, and J. Zaanen for valuable discussions, and in particular S. Sachdev for suggesting a search for plaquette states within the formalism of Refs. . This research was supported by the DFG through SFB 484.