text:
# Spatially Resolved Millimeter Spectroscopy of the Gravitational Lens PKS 1830-211

## 1 Introduction

PKS 1830-211 is a radio bright source determined to be a gravitationally lensed quasar from its ring-like morphology (Jauncey et al., 1991; Subrahmanyan et al., 1990; Rao & Subrahmanyan, 1988) toward which molecular absorption has been detected (Wiklind & Combes, 1996). The object consists of two flat spectrum images connected by a nearly circular, low surface brightness steep spectrum Einstein ring (Jauncey et al., 1991; Kochanek & Narayan, 1992). The two bright sources are separated by $`1^{\prime \prime }`$ and lie to the northeast (NE) and southwest (SW) of the ring center.

Wiklind and Combes (1996) discovered twelve different rotational transitions involving five different molecular species toward PKS 1830-211 in unresolved spectra. The lines had widths of $`30`$ km s-1, depths of $`40`$%, and were redshifted to $`z=0.89`$. Based on standard isotopic ratios, the detection of H13CO$`{}_{}{}^{+}(21)`$ implied that the HCO$`{}_{}{}^{+}(21)`$ was optically thick. This assumption along with published flux ratios of the NE and SW components led to the conclusion that the absorption was only in the southwest component. Spectra of the two main components were obtained with the Berkeley Illinois Maryland Association (BIMA) array by resolving the source in the 3-mm band (Frye et al., 1996). The molecular absorption only appeared in the SW component as was expected, but did not reach the zero flux level. These results seemed to conflict with an analysis done by Wiklind and Combes (1998) on the unresolved spectrum of PKS 1830-211 where the shift in the field phase center at the frequency corresponding to the HCO$`{}_{}{}^{+}(21)`$ transition was consistent with saturated absorption in the SW image.

High-resolution HCO$`{}_{}{}^{+}(21)`$ spectra of the unresolved source taken with the Plateau de Bure interferometer (Wiklind & Combes, 1998) showed that the trough of the absorption is not flat as would be expected from complete saturation. The line profile was also clearly asymmetric with the broader wing toward negative velocities. These results are reasonable for absorption due to a multi-component cloud with possible small scale structure rather than one optically thick cloud.

Galactic models (Wiklind & Combes, 1998), a determination of the quasar redshift, $`z=2.5`$ (Lidman et al., 1998), and lensing models (Nair et al., 1993; Kochanek & Narayan, 1992) combined to form a coherent picture of this gravitational lens system and the Hubble parameter was determined to be $`H_0=65_9^{+15}`$ km s-1 Mpc-1 (Lidman et al., 1998). The rather large uncertainty in the time delay measurement, about 20% (Lovell et al., 1998), is the major source of error in this estimate of $`H_0`$. Continued monitoring of the flux ratio of PKS 1830-211 would allow for a better estimate of the time delay and therefore a better determination of the Hubble parameter. The flux ratio can be easily determined from an unresolved spectrum of the source if the absorption in the SW image is saturated (Wiklind & Combes, 1996). While many lines have been studied in this system, high-resolution spectra of saturated absorption in the resolved SW image have yet to be presented.

This paper presents observations of PKS 1830-211 made with the BIMA interferometer in Hat Creek, California (Welch et al., 1996) which yielded spectra of the two individual components. There have been a couple of important improvements to the BIMA array since the published results of Frye et al. (1996). The most significant upgrade has been the installation of new SIS mixers in 1997 which has halved the receiver noise temperatures (Engargiola & Plambeck, 1998). There have also been new A-array stations added which produce a better imaging beam. The observations presented here were done with higher spectral resolution than the previous BIMA observations and were intended to give a more detailed look at the line profile of the resolved SW image. In § 2 the observations and observing conditions are outlined and the details of the self-calibration procedure used in the data reduction are given. Section 3 presents flux values and an analysis of the spectra. Rotation temperatures, column depths, and the implications of the line profiles are then discussed in § 4.

## 2 Observations and Data Reduction

PKS 1830-211 was observed in the 3-mm band using the “A-” configuration of the BIMA array which is the standard A configuration without the longest baselines at 1.9 km. The correlator was set so that in each sideband there were two spectral windows containing 128 channels with a spectral resolution of 1.2 km s-1 and four wide band windows each containing 32 channels and spanning 100 MHz. The first local oscillator was tuned so that the redshifted ($`21`$) rotational transitions of HCO+ and HCN, at 94.588 GHz and 93.977 GHz respectively, fell into the two high resolution windows. The HCO+ line was Doppler tracked throughout the observations. The standard quasar 3C273 was observed to calibrate the passband, Uranus gave an absolute flux scale, and 1733-130 was used as a secondary calibrator. PKS 1830-211 was observed in 30 minute intervals using 23 second integrations with the intention of self calibration.

Data were used from two separate observing tracks on 1999 December 27 and 28. The observations on the 27th were good with typical system temperatures between 200 K and 300 K and the RMS phase variations measured on 100 m baselines were on the order of 200 $`\mu `$m. The atmospheric conditions on the 28th were even better with slightly lower system temperatures and phase variations of $`80\mu `$m.

PKS 1830-211 is a bright source which can be self calibrated on a record-by-record basis. The data on each day was split into two datasets corresponding to the upper and lower side bands and the phases were first corrected in the lower side band using an iterative self-calibration process. The first iteration used a point source at the observing center as a model, and produced a two-component image. The inverted CLEAN model containing two main components was then used as the self-calibration model in the subsequent iteration. Two self-calibration iterations were needed for the image to converge. The upper side band was then phase calibrated in the same way with the only difference being that the first model used in the self-calibration iteration was the CLEAN model from the lower side band. The spectral windows were then phase calibrated channel by channel using a 65 minute time interval allowing enough signal for a reliable calibration. The bandpass was determined from 3C273 and was applied to the source and then the datasets were scaled to the flux of Uranus measured on December 27. The calibrated datasets were then combined and inverted into a single plane upper side band continuum map and two 128 plane spectral cubes. All reductions were done using the MIRIAD software package.

The two main components were not resolved to a sufficient degree in all of the data due to the position angle of the synthesized beam and therefore some data were discarded. A local sidereal time range of 1700–1900 was used on both nights covering an elevation range from 17 to 27. The resulting synthesized beam dimensions are $`1\mathrm{}.08\times 0\mathrm{}.59`$ FWHM with a position angle of $`11.6^{}`$. This resolution issue and the implications will be considered further in the following section.

## 3 Analysis and Results

The combined datasets yield a NE continuum level of $`1.13\pm 0.1`$ Jy beam-1 and a SW level of $`0.96\pm 0.1`$ Jy beam-1 with the flux levels of the individual datasets agreeing to within 10%. Although this is known to be a highly variable source (van Ommen et al., 1995), these levels are reasonable in comparison to previous results (Wiklind & Combes, 1998; Frye et al., 1996; Wiklind & Combes, 1996). The flux ratio of the two images, $`1.18\pm 0.06`$, is expected to be more accurate than the overall flux scale and agree to within 3% on the two days. The separation angle was determined using the inverted CLEAN model of the combined datasets giving a value of $`0\mathrm{}.99\pm 0\mathrm{}.05`$, consistent with previously published values (Rao & Subrahmanyan, 1988; Subrahmanyan et al., 1990; Nair et al., 1993; Johns et al., 1996; Frye et al., 1996).

Spectra were taken from the two 128 plane data cubes containing the HCO$`{}_{}{}^{+}(21)`$ and HCN$`(21)`$ transitions at the locations of the continuum peaks (see Fig. Spatially Resolved Millimeter Spectroscopy of the Gravitational Lens PKS 1830-211). The line profiles obtained here are in good qualitative agreement with the high resolution (0.5 km s-1) spectrum of HCO$`{}_{}{}^{+}(21)`$ in Fig. 7 of Wiklind and Combes, 1998 (hereafter 7WC). The absorption widths are comparable at about 45 km s-1 FWHM, and the asymmetry in the lines is apparent with the broader side toward negative velocities. There is some residual flux in the trough of the absorption, however there is no apparent slope as seen in 7WC.

The absorption in the SW image reaches the zero level with 5 channels at or below zero flux in the HCO+ window and 3 channels in the HCN window. In reference to 7WC, a numerical average over 20 km s-1 (16 channels) around the centroid of the lines was computed giving an average flux level in the absorption trough. This procedure applied to the high-resolution profile of HCO$`{}_{}{}^{+}(21)`$ in 7WC yielded a value of $``$ 0.04–0.06 Jy due to the slope of flux in the absorption trough assuming that the negative velocity edge of the the narrow profile corresponds to a zero flux level in the SW image. The NE and SW components are known to have sizes at or below the milliarcsecond scale (Johns et al., 1996), however the two images in our dataset did not deconvolve into two point sources. A no noise, no phase error model of our observations revealed flux contamination between the two images at the level of 0.01 Jy beam-1 preventing a double point source deconvolution. Therefore the flux of a saturated line averaged across the narrow trough given the profile of 7WC and our resolution determined from models is expected to be between 0.05–0.07 Jy beam-1. The average flux levels measured in the absorption troughs of the present data are 0.058 and 0.066 Jy beam-1 for the HCO+ and HCN lines respectively.

The data show no evidence of absorption in the NE component in HCO+ however a small absorption feature is noticeable in the HCN window at 0 km s-1. This feature appeared only in the December 28 data and is small but statistically significant reaching the $`3\sigma `$ level. It is interpreted as a small amount of contamination from the SW component.

## 4 Discussion

This dataset is unique in showing direct evidence of spatially and spectrally resolved saturated absorption in PKS 1830-211. The noise in our spectra, with a measured RMS of 0.098 Jy beam-1, is too high to verify the upward slope from negative to positive velocities in the HCO+ window which is seen in 7WC. However, the absorption does reach zero flux levels in multiple channels in each window and the profile is consistent with saturation given the profile in 7WC and the observational models constructed.

These results are contrary to the results obtained by Frye et al. (1996) who also used the BIMA array to resolve the two components, but obtained absorption spectra that were only $`80\%`$ saturated. The data quality and the judiciousness in setting the self-calibration parameters determine the reliability of an iterative self-calibration process. Lower data quality leads to longer, more elaborate calibrations which often require extensive flagging and ultimately challenge the soundness of the result. Starting with a self-calibration model that doesn’t represent the known spatial distribution of flux only makes the process more unstable. Our data were taken in better conditions than the previous BIMA observations and with new SIS receivers, both of which contributing to a higher overall data quality. There were also different A-array antenna stations producing a better beam from which this image based reduction may have benefited. The calibration process involved only two iterations and varying the self-calibration parameters, such as the calibration interval and convergence criteria, produced consistent results. All of these considerations give us much confidence in the robustness of our result.

The residual flux in the absorption trough is thought to be due to small scale structure in the absorbing cloud which translates into a varying covering factor as a function of velocity (Wiklind & Combes, 1998, 1996; Frye et al., 1996). The asymmetry seen in the lines is also evidence of a complex, inhomogeneous absorbing medium. However, in order to get any quantitative information out of these profiles it is necessary to make simplifying assumptions about the equilibrium state of the gas and the structure.

The deep absorption prevents a direct determination of the physical parameters of the absorbing medium, however limits can be determined using a lower limit of the velocity integrated optical depth. The HCO$`{}_{}{}^{+}(21)`$ and HCN$`(21)`$ line profiles presented here with their $`(10)`$ counterparts presented in Menten et al. (1999) and Carilli et al. (1998) give low values for the rotation temperatures, $`6`$ K, consistent with previous analyses (Wiklind & Combes, 1996; Menten et al., 1999). This result implies that HCN and HCO+ are primarily excited by the CMB, so we can assume that the rotation temperature of these molecules are $`T_{rot}5`$ K which is the temperature of the CMB at $`z=0.89`$. This value of $`T_{rot}`$ along with the calculated velocity integrated optical depths give lower limits to the column densities for HCN and HCO+; $`N(`$HCN$`)3\times 10^{14}`$cm -2 and $`N(`$HCO$`{}_{}{}^{+})2\times 10^{14}`$ cm-2.

The lensing galaxy in which this absorption is taking place is likely a massive, early-type spiral (Wiklind & Combes, 1998, 1996). Given that the abundance ratios are somewhat similar to the Milky Way, the column densities derived above imply an H2 column of $`N(`$H$`{}_{2}{}^{})10^{22}`$ but the low average excitation temperatures of the HCN and HCO+ lines imply a moderately low average volume density ($`10^4`$cm-3). Therefore the pathlength through the cloud must be greater than 1 pc, but it is not clear how extended the cloud is in the line of sight direction. The full velocity range over which there is absorption is $`>60`$ km s-1 which is not atypical for observations in the Milky Way (Frye et al., 1996). However according to galactic (Wiklind & Combes, 1998) and lensing (Nair et al., 1993; Kochanek & Narayan, 1992) models, the lensing galaxy of PKS 1830-211 is oriented almost face on ($`i16^{}`$) and the light from the SW image is traversing the lensing galaxy at a galactocentric distance between 1.8 kpc and 3 kpc. Assuming this orientation, it is hard to imagine how the absorption occurs over such a large velocity range (see Menten et al. 1999).

## 5 Summary

The two main components of the gravitational lens PKS 1830-211 have been resolved with the BIMA interferometer. The continuum levels are $`1.13\pm 0.1`$ Jy beam-1 for the NE component and $`0.96\pm 0.1`$ Jy beam-1 for the SW component giving a flux ratio of $`1.18\pm 0.06`$. These levels are expected to vary, but are consistent with other published values.

High-resolution spectra (1.2 km s-1) of the redshifted lines of HCO$`{}_{}{}^{+}(21)`$ and HCN$`(21)`$ were shown to reach the zero flux level in the SW image but with a non zero flux level averaged across the narrow trough. These results are an improvement over the previous BIMA observations of this object (Frye et al., 1996) and are in agreement with what would be expected from saturated absorption by a complex absorbing cloud with multiple components and small scale structure.

Comparing our line profiles with previously published spectra of the HCO$`{}_{}{}^{+}(10)`$ and HCN$`(10)`$ transitions in the SW component (Menten et al., 1999; Carilli et al., 1998) gives a low upper bound on the excitation temperature, $`6`$ K, which translates to fairly high column densities, $`N(`$HCN$`)3\times 10^{14}`$ cm-2 and $`N(`$HCO$`{}_{}{}^{+})2\times 10^{14}`$ cm-2. The high column densities imply a large H2 column but the low excitation temperatures suggests that the volume density of H2 is moderately low. Given the estimated galactocentric distance of the absorbing cloud $`1.8`$ kpc and the inclination angle of the galaxy, $`i16^{}`$, the line widths are remarkably large.

We would like to thank Dick Plambeck for helping with the many questions that arose during this research. Also thanks to Mel Wright for help with generating model observations and for general advice. Michiel Hogerheijde and Geoff Marcy both made helpful contributions. BLF would like to thank Joe Silk and Tom Broadhurst for their ideas and contributions. This research has been funded by NSF grant AST-9981308.


format: text

text:
# Theoretical predictions on the clustering of SCUBA galaxies and implications for small-scale fluctuations at sub-mm wavelengths

## 1 Introduction

The advent of the 450-850 $`\mu `$m Submillimetre Common User Bolometric Array (SCUBA) camera (Holland et al. 1999) has opened the possibility of exploiting the enormous potential of sub-mm astronomy for investigating the early phases of galaxy evolution, taking advantage of the strongly negative K-correction which boosts the visibility of high-$`z`$ galaxies (Smail, Ivison & Blain 1997; Barger et al. 1998; Hughes et al. 1998; Barger, Cowie & Sanders 1999; Eales et al. 1999; Blain et al. 1999a).

A variety of models has been proposed to reproduce the properties of the submillimetre population (e.g. Guiderdoni et al. 1998; Blain et al. 1999b; Devriendt, Guiderdoni & Sadat 1999; Devriendt & Guiderdoni 2000; Takeuchi et al. 2000; Pearson 2000), in particular the dramatic star-formation rate inferred from the data.

Granato et al. (2000 - hereafter Paper I) propose a unified scheme for the formation of the QSO and spheroidal population. According to their model, at the time of the onset of the QSO activity, the host galaxy has already undergone a major star-formation event, with a characteristic time $`T_B`$ ranging from $`0.5`$ to $`2`$ Gyr when going from more to less massive hosts. The QSO onset practically stops the star-formation process, and when eventually the reservoir of gas fueling the Active Galactic Nuclei is exhausted, the host galaxy will evolve in a passive way, with a spectrum becoming quickly red. This approach identifies SCUBA galaxies as the progenitors of QSO hosts and provides an excellent description of (amongst others) the observed 850 $`\mu `$m number counts.

In this paper we present predictions on clustering properties of SCUBA galaxies in the framework of the Granato et al. (2000) model. Such properties are determined by the combination of the evolution of the underlying mass fluctuations and the bias relating galaxy overdensities to mass. Since the bias factor is determined by properties related to the process of galaxy formation, it follows that measurements of clustering can lead to very precious insights on the galaxy population producing the signal.

Unfortunately the number of SCUBA galaxies observed is still too small to allow for any clustering measurement. The model by Granato et al. (2000) however establishes a tight link between this population and Lyman-break galaxies (Steidel et al. 1996; LBGs hereafter), interpreted as lower-mass/lower-submm luminosity sources. Such model can thus be tested against measurements of LBG clustering (Giavalisco et al. 1998). A further test is provided by the recent detection (Daddi et al. 2000) of strong clustering of Extremely Red Objects (EROs), which were found to be the optical counterparts of two of the still few SCUBA sources having reliable optical identifications (Smail et al. 1999).

An important issue related to the above discussion is the effect of clustering on small scale anisotropies at mm/sub-mm wavelengths, where experiments to map the Cosmic Microwave Background either have been (e.g. BOOMERanG, MAXIMA, TopHat, ACBAR, ARCHEOPS) or are about to be carried out (notably the ESA’s Planck mission). Small scale fluctuations of the background intensity due to clustering of unresolved sources are proportional to the amplitude of their two-point correlation function and to their volume emissivity (see, e.g., De Zotti et al. 1996). The deep SCUBA surveys have demonstrated that the bulk of the extragalactic background at $`850\mu `$m comes from a relatively small number of sources (at least ten time less than those accounting for the bulk of the optical/UV background) at substantial redshifts. Combining the (rather limited) direct spectroscopic information (e.g. Frayer 2000) with estimates based on the sub-mm to radio spectral index (Carilli & Yun 1999), it is concluded that the median redshift of the sub-mm population brighter than $`1`$mJy is probably $`2.5`$–3 (Smail et al. 2000). It then follows that most of the background comes from a population of extremely luminous (hence presumably very massive) high-$`z`$ galaxies which probably correspond to the rare, high density peaks in the primordial matter distribution. The standard structure formation scenario predicts that these peaks are strongly clustered (Kaiser 1984; Baron & White 1987). Also, according to the currently favoured view, SCUBA sources probably correspond to the phase of intense star formation in large spheroidal galaxies, a process that appears to be completed for $`z>1`$; as a consequence, the relevant redshift range is limited and dilution of the clustering signal is relatively small.

The amplitude of intensity fluctuations due to clustering is correspondingly expected to be large. Scott & White (1999) worked out a preliminary estimate adopting for sub-mm sources the angular correlation function determined for $`z3`$ LBGs (which however are mostly undetectable by SCUBA; see Chapman et al. 2000). Here we provide quantitative estimates of this effect using a technique developed by Peacock & Dodds (1996), Matarrese et al. (1997) and Moscardini et al. (1998), which takes into account the effects of the past light-cone, the linear and non-linear growth of clustering and the evolution of the bias factor. The epoch-dependent mass distribution of sources is modelled following Granato et al. (2000).

The plan of the paper is as follows. Section 2 will be devoted to the analysis of the clustering of SCUBA galaxies and comparisons with the observed clustering properties of LBGs and EROs, while in Section 3 we use the findings of Section 2 to make predictions on the temperature fluctuations produced by unresolved SCUBA-like sources and to discuss their implications for experiments aimed at mapping the CMB at sub-mm wavelengths. Section 4 summarizes our conclusions.

Throughout the paper we will assume a cosmological model with $`h_0=0.7`$, $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Lambda }=0.7`$, and adopt the normalization of the rms mass fluctuation $`\sigma _8=1`$ determined by the 4-year COBE data (Bunn & White 1997).

## 2 Clustering of SCUBA sources

The angular two-point correlation function $`w(\theta )`$ is defined as the excess probability, $`\delta P`$, with respect to a random Poisson distribution, of finding two sources in the solid angles $`\delta \mathrm{\Omega }_1`$ $`\delta \mathrm{\Omega }_2`$ separated by an angle $`\theta `$, and it is defined by

$`\delta P=n^2\delta \mathrm{\Omega }_1\delta \mathrm{\Omega }_2\left[1+w(\theta )\right]`$ (1)
where $`n`$ is the mean number density of objects in the particular catalogue under consideration.

The theoretical expression for $`w(\theta )`$ can be derived from its spatial counterpart $`\xi `$ by projection via the relativistic Limber equation (Peebles 1980):

$`w(\theta )=2{\displaystyle \frac{_0^{\mathrm{}}_0^{\mathrm{}}N^2(z)b_{\mathrm{eff}}^2(M_{\mathrm{min}},z)(dz/dx)\xi (r,z)𝑑z𝑑u}{\left[_0^{\mathrm{}}N(z)𝑑z\right]^2}},`$ (2)
where $`x`$ is the comoving radial coordinate, $`r=(u^2+x^2\theta ^2)^{1/2}`$ (for a flat universe and in the small angle approximation: $`\theta 1`$rad), and $`N(z)`$ is the number of objects within the shell ($`z,z+dz`$).

The mass-mass correlation function $`\xi (r,z)`$ to be inserted in eq.(2) has been obtained following the work by Peacock & Dodds (1996, but also see Magliocchetti et al., 2000), extended by Matarrese et al. (1997) and Moscardini et al. 1998) to redshifts $`z>0`$, which provides an analytical way to derive the trend of $`\xi (r,z)`$ both in the linear and non-linear regime. Note that $`\xi (r,z)`$ only depends on the underlying cosmology and on the normalization of $`\sigma _8`$. The relevant properties of SCUBA galaxies are included in the redshift distribution of sources $`N(z)`$, where $`N(z)dz`$ gives the number of objects in the shell $`(z,z+dz)`$ (given in Paper I and shown in Fig. 1 for fluxes respectively greater than 1, 10, 50 mJy), and in the bias factor $`b_{\mathrm{eff}}(M_{\mathrm{min}},z)`$.

The effective bias factor $`b_{\mathrm{eff}}(M_{\mathrm{min}},z)`$ of all haloes with masses greater than some threshold $`M_{\mathrm{min}}`$ is obtained by integrating the quantity $`b(M,z)`$ \- representing the bias of individual haloes of mass $`M`$ \- opportunely weighted by the number density $`n_{SCUBA}(M,z)`$ of SCUBA sources:

$`b_{\mathrm{eff}}(z)={\displaystyle \frac{_{M_{\mathrm{min}}}^{\mathrm{}}𝑑Mb(M,z)n_{\mathrm{SCUBA}}(M,z)}{_{M_{\mathrm{min}}}^{\mathrm{}}𝑑Mn_{\mathrm{SCUBA}}(M,z)}}.`$ (3)
Note that as $`n_{\mathrm{SCUBA}}`$ can be thought as the fraction of haloes hosting a galaxy in the process of forming stars, its expression can also be written as $`n_{SCUBA}(M,z)=n(M,z)T_B/t_h`$, where $`n(M,z)dM`$ is the mass spectrum of haloes with masses between $`M`$ and $`M+dM`$ (Press & Schechter, 1974; Sheth & Tormen, 1999), $`T_B`$ is the duration of the star-formation burst and $`t_h`$ is the life-time of the haloes in which these objects reside (see Martini & Weinberg, 2000). The duration of the star-formation phase $`T_B`$ varies with the mass of the spheroid undergoing the burst (see PaperI), assuming a value $`T_b=0.5`$ Gyr in the case of massive objects (i.e. with bulge masses $`M_{sph}>\mathrm{1.5\hspace{0.33em}10}^{11}M_{}`$) and increasing when going to lower masses.
We take the expression for $`b(M,z)`$ to be plugged into eq.(3) from Jing (1998):

$`b(M,z)=(1+{\displaystyle \frac{1}{\delta _c}}\left[{\displaystyle \frac{\delta _c^2(z)}{\sigma ^2(M)}}1\right])\left({\displaystyle \frac{\sigma ^4(M)}{2\delta _c^4(z)}}+1\right)^{\frac{(62n_{\mathrm{eff}})}{100}}`$ (4)
with $`\delta _c(z)=\delta _c/D(z)`$, $`D(z)`$ being the linear growing factor and $`\delta _c`$ the critical density for collapse of a homogeneous spherical perturbation ($`\delta _c=1.686`$ in an Einstein-de Sitter universe). $`\sigma ^2(M)`$ is the variance of density fluctuations on the mass scale $`M`$, and $`n_{\mathrm{eff}}=36(d\mathrm{ln}\sigma /d\mathrm{ln}M)`$. Note that the assumption of a linear bias is justified by the redshift distribution of the sources, all appearing for $`z>0.75`$ (see e.g. Somerville et al., 2001).

According to Granato et al. (2000), sources showing up in the SCUBA counts can be broadly divided into three categories: low-mass (masses in the range $`M_{\mathrm{sph}}10^910^{10}M_{}`$, duration of the star formation burst $`T_B2`$ Gyr, and typical fluxes $`S<1`$ mJy), intermediate-mass ($`M_{\mathrm{sph}}10^{10}`$$`10^{11}M_{}`$ and $`T_B1`$ Gyr) and high-mass ($`M_{\mathrm{sph}}>10^{11}M_{}`$, $`T_B0.5`$ Gyr, dominating the counts at fluxes $`S>510`$ mJy). Note that by $`M_{\mathrm{sph}}`$ we denote the mass in stars at completion of the star formation process.
In order to evaluate the bias factor in eq.(3) we then consider two extreme cases for the ratio between the mass in stars and the mass of the host dark halo: $`M_{\mathrm{halo}}/M_{\mathrm{sph}}=100`$ and $`M_{\mathrm{halo}}/M_{\mathrm{sph}}=10`$. Note that $`M_{\mathrm{halo}}/M_{\mathrm{sph}}=10`$ roughly corresponds to the ratio $`\mathrm{\Omega }_0/\mathrm{\Omega }_{\mathrm{baryon}}`$ between total and baryon density, where we adopted for the latter quantity the standard value from primordial nucleosynthesis; this corresponds to having assumed all the baryons to be locked into stars and, as a consequence, gives a conservative lower limit to the quantity $`M_{\mathrm{halo}}/M_{\mathrm{sph}}`$. $`M_{\mathrm{halo}}/M_{\mathrm{sph}}=100`$ is instead related to $`\mathrm{\Omega }_0/\mathrm{\Omega }_{}`$, $`\mathrm{\Omega }_{}`$ being the present mass density in visible stars (Persic & Salucci, 1992; Fukugita et al. 1998). The likely value should be $`M_{\mathrm{sph}}/M_{\mathrm{halo}}=13`$ % (cf. Paper I).

Armed with the above results we can finally evaluate the two-point correlation function in eq. (2) for different $`M_{\mathrm{halo}}/M_{\mathrm{sph}}`$ ratios and different flux cuts. Figure 2 presents our predictions for $`w(\theta )`$, respectively for a flux cut of 50 (solid line), 10 (dashed line) and 1 (dotted line) mJy. Higher curves of each kind correspond to the case $`M_{\mathrm{halo}}/M_{\mathrm{sph}}=100`$, while lower curves refer to the $`M_{\mathrm{halo}}/M_{\mathrm{sph}}=10`$ one.

The highest clustering amplitude is found for the brightest sources ($`S50`$ mJy). This is because they are associated to the most massive dark halos and are therefore highly biased tracers of the dark matter distribution. In addition, according to the model by Granato et al. (2000), they have a rather narrow redshift distribution (see Fig. 1), so that the dilution of the clustering signal is minimum.
The very sharp drop of all the curves at $`\theta 1^{}`$ is due to the absence of nearby ($`z<1`$, see Fig. 1) objects. This reflects the notion that the actively star-forming phase in spheroids is completed at $`z>1`$.

The model introduced in PaperI establishes a strong link between SCUBA sources and LBGs where this latter population is predicted to be the low- to intermediate-mass tail of primeval spheroidal galaxies, with $`T_B12`$ Gyr and a star-formation rate ranging from a few to a hundred $`M_{}\text{yr}^1`$.

Measurements of the angular correlation function of LBGs (Giavalisco et al. 1998) then offer a way to test such model. In Fig. 3 we plotted the expected $`w(\theta )`$ for sources with $`S1`$ mJy (corresponding to $`M_{\mathrm{sph}}>10^{10}M_{}`$) within the redshift range $`2.5z3.5`$, covered by the Giavalisco et al. (1998) sample. As in Fig. 2, the higher curve is for $`M_{\mathrm{halo}}/M_{\mathrm{sph}}=100`$, the lower one for $`M_{\mathrm{halo}}/M_{\mathrm{sph}}=10`$. The data points show the Giavalisco et al. (1998) measurements. Despite the large errors, it is clear that the predicted trend for $`w(\theta )`$ correctly reproduces the data, a better fit obtained for $`M_{\mathrm{halo}}/M_{\mathrm{sph}}=100`$. This result, apart from providing a further proof of the tight link between SCUBA galaxies and LBGs, also confirms the expectations of Paper I for a small fraction (on the order of a few percent) of the total mass to be confined into stars.

Strong clustering of EROs has been recently detected by Daddi et al. (2000). This result may again provide a relevant test for our model since it seems likely for a substantial fraction of SCUBA galaxies to be EROs with $`RK6`$: at least two out of nine ($`>20\%`$) sources in the Cluster Lens Survey belong to this class (Smail et al. 1999), while other submillimeter surveys have discovered additional EROs as counterparts of SCUBA sources (Ivison et al. 2000). For their ERO sample with $`K_s<19.2`$ and $`RK_s5.3`$, Daddi et al. (2000) find $`w(\theta )0.07`$ for $`\theta 0^{}.1`$, not far from the predictions in Fig. 2.

## 3 Small Scale Fluctuations from Discrete Sources

### 3.1 The Intensity Correlation Function

An issue intimately connected with the analysis of galaxy clustering performed in the previous Section is the study of the contribution of unresolved sources (i.e. sources with fluxes fainter than some detection limit $`S_d`$) to the background intensity. Its general expression is:

$`I={\displaystyle _0^{S_d}}{\displaystyle \frac{dN}{dS}}S𝑑S`$
$`={\displaystyle \frac{1}{4\pi }}{\displaystyle _{L_{\mathrm{min}}}^{L_{\mathrm{max}}}}𝑑\mathrm{log}LL{\displaystyle _{z(S_d,L)}^{z_{\mathrm{max}}}}𝑑z\mathrm{\Phi }(L,z){\displaystyle \frac{K(L,z)}{(1+z)^2}}{\displaystyle \frac{dx}{dz}}`$ (5)
(see e.g. De Zotti et al. 1996), where $`dN/dS`$ denotes the differential number counts, $`L_{\mathrm{max}}`$ and $`L_{\mathrm{min}}`$ are respectively the maximum and minimum local luminosity of the sources, $`K(L,z)=(1+z)L[\nu (1+z)]/L(\nu )`$ is the K-correction, $`z_{\mathrm{max}}`$ is the redshift when the sources begin to shine, $`z(S_d,L)`$ is the redshift at which a source of luminosity $`L`$ is seen with a flux equal to the detection limit $`S_d`$, $`\mathrm{\Phi }(L,z)`$ is the luminosity function (i.e. the comoving number density of sources per unit $`d\mathrm{log}L`$), and $`x`$ is the comoving radial coordinate.

The intensity fluctuation $`\delta I`$ due to inhomogeneities in the space distribution of unresolved sources is then given by eq. (5), with the quantity $`\mathrm{\Phi }(L,z)`$ replaced by $`\delta \mathrm{\Phi }(L,z)`$. It is easily shown that the angular correlation of such intensity fluctuations

$`C(\theta )=\delta I(\theta ^{},\varphi ^{})\delta I(\theta ^{\prime \prime },\varphi ^{\prime \prime }),`$ (6)
where $`(\theta ^{},\varphi ^{})`$ and $`(\theta ^{\prime \prime },\varphi ^{\prime \prime })`$ define two positions on the sky separated by an angle $`\theta `$, can be expressed as the sum of two terms $`C_P`$ and $`C_C`$, the first one due to Poisson noise (i.e. fluctuations given by randomly distributed objects), and the second one owing to source clustering.

Scott & White (1999) have shown that if dusty galaxies cluster like LBGs (Giavalisco et al. 1998), anisotropies due to clustering dominate the Poisson ones at all angular scales. According to our model, LBGs are the optical counterparts of the low- to intermediate- mass tail of SCUBA galaxies; as illustrated in Figure 2, this makes the clustering signal of LBGs a factor $`10`$ lower than the one obtained for SCUBA galaxies; therefore we can safely assume the Poisson term $`C_P`$ to be negligible in the following analysis and only concentrate on temperature fluctuations caused by the $`C_C`$ term (hereafter simply called $`C`$).

By making use of the quantities defined in Section 2 and of eq. (5), the clustering term in eq. (6) takes the form:

$`C(\theta )=\left({\displaystyle \frac{1}{4\pi }}\right)^2{\displaystyle _{z_{(L_{\mathrm{min}},S_d})}^{z_{\mathrm{max}}}}𝑑zb_{\mathrm{eff}}^2(z){\displaystyle \frac{j_{\mathrm{eff}}^2(z)}{(1+z)^4}}\left({\displaystyle \frac{dx}{dz}}\right)^2{\displaystyle _0^{\mathrm{}}}𝑑u\xi (r,z),`$ (7)
where

$`j_{\mathrm{eff}}={\displaystyle _{L_{\mathrm{min}}}^{\mathrm{min}[L_{\mathrm{max},L(S_d,z)}]}}\mathrm{\Phi }(L,z)K(L,z)L𝑑\mathrm{log}L`$ (8)
is the effective volume emissivity (see Toffolatti et al. 1998).

Following the prescriptions of Section 2, we have then evaluated $`C(\theta )`$ in eq.(7) separately for the three cases of low-, intermediate- and high-mass objects, by plugging in eq. (8) the appropriate expressions for the luminosity function. Note that in this case, the quantity $`b^2(M_{\mathrm{min}},z)`$ should indeed be read as $`b^2(M_{\mathrm{min}},M_{\mathrm{max}},z)`$, where $`M_{\mathrm{max}}`$ is the maximum halo mass corresponding to the maximum visible bulge mass (i.e. upper limit for the mass locked into stars whose values are given in Section 2). This corresponds in eq.(3) to a replacement in the upper limit of the integrals of $`\mathrm{}`$ with $`M_{\mathrm{max}}`$.

Figure 4 shows predictions for the correlation signal due to intensity fluctuations (in units ergs2 cm-4 Hz-2 sr-2) in the case of SCUBA galaxies. Dotted, dashed and solid lines are respectively for spheroids of $`M_{sph}^{min}=10^9,\mathrm{\hspace{0.33em}10}^{10}`$ and $`10^{11}M_{}`$. In each plot the two solid curves correspond to different detection limits, the 100 mJy case represented by the higher curve and the 10 mJy case by the lower one. Note that, due to the steepness of the bright end of source counts, the results are independent of $`S_d`$ if $`S_d100`$ mJy. Results are shown for $`M_{\mathrm{halo}}/M_{\mathrm{sph}}=10`$ (left-hand panel) and $`M_{\mathrm{halo}}/M_{\mathrm{sph}}=100`$ (right-hand panel). Note once again the steep drop of all the curves around $`\theta 1^{}`$ caused by the absence of nearby sources of this kind and the very small contribution given by the low-mass objects, mostly due to their small bias function and their modest contribution to the effective volume emissivity.

In order to evaluate the total contribution of the clustering to intensity fluctuations one has to add up all the values of $`C(\theta )`$ in eq. (7) obtained for the different mass intervals and also to take into account the cross-correlation terms between objects of different masses, so that the final expression for $`C(\theta )`$ is given by

$`C^{TOT}(\theta )={\displaystyle \underset{i,j=1}{\overset{3}{}}}\sqrt{C_i(\theta )C_j(\theta )},`$ (9)
where the indexes $`i`$,$`j`$ stand for high, intermediate and low masses.

### 3.2 Power Spectrum of Temperature Fluctuations

The angular power spectrum of the intensity fluctuations \[eqs. (7)–(9)\] can be obtained starting from the well known expression for the expansion of temperature fluctuations in spherical harmonics

$`{\displaystyle \frac{\delta T}{T}}(\theta ,\varphi )={\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=l}{\overset{l}{}}}a_l^mY_l^m(\theta ,\varphi ).`$ (10)
If the fluctuations are a stationary process it can be shown (Peebles 1993) that the previous expression does not depend on the index $`m`$ (i.e. on the $`\varphi `$ angle), so that we can write:

$`C_l={\displaystyle \frac{1}{2l+1}}{\displaystyle \underset{m=l}{\overset{l}{}}}|a_l^m|^2=|a_l^0|^2,`$ (11)
with

$`a_l^0={\displaystyle _0^{2\pi }}{\displaystyle _0^\pi }{\displaystyle \frac{\delta T(\theta )}{T}}\sqrt{{\displaystyle \frac{2l+1}{4\pi }}}P_l(\mathrm{cos}\theta )\mathrm{sin}\theta d\theta d\varphi `$ (12)
(see Toffolatti et al. 1998), where $`P_l(\mathrm{cos}\theta )`$ are the Legendre polynomials, and $`\delta T/T`$ is expressed as

$`{\displaystyle \frac{\delta T}{T}}(\theta )=\left({\displaystyle \frac{\mathrm{\Delta }T}{T}}\right)^2^{1/2}`$ $`=`$ $`{\displaystyle \frac{\lambda ^2\sqrt{C^{TOT}(\theta )}}{2k_bT}}\left[\mathrm{exp}\left({\displaystyle \frac{h\nu }{k_bT}}\right)1\right]^2`$ (13)
$`\times `$ $`\mathrm{exp}\left({\displaystyle \frac{h\nu }{k_bT}}\right)/\left({\displaystyle \frac{h\nu }{k_bT}}\right)^2,`$
which relates intensity and temperature fluctuations (T=2.726 K, Mather et al. 1994).

Figure 5 shows the predicted values for the quantity $`\delta T_l=\sqrt{l(l+1)C_l/2\pi }`$ (in units of K) at 353 GHz ($`850\mu `$m) - the central frequency of one of the channels of the High Frequency Instrument (HFI) of the ESA’s Planck mission - as a function of the multipole $`l`$ up to $`l=1000`$. Results are plotted for two values of the source detection limit ($`S_d=100`$ and 10 mJy; note that the source detection limit for this Planck channel is likely to be $`>100`$mJy) and the usual two values of $`M_{\mathrm{halo}}/M_{\mathrm{sph}}`$. Also shown, for comparison, is the power spectrum of primary (CMB) anisotropies (solid line) predicted by a standard Cold Dark Matter model for a $`\mathrm{\Lambda }`$CDM ($`\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_0=0.3`$, $`h_0=0.7`$), computed with the CMBFAST code developed by Seljak & Zaldarriaga (1996).

At this frequency our model predicts fluctuations of amplitude due to clustering comparable to (and possibly even larger than) those obtained for primary CMB anisotropies at $`l>50`$. As illustrated by Fig. 4, most of the clustering signal comes from massive galaxies with fluxes $`S>1`$ mJy, which lie at substantial redshifts (see Fig. 1) and are therefore highly biased tracers of the underlying mass distribution. Also the strongly negative K-correction increases their contribution to the effective volume emissivity \[eq. (8)\] and therefore to the fluctuations.
We therefore expect this effect to give a substantial contribution to the power spectrum extracted e.g. from the BOOMERanG map (de Bernardis et al. 2000) at 400 GHz. For $`l>200`$ its contribution may be comparable to that due to inhomogeneous emission by the interstellar dust.

## 4 Conclusions

We have worked out predictions for the clustering properties of SCUBA-selected galaxies in the framework of an astrophysically grounded model, relating the formation of QSOs and spheroids (ellipticals, S0 galaxies and bulges of spirals). The theoretical angular correlation function has been derived for different bias functions, corresponding to different values of the ratio $`M_{\mathrm{halo}}/M_{\mathrm{sph}}`$ between the mass of a spheroid locked in stars and the mass of its host halo.

SCUBA sources are predicted to be strongly clustered, with a clustering strength which increases with increasing mass (or equivalently with increasing luminosity), as a consequence of the fact that they are very massive and shining at substantial redshifts. Since the clustering amplitude strongly depends on the quantity $`M_{\mathrm{halo}}/M_{\mathrm{sph}}`$, future measurements of the angular correlation function $`w(\theta )`$ will be able to discriminate amongst different models of SCUBA galaxies and in particular to determine the amount of baryonic mass actively partaking the process of star formation.

Under the hypothesis of low- to intermediate-mass primeval spheroids to show up in the optical band as LBGs (as argued by Granato et al. 2000), we have then compared our predictions with the clustering measurements obtained by Giavalisco et al. (1998) for a wide sample of LBGs at $`z3`$. The agreement is good for a mass ratio between the dark halo and the stellar component in the range 10–100; high values for this ratio seem to be favoured.
Also, the predicted amplitude of the angular correlation function of SCUBA galaxies is consistent with the one determined by Daddi et al. (2000) for EROs with $`RK_s5.3`$ and $`K_s19.2`$. This would support the observational evidence for a substantial fraction of SCUBA galaxies to be identified with EROs (Smail et al. 1999; Ivison et al. 2000).

We have also considered the effect of clustering on maps of the Cosmic Microwave Background at sub-mm wavelengths, with special attention to the 353 GHz ($`850\mu `$m) channel of the High Frequency Instrument (HFI) of the ESA’s Planck mission. The possibility of large fluctuations due to clustering at this (and also shorter) wavelength was first pointed out by Scott & White (1999) who adopted a phenomenological approach, assuming for SCUBA galaxies the clustering amplitudes observed for LBGs (Giavalisco et al. 1996). Our treatment, leading to a derivation of the angular autocorrelation function in the framework of the currently standard hierarchical clustering scenario, indicates that massive SCUBA galaxies (expected to dominate the population of sources sampled at fluxes $`S>10`$ mJy) are likely to be even more strongly clustered than LBGs. The corresponding fluctuations are then expected to be, at $`350`$GHz, of an amplitude comparable with those due to primary CMB anisotropies for multipole numbers $`l>50`$. This would imply that important information on the clustering properties of faint sub-mm galaxies (and hence on their physical properties such as their mass and/or the amount of baryons involved in the star-formation process) resides in the BOOMERanG maps at 400 GHz where, however, the dominant signal is expected to come from interstellar dust emission. On the other hand, if the ratio $`M_{\mathrm{halo}}/M_{\mathrm{sph}}`$ is large, we expect the clustering signal to be comparable with dust emission fluctuations at $`l>200`$. It might therefore be possible to separate these two contributions and possibly exploit the different shapes of their power spectra.

As stressed by Scott & White (1999), the amplitude of the clustering signal drops rapidly - in comparison to the one arising from primary anisotropies - with decreasing frequency. In the nearest Planck/HFI channel (centered at 217 GHz) its contamination on the CMB anisotropy maps should already be negligible.

## ACKNOWLEDGMENTS

Mauro Giavalisco is warmly thanked for providing us with the data shown in Figure 3. We are grateful to Pierluigi Monaco for useful clarifications on the relationship between SCUBA galaxies and quasars. We also thank ASI and Italian MURST for financial support.


text:
# Ferroelectrics with Low-energy Electronic Excitations

## I  Introduction

In the last few years there has been an enormous revival in the study of ferroelectric and ferroelectric materials ,, . In a large part it is due to the investigation of complex materials with very large number of atoms per unit cell and substantial disorder. The disorder also introduces low-energy electronic states. Such materials display what has come to be known as ‘relaxor-ferroelectric’ properties. This term appears to cover all properties, not seen in simple ferroelectrics- for example $`BaTiO_3`$. Generally these materials do not have a ferroelectric transition but do have a very large dielectric-constant, the polarization displays complicated hysteresis loops and has a strong frequency and history dependence at low energies.

Various ideas for the form of disorder and its effect on the ferroelectric transition have been proposed ,,. But to our knowledge, the important effect of low-energy particle-hole excitations due to shallow electronic levels has not been treated. Ferroelectricity (like superconductivity) is rather special because the order parameter, the polarization, couples to the low-energy electronic excitation through the electromagnetic-fields. In this paper we study such effects. The effects of disorder are neglected here by adopting a virtual-crystal aproximation. In fact both effects, disorder and low-energy electronic excitations are important; we hope to treat them together in the future.

This paper is organised as follows. In Section 2, we present a general Lagrangian for a polarization instability coupled to electronic excitations through the electromagnetic fields. The generalization to the disordered case is also indicated. In section 3, we consider the simple case that the low energy excitations are due to a low density of free electrons. While it is common-place knowledge that free-electrons will screen the electric field so that macroscopic fields cannot arise, it is worthwhile to precisely derive the conditions for this and see the condition for the formation of ferroelectric domains at low enough charge density. The domain walls in this case may be charged. We also derive the low energy charge and current response functions for this case. In Sec. 4 we take up the more interesting problem of localised electrons in a Coulomb glass phase and derive the charge response functions.

## II The Lagrangian for soft optic modes in presence of low-energy electronic excitations

We consider the following Lagrangian density (at imaginary times), $`L=L_{ferro}+L_{EM}+L_{int}+L_e`$, where

$$L_{ferro}=\frac{1}{2\chi }\left[(\frac{\stackrel{}{P}}{\tau })^2+\omega _0^2\stackrel{}{P}^2+\gamma (\stackrel{}{P})^2\right]+\frac{1}{4}u\stackrel{}{P}^4,$$
(2)
is a Langrangian describing the polarization vector $`\stackrel{}{P}`$. For $`\omega _0^2=\alpha (TT_c)`$ Eq. (1a) describes an instability of the system (in the absence of coupling to electric fields) when $`TT_c`$.

$$L_{EM}=\frac{1}{8\pi }\left[(\frac{1}{c}\frac{\stackrel{}{A}}{\tau })^2+(\varphi )^2\right]+\frac{1}{8\pi }(\times \stackrel{}{A})^2$$
(3)
is the usual Lagrangian for EM fields in terms of vector potential $`\stackrel{}{A}`$ and scalar potential $`\varphi `$ in the Coulomb gauge, $`.\stackrel{}{A}=0`$. In the microscopic theory, $`\omega _0`$ is the frequency of the zone-center transverse optic phonon. The polarization vector $`\stackrel{}{P}`$ couples to the electric field $`\stackrel{}{E}`$ through the usual Lagrangian $`L_{int}`$, where

$$L_{int}=\stackrel{}{P}.\stackrel{}{E}=\stackrel{}{P}.\left[\frac{1}{c}\frac{\stackrel{}{A}}{\tau }+\varphi \right].$$
(4)
In the presence of charges, the EM potentials pick up additional terms in their Lagrangian. In the Coulomb gauge they are

$$L_e=\frac{1}{2}\left[\chi _L\varphi ^2+(\frac{1}{c})^2\chi _T\stackrel{}{A}^2\right],$$
(5)
where $`\chi _L`$ and $`\chi _T`$ are the longitudinal and transverse response functions of the electrons to scalar and vector potentials $`\varphi `$ and $`\stackrel{}{A}`$, respectively. In Eqs. (1), $`\stackrel{}{P},\varphi `$ and $`\stackrel{}{A}`$ depend on the position $`𝐫`$. For the translational invariant case, this may be directly Fourier-transformed. In the case of disorder, they are random functions of position and the response functions $`\chi _L`$ and $`\chi _T`$ are functions of two co-ordinates. In what follows, we shall assume an effective translational invariant approximation, where disorder is considered averaged over at the outset.

In the Coulomb gauge, the Lagrangian can be separated into longitudinal and transverse parts, $`L=L_L+L_T+L_u`$. We obtain after fourier transforming,

$$L_L=\frac{1}{2}(v(q)^1\chi _L)\varphi ^2\frac{1}{2\chi }[(i\omega )^2(\omega _0^2+\gamma ^2q^2)]\stackrel{}{P}_L.\stackrel{}{P}_Li(\stackrel{}{P}_L.\stackrel{}{q})\varphi ,$$
(7)
where $`v(q)=4\pi /q^2`$ is the Coulomb potential, and

$$L_T=\frac{1}{8\pi c^2}[(i\omega )^2c^2q^2+4\pi \chi _T]\stackrel{}{A}^2\frac{1}{2\chi }[(i\omega )^2(\omega _0^2+\gamma ^2q^2)]\stackrel{}{P}_T.\stackrel{}{P}_Ti\stackrel{}{P}_T.(\frac{i\omega \stackrel{}{A}}{c}),$$
(8)
where we have divided the polarization vector $`\stackrel{}{P}`$ into longitudinal and transverse parts $`\stackrel{}{P}=\stackrel{}{P}_L+\stackrel{}{P}_T`$, and

$$L_u=\frac{u}{4}(\stackrel{}{P}_L+\stackrel{}{P}_T)^4,$$
(9)
is the fourth-order term in $`L_{ferro}`$. Notice that $`L_u`$ mixes the longitudinal and transverse parts of $`\stackrel{}{P}`$.

To study the ferroelectric instability, we take the $`i\omega =0`$ limit and first minimize the free energy with respect to the $`\varphi `$ and $`\stackrel{}{A}`$ fields. We obtain

$`\varphi (q)`$ $`=`$ $`{\displaystyle \frac{i\stackrel{}{q}.\stackrel{}{P}_L}{v(q)^1\chi _L(q,0)}}`$ (10)
$`{\displaystyle \frac{\stackrel{}{A}}{c}}`$ $`=`$ $`4\pi i{\displaystyle \frac{i\omega \stackrel{}{P}_T}{(i\omega )^2c^2q^2+4\pi \chi _T(q,0)}}0.`$ (11)
Substituting these back into Eq. (II), we obtain an effective free energy in terms of $`\stackrel{}{P}`$ only,

$$f=\frac{1}{2\chi }(\omega _o^2+\gamma ^2q^2)[(\stackrel{}{P}_L)^2+(\stackrel{}{P}_T)^2]+\frac{1}{2}\frac{(\stackrel{}{q}.\stackrel{}{P}_L)^2}{[v(q)^1\chi _L]}+\frac{u}{4}\stackrel{}{P}^4.$$
(12)
In an infinite medium $`P_T`$ decouples from EM field completely in the static limit, as is required by Faraday’s Law that $`\dot{\stackrel{}{B}}=\times \stackrel{}{E}0`$ when a transverse electric field exists. This is not true in a finite sample, where a macroscopic $`P_T`$, if it exists is discontinuous at the surface. A surface charge - density $`\sigma =\stackrel{}{P}.\widehat{n}`$ is then generated leading to an electric field.

## III Free electrons

For free electrons, we may use the fact that

$$(v(q)^1\chi _L(q.0))^1\frac{4\pi }{q^2+q_{sc}^2}$$
(13)
in the small $`q`$ limit. In Eq.(5), $`q_{sc}`$ is the static screening length, given in terms of density by the Thomas-Fermi length at low temperatures and by the Debye length in the classical regime at high temperatures, i.e.

$$q_{sc}^2q_{TF}^2=6\pi n_0e^2/E_f=4(3/\pi )^{1/3}n_0^{1/3}/a_0$$
(14)
for $`T<<E_f`$. Here $`n_0`$ is the density of electrons, $`E_f`$ their Fermi-energy and $`a_0`$ the Bohr radius, and

$$q_{sc}^2q_D^2=4\pi n_0e^2/kT.$$
(15)
at high temperature. For most problems of interest in ferroelectrics, the density is low enough that near room temperature the Debye Screening is applicable.

In the absence of free electrons, an instability developes in $`\stackrel{}{P}_T`$ field in the $`\stackrel{}{q}0`$ limit as $`TT_c`$ (or $`\omega _o0`$) always before any possible instability in $`\stackrel{}{P}_L`$. This instability occurs at $`T_c`$, where the transverse optic phonon frequency $`0`$. The magnitude of $`T_c`$ will in general change in the presence of free electrons due to microscopic effects not treated here. The presence of free electrons bring the transition in $`\stackrel{}{P}_T`$ and $`\stackrel{}{P}_L`$ to the same temperature $`T_c`$ because of screening. i.e. in infinite samples, no uniform electric-field can be produced due to a macroscopic $`\stackrel{}{P}`$ because the electric field sets up by the ions oscillation is screened. The electric field in the ordered state of the system with $`<\stackrel{}{P}_T>0`$ is given by

$`\stackrel{}{E}=<\varphi >=\underset{q0}{lim}\left({\displaystyle \frac{q^2\sigma }{v(q)^1\chi _L}}\right).`$Here $`\sigma `$ is the surface charge density at the surface of the sample at which $`\stackrel{}{P}_T`$ has a normal component $`\stackrel{}{P}_T.\widehat{n}0`$. $`\stackrel{}{E}`$ goes to zero as long as $`\chi _L(q,0)`$ is finite, i.e. $`<\stackrel{}{P}_T>0`$ does not lead to uniform electric field in the presence of free electrons. This does not necessarily mean that ferroelectricity is destroyed by arbitrarily small density of electrons. As shown below a domain structure of electric-fields is still possible provided the charge density is low enough.

At temperatures $`T<T_c`$, it is easy to minimize the energy to show that $`<\stackrel{}{P}>=\stackrel{}{P}_o=\sqrt{\frac{\omega _o^2}{\chi u}}`$. The Langrangian can be written as $`L=L(P_o)+L(P^{})`$, where $`\stackrel{}{P}^{}`$ are fluctuations around $`\stackrel{}{P}=\stackrel{}{P}_o`$. We find that to second order in $`P^{}`$, $`L(P^{})=L_L(P_L^{})+L_T(P_T^{})`$, where

$$L_L(P_L^{})=\frac{1}{2}\left[\frac{1}{\chi }((i\omega )^2(2\omega _o^2+\gamma ^2q^2))\frac{4\pi q^2}{q^2+q_{sc}^2}\right](\stackrel{}{P}_L^{})^2$$
(17)
and

$$L_T(P_T^{})=\frac{1}{2}\left[\frac{1}{\chi }((i\omega )^2(2\omega _o^2+\gamma ^2q^2))\frac{4\pi (i\omega )^2}{(i\omega )^2c^2q^2+4\pi \chi _T}\right](\stackrel{}{P}_T^{})^2.$$
(18)
Notice that the only difference between $`T>T_c`$ and $`T<T_c`$ in the quadratic fluctuations in $`\stackrel{}{P}`$ is that $`\omega _o^2`$ is replaced by $`2\omega _o^2`$ as temperature decreases from $`T>T_c`$ to $`T<T_c`$.

## IV domain wall formation

For usual (insulating) ferroelectrics with finite dimensions, domains are formed to minimize the electric field energy $`d^3x\frac{\stackrel{}{E}^2}{8\pi }`$, coming from effective surface charges $`\sigma (\stackrel{}{x})\widehat{n}.\stackrel{}{P}(\stackrel{}{x})`$, where $`\widehat{n}`$ is a surface unit vector. In this section the criteria of domain formation in the presence of free charges is studied. For simplicity we shall consider a ferroelectric with a finite width $`2L`$ in $`z`$-direction, and with infinite extent in $`xy`$ plane. We shall assume a transverse polarization field $`\stackrel{}{P}_T`$ with $`\stackrel{}{P}_T`$ pointing in $`z`$-direction, and with magnitude function of $`x`$ and $`y`$ coordinates, i.e. $`\stackrel{}{P}=\widehat{z}P(x,y)`$ for $`L<z<L`$. For a finite $`L`$ surface charges $`\sigma (x,y)=(\delta (zL)\delta (z+L))P(x,y)`$ will be set up, and the total energy of the system is given in momentum space by

$`F`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{q}}{}}\left[{\displaystyle \frac{1}{2}}(v(q)^1\chi _L)|\varphi (\stackrel{}{q})|^2\sigma (\stackrel{}{q})\varphi (\stackrel{}{q})+{\displaystyle \frac{1}{2\chi }}(\omega _0^2+\gamma ^2q^2)|\stackrel{}{P}(\stackrel{}{q})|^2\right]`$ (20)
$`+{\displaystyle \frac{u}{4}}{\displaystyle _L^L}𝑑z{\displaystyle 𝑑x𝑑y\stackrel{}{P}^4},`$
where $`\sigma (\stackrel{}{q})=(2i)sin(q_zL)P(q_x,q_y)`$ is the fourier transform of the surface charge density. Notice that $`\stackrel{}{q}.\stackrel{}{P}_T(\stackrel{}{q})=0`$ in the sample. Minimizing the energy with respect to $`\varphi `$ we obtain

$`F`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{q}}{}}\left[{\displaystyle \frac{1}{2}}{\displaystyle \frac{(2sin(q_zL)|P(q_x.q_y)|)^2}{v(q)^1\chi _L}}+{\displaystyle \frac{1}{2\chi }}(\omega _o^2+\gamma ^2q^2)|P(q_x,q_y)|^2\right]`$ (22)
$`+{\displaystyle \frac{u}{4}}{\displaystyle _L^L}𝑑x{\displaystyle 𝑑x𝑑y\stackrel{}{P}^4}.`$
Summing over $`q_z`$, and using the results that

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑q_z{\displaystyle \frac{sin^2(q_zL)}{q_z^2+Q^2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑q_z{\displaystyle \frac{1}{2(q_z^2+Q^2)}}={\displaystyle \frac{\pi }{2Q}},`$when $`L\mathrm{}`$, where $`Q^2`$ is a positive number, we obtain

$`F`$ $``$ $`(L){\displaystyle \underset{q_x,q_y}{}}\left[{\displaystyle \frac{\pi }{L}}{\displaystyle \frac{\stackrel{}{P}(\stackrel{}{q})^2}{\sqrt{q^2+q_{sc}^2}}}+{\displaystyle \frac{1}{\chi }}(\omega _o^2+\gamma ^2q^2)\stackrel{}{P}(\stackrel{}{q})^2\right]`$ (24)
$`+{\displaystyle 𝑑x𝑑y\frac{uL}{2}\stackrel{}{P}^4},`$
where $`\stackrel{}{q}=(q_x,q_y)`$ in Eq. (24). Minimizing the free energy with respect to $`P(\stackrel{}{q})`$ and $`\stackrel{}{q}`$ we see that the system has now an instability at finite $`q`$ which indicates the instability of the system to domain formation. Expanding the quadratic terms in $`\stackrel{}{P}`$ in Eq. (24) in powers of $`q^2`$ at small $`q`$, we find that the system has an instability towards domain formation when

$$\frac{\gamma ^2}{\chi }\frac{2\pi }{q_{sc}^3L}<0,$$
(25)
or when $`q_{sc}L_D^1=(\frac{\pi \chi }{2\gamma ^2L})^{\frac{1}{3}}`$. It is easy to show from Eq. (24) that in the absence of free electrons ($`q_{sc}=0`$) $`L_D^1`$ gives the characterestic size of the domains. The criteria for formation of domains (25) can thus be interpreted as simply saying that the screening length of free electrons should be longer than the natural domain size in order for domains to be formed when free electrons are present.

The presence of a finite charge density has important consequences for the switching and hysteresis properties of domains by reversing external-fields. Domain reversal must be accompanied now by motion of charges. The time scale for this is now limited by the time for the charges to diffuse over the characterestic size of the domains.

## V dielectric functions

To obtain the dielectric function at $`T>T_c`$, we neglect the $`L_u`$ term and integrate out the $`\stackrel{}{P}`$ fields first in the Lagrangian (II). As a result, we obtain two effective Lagrangain $`L_{eff}(\varphi )`$ and $`L_{eff}(\stackrel{}{A})`$ for the electromagnetic potentials. The dielectric functions can be obtained by identifying

$`L_{eff}(\varphi )`$ $`=`$ $`{\displaystyle \frac{ϵ_L(q,i\omega )}{2v(q)}}\varphi ^2,`$ (26)
$`L_{eff}(\stackrel{}{A})`$ $`=`$ $`{\displaystyle \frac{ϵ_T(q,i\omega )[(i\omega )^2c^2q^2]}{8\pi }}A_\mu ^2,`$ (27)
where we have again considered the Coulomb gauge $`.\stackrel{}{A}=0`$. After some algebra we obtain

$$ϵ_L(q,i\omega )=1v(q)\chi _L(q,i\omega )\frac{4\pi \chi }{(i\omega )^2\omega _o^2(q)},$$
(29)
and

$$ϵ_T(q,i\omega )=1\frac{4\pi \chi _T(q,i\omega )}{(i\omega )^2}\frac{4\pi \chi }{(i\omega )^2\omega _o(q)^2},$$
(30)
where $`\omega _o(q)^2=\omega _o^2+\gamma ^2q^2`$. A similar expression can also be derived for $`T<T_c`$. The only difference is that $`\omega _o^2`$ is replaced by $`2\omega _o^2`$ when $`T<T_c`$. Notice that $`\omega _o(q)^2`$ is positive definite at all temperatures $`TT_c`$.

The longitudinal electromagnetic collective mode dispersions in the system can be obtained by solving the equation $`ϵ_L(q,\omega (q))=0`$. In the small $`q`$ limit and with $`\omega >>q`$, $`v(q)\chi _L(q,\omega )\frac{\omega _P^2}{\omega ^2}`$, where $`\omega _P`$ is the electronic plasma frequency, and we obtain the secular equation,

$$1\frac{\omega _P^2}{\omega ^2}\frac{4\pi \chi }{\omega ^2\omega _o(q)^2}=0.$$
(31)
(This equation is not modified in the classical limit for long-wavelengths.) Solving the equation we find that there are in general two solutions $`\omega _1`$ and $`\omega _2`$ in the equation. At small $`q`$ and for $`\omega _o(q)^2<<\omega _P^2`$, we obtain

$`\omega _1^2`$ $``$ $`\omega _P^2+4\pi \chi `$ (32)
$`\omega _2^2`$ $``$ $`\left({\displaystyle \frac{\omega _P^2}{\omega _P^2+4\pi \chi }}\right)\omega _o(q)^2.`$ (33)
The transverse electromagnetic collective mode dispersions are obtained by solving $`\omega ^2ϵ_T(q,\omega )c^2q^2=0`$. We obtain after some simple algebra a secular equation very similar to Eq. (31), except that $`\omega _P^2`$ is replaced by $`\omega _P^2+c^2q^2`$, with similar solutions

$`\omega _1^2`$ $``$ $`\omega _P^2+4\pi \chi +c^2q^2`$ (34)
$`\omega _2^2`$ $``$ $`\left({\displaystyle \frac{\omega _P^2+c^2q^2}{\omega _P^2+c^2q^2+4\pi \chi }}\right)\omega _o(q)^2.`$ (35)
Notice that $`\omega _1`$ gives the usual plasmon and polariton modes, respectively for longitudinal and transverse EM fields in metals, with modifications coming from polarization field, whereas $`\omega _2`$ gives the collective modes in longitudinal and transverse polarization fields $`\stackrel{}{P}_L`$ and $`\stackrel{}{P}_T`$, respectively, with modifications coming from metallic behaviour.Notice that at small $`q`$ and when $`\omega _o^20`$, $`\omega _2`$ becomes very small and will be in the hydrodynamic regime where $`\chi _L`$ would have a diffusive form

$`\chi _L(q,\omega )={\displaystyle \frac{dn}{d\mu }}({\displaystyle \frac{Dq^2}{i\omega Dq^2}}).`$In this regime, the dispersion $`\omega _2`$ will be modified. In particular, a finite damping of the collective mode will occur because of coupling to the electronic degrees of freedom.

## VI ferroelectrics with localized electrons

We shall consider here ferroelectrics in the presence of strongly localized electrons in the Coulomb glass phase. In the dilute particle-hole pairs approximation (see Appendix), the localized electrons have only longitudinal response to external electric fields, and with response functions $`\chi _{CG}(q,\omega )`$ given in the Appendix. Their effects on the ferroelectrics can be studied by replacing $`\chi _L`$ for free electrons by $`\chi _{CG}`$. First we consider $`\omega =0`$ and study ferroelectric instability and domain formation at different temperature ranges.

### A ferroelectric instability and domain formation

As a function of temperature, the response function $`\chi _{CG}(q,0)`$ has three different regimes. First we consider high temperature $`k_BT>>\mathrm{\Delta }`$ where the Coulomb glass effect is unimportant and the response of the electrons are dominated by phonon-assisted hopping. This is probably the regime of experimental interests in ferroelectrics. Using Eq. (A23) for $`Re\chi _{CG}(q,0)`$ in Eq.(10) and repeating the analysis afterward we find that the free energy is given by Eq.(11) with $`\sqrt{q^2+q_{sc}^2}`$ replaced by $`q\sqrt{ϵ}`$ in first term of Eq.(11), where $`ϵϵ_o+(\frac{4\pi ^4g^2Aa_o^5}{1440})(\frac{A}{k_BT})^{\frac{1}{4}}`$, where $`ϵ_o`$ includes all the non-singular contributions to $`ϵ`$. The characteristic domain size $`L_D`$ is given by

$$L_D(\frac{2\gamma ^2L\sqrt{ϵ}}{\pi \chi })^{\frac{1}{3}},$$
(36)
which differs from the pure ferroelectric case only by the factor $`ϵ`$. Notice that $`ϵ`$ is expected to be large in the Coulomb or fermi glass regime and we expect that the effect of introducing localized electrons in the system is to enlarge the domain size. Notice that for small enough $`L`$, we may enter a regime where $`L_D/a_o<<(r_T/a_o)^{\frac{5}{2}}`$. In this regime, we should use Eq. (A22) for $`Re\chi _{CG}(q,0)`$. In this regime, we find that the enhancement in domain size by the factor $`ϵ`$ is lost, and the analysis is similar to the case of metals except that the screening length $`q_{sc}^1`$ is given by

$`q_{sc}^2{\displaystyle \frac{4\pi ^4e^2g^2}{144}}a_o^3(k_BT).`$Notice that in contrast to metals, the screening length increases as temperature is lowered.

As temperature is lowered we find that for large enough $`L`$ the domain size is still given by Eq. (36) except that the value of $`ϵ`$ goes thought two more regimes. At temperature $`\mathrm{\Delta }>k_BT`$ but $`\omega _T>I(r_o)`$, we find

$`ϵϵ_o+{\displaystyle \frac{\pi g^2e^4a_o^4}{6}}({\displaystyle \frac{T_o}{T}})^2,`$whereas at $`\omega _T<I(r_o)`$, we find

$`ϵϵ_o+{\displaystyle \frac{n^2}{10\pi }}ln({\displaystyle \frac{T_o}{T}}).`$Notice that our theoretical prediction of domain size variation as a function of temperature can be tested experimentally. Note that for small enough $`L`$, the effect of screening will become stronger as in the high temperature case $`k_BT>\mathrm{\Delta }`$, and the enhancement factor $`ϵ`$ in domain wall size is lost.

It should also be kept in mind that the Coulomb glass has hysteretic low frequency response. One should therefore expect a slow and sluggish transition over a wide temperature range starting from about the $`T_c`$ to roughly when $`\omega _0^2+4\pi \chi 0`$.

### B (longitudinal) collective modes

At $`T>T_c`$, the longitudinal collective modes of the system can be obtained by solving

$$1v(q)\chi _{CG}\frac{4\pi \chi }{\omega ^2\omega _o(q)^2}=0.$$
(37)
In the $`\omega >>q0`$ limit, we obtain the secular equation

$$\omega ^2=\omega _o(q)^2+\frac{4\pi \chi }{1v(q)\chi _{CG}(q0,\omega )},$$
(38)
where $`\chi _{CG}(q0,\omega )`$ are given by Eqs.(A6a), (A6b), at temperature range $`k_BT<\mathrm{\Delta }`$ at $`\omega <I(r_o)`$ and $`\omega >I(r_o)`$, respectively. At temperatures $`k_BT>>\mathrm{\Delta }`$, $`\chi _{CG}(q0,0)`$ is given by Eq.(A9). Although the quantitative behaviours of the collective modes are different, logarithmic corrections appear in the pure dielectric collective mode frequency $`\omega _o(q)`$ in all the above cases. The logarithmic correction vanishes in the ferroelectric instability limit $`\omega \omega _o(0)0`$, in agreement with the static instability analysis which indicates that the ferroelectric instability occurs exactly when $`\omega _o(0)0`$. Notice, however that a collective mode solution with $`\omega ^2\omega _o(q)^2+\frac{4\pi \chi }{ϵ}`$ exists even up to the instability limit $`\omega _o(0)0`$. The polarization collective modes are in general much broadened by the Coulomb glass particle-hole excitation spectrum because of the very slow vanishing of $`Im\chi _{CG}(q,\omega )`$ at small $`\omega `$. Notice that at $`T<T_c`$, we may obtain the collective modes by replacing $`\omega _o^2`$ by $`2\omega _o^2`$ in Eq. (37). This is valid as long as we are interested at responses of the system at wave vectors $`q<L_D^1`$.

## A Electromagnetic response of Coulomb glass

In this appendix we summarize our theoretical understanding for response functions of Coulomb glass systems. (Notice that we have set $`e^2=1`$ in the main text. Here $`e^2`$ is kept to make the equations more clear). The electromagnetic responses of Coulomb glass systems are rather complicated, and can be divided into many different regimes, depending on temperature, frequency and wavevector. A detailed discussion of these can be found in ref\[\].

The starting model for responses of Coulomb glass systems is the effective two-site Hamiltonian by Shklovskii and Efros. The Hamiltonian is derived under the assumption that the system is in the strongly disordered regime where the electrons are strongly localized. In this regime tunneling events of electrons from one site to another is rare and can be treated in a dilute particle-hole pairs approximation where interaction between different tunneling events are neglected and the tunneling of electron between two sites (particle-hole pair) is described by a two-site Hamiltonian

$$H_{ij}=E_in_i+E_jn_j+\frac{e^2}{r_{ij}}n_in_j+I(r_{ij})(a_i^+a_j+a_j^+a_i),$$
(A2)
where $`n_i`$ is the electron occupation number ($`=0,1`$) at site $`i`$, $`r_{ij}`$ is the distance between the two sites, $`I(r)I_oe^{r/a_o}`$ is the effective tunneling matrix element between the two sites, where $`a_o`$ lattice spacing and $`a(a^+)`$ are electron annihilation and creation operators. the on-site energy $`E_i`$ is given by

$$E_i=ϵ_i+\underset{li,j}{}\frac{n_le^2}{|r_lr_i|},$$
(A3)
where $`A<ϵ_i<A`$ is the random on-site energy with bandwidth $`2A`$ and the second term comes from Coulomb interaction. Notice that electron spins are ignored in the theory. $`H_{ij}`$ can be diagonalized easily to obtain the following energy spectrum,

$`E_1^\pm `$ $`=`$ $`{\displaystyle \frac{1}{2}}(E_i+E_j)\pm {\displaystyle \frac{\mathrm{\Gamma }}{2}}`$ (A4)
$`E_2`$ $`=`$ $`E_i+E_j+{\displaystyle \frac{e^2}{r_{ij}}}`$ (A5)
where $`E_1^+`$ and $`E_1^{}`$ are the two possible energy states when the two-site system is occupied by a single electron, and $`E_2`$ is the energy when both sites are occupied, $`\mathrm{\Gamma }=\sqrt{(E_iE_j)^2+4I(r_{ij})^2}`$. Notice that for the pair of sites to be occupied by a single electron and participate in the electromagnetic response, we must have $`E_1^{}<\mu <E_2E_1^{}`$, where $`\mu `$ is the chemical potential. The whole system is described as a diluted gas of independent (i.e. non-interacting), randomly distributed particle-hole pairs and the total response of the system to external electromagnetic fields is given by by sum of the response of individual pairs. For example, the density-density response function of a pair in the two-site Hamiltonian approximation is found to be (at temperature $`T=0`$)

$$\chi _p(\stackrel{}{q},i\omega )=2\frac{I(r_{ij})^2}{\mathrm{\Gamma }(r_{ij})^2}(1cos(\stackrel{}{q}.\stackrel{}{r}_{ij}))[\frac{1}{i\omega \mathrm{\Gamma }(r_{ij})}\frac{1}{i\omega +\mathrm{\Gamma }(r_{ij})}],$$
(A7)
and the average density-density response function of the whole system is

$`\chi _{CG}(\stackrel{}{q},i\omega )`$ $`=`$ $`<{\displaystyle \underset{p}{}}\chi _p(\stackrel{}{q},i\omega )>`$ (A8)
$``$ $`2{\displaystyle _A^A}dE_1{\displaystyle _{E_1}^A}dE_2{\displaystyle }d^3r{\displaystyle \frac{I(r)^2}{\mathrm{\Gamma }(r)^2}}(1cos(\stackrel{}{q}.\stackrel{}{r}))F(E_1,E_2,r){\displaystyle \frac{2\mathrm{\Gamma }(r)}{(i\omega )^2\mathrm{\Gamma }(r)^2}},`$ (A9)
where $`F(E_1,E_2,r)`$ is the probability that a particle-hole pair with on-site energies $`E_1`$ and $`E_2`$ is found separated by distance $`r`$. For $`E_1,E_2`$ larger than the Coulomb gap energy $`\mathrm{\Delta }(e^2/a_o)^{3/2}/A^{1/2}`$, $`Fg^2`$ where $`g=n/Aa_o^3`$ is the density of states above the Coulomb gap, $`n`$ is the average number of electrons per site. For $`E_1,E_2`$ much smaller than $`\mathrm{\Delta }`$, $`F\rho (E_1)\rho (E_2)`$, where

$`\rho (E)={\displaystyle \frac{3n}{\pi e^6}}E^2.`$Notice that a two-site system behaves as a electric dipole and has only longitudinal response to external electric fields. Evaluating the integrals, we find for $`I(r_o)<<\omega <<A`$,

$$Im\chi _{CG}(q,\omega )a_o\pi g^2e^2r_\omega (1\frac{sin(qr_\omega )}{qr_\omega }),$$
(A11)
where $`r_\omega =a_oln(\frac{2I_o}{\omega })`$, $`r_o=e^2/\mathrm{\Delta }`$. For $`\omega <<I(r_o)`$, we obtain

$$Im\chi _{CG}(q,\omega )\frac{3n^2a_o}{10\pi e^2}\frac{1}{r_\omega ^3}(1\frac{sin(qr_\omega )}{qr_\omega }).$$
(A12)
The real part of the response function can be obtained using the Kramers-Kronig relation, we obtain for $`qr_o<<1`$,

$$Re\chi _{CG}(q,0)\frac{n^2}{20\pi ^2e^2}q^2ln(qr_o),$$
(A14)
and for $`qr_o>>1`$,

$$Re\chi _{CG}(q,0)\frac{g^2}{2}e^2r_o^2=\frac{g^2}{2}r_o^3\mathrm{\Delta }.$$
(A15)
For $`I(r_o)>\omega `$ and $`qr_\omega <<1`$, we obtain

$$Re\chi _{CG}(q,\omega )\frac{n^2}{20\pi ^2e^2}q^2ln(ln(\frac{\omega }{2I_o})),$$
(A17)
and for $`\omega >I(r_o)`$,

$$Re\chi _{CG}(q,\omega )\frac{g^2e^2}{24}q^2r_\omega ^4.$$
(A18)
At finite temperatures $`T<\mathrm{\Delta }`$, eq. (A12) for $`Im\chi _{CG}(0,\omega )`$ is cutoff at a low frequency $`\omega _Te^{(T_o/T)^{\frac{1}{2}}}`$ by the thermal mechanism of Mott-variable range hopping (in the presence of Coulomb gap) where $`T_oe^2/a_o`$. Correspondingly, at temperature $`\omega _T<I(r_o)`$ $`Re\chi _{CG}(q,0)`$ is replaced by

$$Re\chi _{CG}(q,0)\frac{n^2}{20\pi ^2e^2}q^2ln(\frac{a_o}{r_{\omega _T}})\frac{n^2}{40\pi ^2e^2}q^2ln(\frac{T}{T_o}),$$
(A19)
at $`\omega <\omega _T,(q=0)`$ and at $`qr_o<r_o/r_{\omega _T},(\omega =0)`$. At temperature $`\omega _T>I(r_o)`$, a corresponding analysis gives

$$Re\chi _{CG}(0,0)\frac{g^2e^2}{24}q^2r_{\omega _T}^4\frac{g^2e^2a_o^4}{24}q^2(\frac{T_o}{T})^2,$$
(A20)
for the same low frequency range and at small wavevector range $`qr_o<(r_o/r_{\omega _T})^2`$.

At still higher temperatures $`k_BT>>\mathrm{\Delta }`$ where the Coulomb gap becomes unimportant and phonon-assisted hopping becomes dominant, we enter another regime where $`\chi _{CG}`$ become different again. At frequencies $`qr_\omega <<1`$ and $`\omega >\overline{\omega }_T`$, where $`\overline{\omega }_T\omega _Pe^{(A/k_BT)^{\frac{1}{4}}}`$, where $`\omega _P`$ is a characteristic phonon frequency, we obtain

$$Re\chi _{CG}(q,\omega )\frac{\pi ^3g^2k_BT}{1440}q^2r_\omega ^5,$$
(A21)
and for $`\omega =0`$,

$$Re\chi _{CG}(q,0)\frac{\pi ^3g^2}{144}a_o^3(k_BT).$$
(A22)
For frequencies $`\omega <\overline{\omega }_T`$ and for wavevectors $`qa_o<<(a_o/r_T)^{\frac{5}{2}}`$, where $`r_T(\frac{A}{k_BT})^{\frac{1}{4}}a_o`$ is the characteristic thermal hopping distance, $`Re\chi _{CG}`$ saturates at a value

$$Re\chi _{CG}(q,0)\frac{\pi ^3g^2k_BT}{1440}q^2a_o^5(\frac{A}{k_BT})^{\frac{5}{4}}.$$
(A23)
Notice that the Mott variable range hopping introduce characteristic cutoff length scales $`r_{\omega _T})`$ and $`r_T`$ at different temperature regimes. As a result the dielectric constant $`ϵ(q,0)=1v(q)\chi _{CG}(q0,0)`$ becomes finite but large at low temperatures although the temperature dependence of its value goes through several different regimes.