text:
# Finite-size calculations of spin-lattice relaxation rates in Heisenberg spin-ladders

## I Introduction

Spin-ladder systems, in particular the two-leg, $`S=1/2`$, antiferromagnetic variety, have been the subject of considerable theoretical and experimental investigation.dagotto96  Spin ladders are appealing because they are one-dimensional systems and thus can be effectively investigated using many powerful theoretical tools, while offering a wider parameter space of “simple,” and potentially experimentally realizable, Heisenberg Hamiltonians than spin chains. The simplest Heisenberg ladder Hamiltonian has the form

$$=\underset{n}{}J_{}(𝑺_{n,1}𝑺_{n+1,1}+𝑺_{n,2}𝑺_{n+1,2})+J_{}𝑺_{n,1}𝑺_{n,2}$$
(1)
which offers a dimensionless parameter $`J_{}/J_{}`$ which is in principle tunable by chemistry or pressure. In addition, compounds containing weakly coupled Cu2O3 ladders are appealing because of possible connections with cuprate superconductivity.

The present work was motivated by the nuclear spin-lattice relaxation measurements in La6Ca8Cu24O41, an undoped ladder compound, by Imai et al.imai98  These measurements were carried out for all of the nuclear sites on the ladder, namely the copper, the “rung” oxygen, and the “ladder” (or “chain”) oxygen, over a wide temperature range, from low temperatures up to nearly $`900\mathrm{K}`$. Because the principal exchange interactions in cuprates are so large, on the order of $`1000\mathrm{K}`$, it is quite challenging to do experimental work at temperatures significantly greater than the spin gap ($`\mathrm{\Delta }500\mathrm{K}`$).

The experimental results (see Figure 1(c) of Imai et al.imai98 ) have the following noteworthy features. At temperatures below about $`425\mathrm{K}`$, the relaxation rates for all three sites follow a common (activated) temperature dependence up to a scale factor. However, on increasing $`T`$ the copper $`1/T_1`$ (which we will refer to as $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$) exhibits a rather sharp departure from that of the two oxygen sites ($`1/{}_{}{}^{\mathrm{O}(1)}T_{1}^{}`$ and $`1/{}_{}{}^{\mathrm{O}(2)}T_{1}^{}`$ for ladder and rung, respectively). There seems to be a nearly discontinuous decrease in the derivative of $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$; moreover, above $`425\mathrm{K}`$ the $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$ data appear nearly linear with an almost vanishing intercept. The relaxation rates for the two oxygen sites, in contrast, exhibit no particular features in the vicinity of $`425\mathrm{K}`$.

Several aspects of the wavevector dependence of the low-frequency spin susceptibility can be gleaned directly from the data.

One can express the spin-lattice relaxation rate in terms of the dynamic structure factor for the Cu2+ spins

$$\frac{1}{{}_{}{}^{n}T_{1}^{}}𝑑𝒒H_n(𝒒)S(𝒒,\omega _n)$$
(2)
where $`H_n`$ is the hyperfine form factor associated with nucleus $`n`$, $`\omega _n`$ is the NMR frequency (which we will take to be zero in everything that follows), and $`S`$ is the structure factor. The proportionality constants can be neglected for our present purposes. The spin correlations are isotropic, so there is no need to consider the various components, $`S^{xx}`$ and so forth, individually. The hyperfine interactions are not isotropic, so the orientation of the magnetic field in the NMR experiment does affect the results; however, all of the results of present interest can be obtained with a single field orientation, which then specifies $`H_n(𝒒)`$ uniquely. The largest hyperfine couplings are between a given nuclear site and the closest spins; at that level of approximation, and taking the intra- and inter-chain lattice constants to be of unit length, one has

$`H_{\mathrm{Cu}}=A^2,H_{\mathrm{O}(1)}=4C^2\mathrm{cos}^2(q_x/2),`$
$`H_{\mathrm{O}(2)}=4F^2\mathrm{cos}^2(q_y/2)+D^2`$ (3)
where $`C`$, $`F`$, and $`D`$ are the hyperfine couplings identified in Fig. 1(a) of Imai et al.,imai98  $`A`$ is the on-site hyperfine interaction for copper, and we have elided the orientation dependence of the hyperfine interactions (so, for example, $`A^2`$ should really be $`A_x^2+A_y^2`$ if the static field is along the $`z`$ axis).

The essential difference between copper and oxygen sites is that in the latter the hyperfine interaction in the vicinity of $`𝒒=(\pi ,\pi )`$ is much smaller than in the vicinity of $`𝒒=(0,0)`$. If, at all temperatures of experimental relevance, $`S(𝒒,0)`$ had most of its weight in the vicinity of $`𝒒=(0,0)`$, then then all three relaxation rates would have tracked one another. The marked decrease of $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$ relative to the other two relaxation rates at $`425\mathrm{K}`$ indicates that this cannot be the case, and in fact suggests that at temperatures below $`425\mathrm{K}`$ the ratio of the spectral weight near $`(\pi ,\pi )`$ to that near $`(0,0)`$ is roughly constant and of order unity, while above $`425\mathrm{K}`$ the ratio falls markedly. (The decrease is crucial. If there were an increase in $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$ relative to the oxygen rates with increasing $`T`$, one could ascribe that to a turn-on of $`S((\pi ,\pi ),0)`$ for $`T\mathrm{\Delta }`$ but $`S((\pi ,\pi ),0)`$ might have been negligible compared to $`S((0,0),0)`$ at lower temperatures.)

Why the emphasis on $`𝒒=(0,0)`$ and $`(\pi ,\pi )`$? In gapped systems such as spin ladders, the low-energy spin fluctuations are Raman processes, and at low temperatures one needs to consider only the lowest energy magnons, namely those near $`𝒒=(\pi ,\pi )`$. Spin fluctuations near $`(0,0)`$ are associated with two-magnon processes, and those near $`(\pi ,\pi )`$ with three-magnon processes, and on the face of it one would be justified in neglecting the three-magnon processes entirely at low temperatures: see Ref. ivanov99,  and references cited therein. However, as we have just seen, this appears to be inconsistent with the experimental data for La6Ca8Cu24O41, and it is also inconsistent with the quantum Monte Carlo calculations of spin-lattice relaxation in a particular Heisenberg ladder ($`J_{}/J_{}=1`$) by Sandvik, Dagotto, and Scalapino,sandvik96  at least at temperatures greater than half the magnon gap.

An extensive theoretical treatment of spin dynamics in gapped one-dimensional Heisenberg models, including spin ladders, has been presented by Damle and Sachdev.damle98  Their analysis of $`S(𝒒,\omega 0)`$ was restricted to $`𝒒`$ near $`(0,0)`$, but they did find the quite interesting result that the activation energy for $`1/T_1`$ is larger, by a factor of $`3/2`$, than the activation energy for the uniform static susceptibility (which is simply the spin gap). An analysis of $`S(𝒒,\omega 0)`$ for $`𝒒`$ near $`(\pi ,\pi )`$, for systems with $`J_{}J_{}`$ has been presented by Ivanov and Lee.ivanov99  Their results are suggestive of a fairly sharp crossover from low- to high-temperature regimes at $`T\mathrm{\Delta }`$, and also indicate that the $`(\pi ,\pi )`$ contribution to $`1/T_1`$ “overshoots” its $`T=\mathrm{}`$ value and thus decreases as $`T\mathrm{}`$.

In the present work, we have applied exact diagonalization to evaluate spin-lattice relaxation rates, following the method of Sokol, Gagliano and Bacci.sokol93  We have considered three different ladder Hamiltonians, namely $`J_{}/J_{}=0.5`$, 1.0, and 2.0, and have obtained $`1/T_1`$ for Cu, O(1), and O(2) sites taking the simplest conceivable hyperfine couplings, namely $`A=C=F=1`$, with all other interactions neglected. All of the calculations were for rather small systems, $`2\times 6`$, such that exact diagonalization could be carried out in an extremely straightforward manner.

It was noted above that calculations of spin lattice relaxation rates for spin ladders have already been carried out by means of large scale quantum Monte Carlo,sandvik96  but those calculations were limited to the Cu sites. Our goal is somewhat different than that of Sandvik, Dagotto, and Scalapino’s work. We are not trying to fit the data in detail, rather we want to see what can be learned from modest numerical calculations. One reason not to fit the data is that to get the gap correct to 10% by exact diagonalization for $`J_{}/J_{}=0.5`$ would require a system at least $`2\times 12`$. Another is that we do not treat the spin diffusion contribution to the relaxation rates correctly: our calculations effectively introduce an artificial cut-off so that we obtain a finite spin-lattice relaxation rate. Finally, the precise form of the spin Hamiltonian for the cuprate ladder compounds is still subject to argument. Although the Knight-shift results of Imai et al.imai98  appear to be consistent with the simple spin-ladder Hamiltonian of Eq. (1) for $`J_{}/J_{}0.5`$, it has been suggested by Brehmer et al.brehmer99  that instead $`J_{}/J_{}1`$ and in addition there is a modest amount of plaquette “ring exchange” in the Hamiltonian. A quantum-chemical analysis of the exchange interactions in various cupratesmizonu98  provides some support for the latter proposal, since it concludes that $`J_{}/J_{}0.9`$ for Sr14Cu24O41 (which is a lightly-self-doped version of the undoped La6Ca8Cu24O41 compound).

To be precise, the goals of our calculation are as follows. First, we want to verify that $`S(𝒒,0)`$ has significant weight near $`𝒒=(\pi ,\pi )`$ as well as near $`(0,0)`$ and see if there are any noticeable trends with varying $`J_{}/J_{}`$. Second, we want to explore the crossover from low to high temperature behavior in $`1/T_1`$: can we see anything like the experimental results, or like the theoretical results of Ivanov and Lee? Third, we want to keep our eyes open for any unanticipated patterns that might emerge in the numerical results.

## II Method of Calculation and Results

The finite-size calculations of spin-lattice relaxation rates are carried out following Sokol, Gagliano, and Bacci.sokol93  Rather than repeating their discussion of the method let us make a few remarks. We take $`J_{}`$ as the unit of energy.

The first step in the calculation is a complete diagonalization of the Hamiltonian and evaluation of matrix elements for certain local spin operators (depending on which nuclear site one is interested in). For the $`2\times 6`$ lattices all of the calculations could be done using the simplest possible representations of the states in terms of local $`S^z`$ values; it was not even necessary to use translational invariance to classify states by wave vector.

The second step is the construction of an auxiliary function which Sokol et al. refer to as $`I(\omega )`$. This is implicitly dependent on $`T`$ and the hyperfine couplings. We considered temperatures ranging from 0.3 to 50. Typically we constructed $`I(\omega )`$ at intervals of 0.02 in $`\omega `$ up to at least $`\omega =0.6`$.

Finally, one needs to estimate the zero-frequency derivative of $`I(\omega )`$, because $`1/T_1`$ is is proportional to $`T(I/\omega )|_{\omega =0}`$. At high temperatures $`I(\omega )`$ is quite smooth, but at temperatures comparable to the gap significant structure develops (see Fig. 1). In order to avoid introducing spurious temperature dependences into $`1/T_1`$ it is important to use a consistent procedure for extracting the derivative from the data. What we did was to fit a zero-intercept line through all the data points up to a cutoff $`\omega _{\mathrm{max}}`$, weighting all points equally in the fit. We did all of the calculations using both $`\omega _{\mathrm{max}}=0.5`$ and 0.3. While there are noticeable differences in the results using these two cutoffs, as shown in Fig. 2, our conclusions turn out the same no matter which is chosen. The use of a much smaller cutoff, which might seem to be preferred on the grounds that one is really looking for a zero-frequency derivative, is not beneficial. The structure that develops in $`I(\omega )`$ as $`T`$ is lowered, making it look like a Devil’s staircase, is a finite-size artifact and must be averaged over, using a suitably large $`\omega _{\mathrm{max}}`$, to obtain results that are representative of the thermodynamic limit.

We now turn to the results of the calculations for the three nuclear sites and three values of $`J_{}`$ considered (0.5, 1.0, and 2.0). In every case we take $`\omega _{\mathrm{max}}=0.5`$. In Fig. 3 we present results on a linear temperature scale, for $`T2`$. The behavior of the spin-lattice relaxation rate at high temperatures is a bit surprising: comparing the plots in Fig. 3(a) through (c) it is apparent that while the Cu and O(1) rates decrease strongly as $`J_{}`$ increases, the trend for the O(2) rate is different. This is made more explicit in Fig. 4, where we show $`1/T_1`$ for all three sites as a function of $`J_{}`$ at $`T=50`$ (effectively infinite temperature). In contrast, at low temperatures $`1/T_1`$ decreases with increasing $`J_{}`$ at all sites, as one would expect since the spin gap is an increasing function of $`J_{}`$.

## III Discussion and Conclusions

It is evident that for $`J_{}=0.5`$ and 1.0, $`1/T_1`$ for all three sites is nearly equal for temperatures below the spin gap. (Of course we do not claim that this holds to arbitrarily low temperatures, just that it seems correct for temperatures as low as we dare to estimate $`1/T_1`$.) Because of our choice of hyperfine interactions, this suggests that in such cases the weight in $`S(𝒒,0)`$ for $`𝒒(\pi ,\pi )`$ is approximately three times that for $`𝒒(0,0)`$. This is in quantitative agreement with the results of Sandvik et al.sandvik96  at $`J_{}=1.0`$. However, the story is rather different at $`J_{}=2.0`$, where the spin-lattice relaxation rates for all three sites, including the two oxygen sites, are significantly different even at $`T=\mathrm{\Delta }/2`$. In the strong-coupling limit, then, the simple picture for $`S(𝒒,0)`$ in which its weight is concentrated at $`(0,0)`$ and $`(\pi ,\pi )`$ does not work even for temperatures that are a modest fraction of $`\mathrm{\Delta }`$.

What can we say about the low-to-high temperature crossover in the spin-lattice relaxation rates? First of all, the sort of behavior seen experimentally, in which $`1/T_1`$ for the oxygen sites track each other closely while $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$ splits off, appears to be a special feature of $`J_{}1`$ in the present calculations; it is not at all generic and does not hold for the putative experimental value $`J_{}0.5`$. Second, in no case does $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$ exhibit any sort of sharp “break” as seen experimentally; nor does $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$ exhibit linear-in-$`T`$ behavior (with zero intercept, or otherwise) in the high temperature regime, even over a restricted temperature range (say $`\mathrm{\Delta }`$ to $`2\mathrm{\Delta }`$). Finally, in no case does $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$ exhibit an “overshoot” during the crossover: the spin-lattice relaxation rate associated with all sites monotonically increases with $`T`$.

Our calculations thus suggest that there are quite a few open problems in this field. Almost none of the prominent experimental facts concerning $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$ in La6Ca8Cu24O41 are reproduced in our finite-size calculations. Furthermore, the work of Ivanov and Leeivanov99  does not seem to have much to say about our results, either. Their calculation is controlled only in the $`J_{}1`$ regime, so we should only look at the $`J_{}=0.5`$ data. Here we have no evidence of overshoot in $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$, and no reason to believe that one can just examine the spectral weight near $`(\pi ,\pi )`$ since $`1/{}_{}{}^{\mathrm{O}(2)}T_{1}^{}`$ “peels off” from $`1/{}_{}{}^{\mathrm{O}(1)}T_{1}^{}`$ in a manner not very different from $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$.

At this point we face several alternatives. It is possible that our results are simply unreliable, because we are considering systems that are too small (especially for $`J_{}=0.5`$) and our procedure for estimating $`dI(\omega )/d\omega `$ is flawed. We cannot rule this out, but we strongly suspect that the trends in the results as a function of $`J_{}`$ are robust. It is possible that the spin Hamiltonian for the ladders in La6Ca8Cu24O41 is more complicated than the model we have considered. Whether the Hamiltonian of Brehmer et al.brehmer99  can reproduce the spin-lattice relaxation data requires another calculation. Another possibility that must be considered, given the remarkably sharp feature in $`1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}`$ found in the experimental data, is that La6Ca8Cu24O41 undergoes, by coincidence, a subtle structural transition at $`425\mathrm{K}`$. This could introduce an anomalously strong $`T`$-dependence to the hyperfine interactions, though why the effect should be so much stronger in $`H_{\mathrm{Cu}}(𝒒)`$ than $`H_{\mathrm{O}(1)}(𝒒)`$ and $`H_{\mathrm{O}(2)}(𝒒)`$ is difficult to envision.

Let us now turn to the results of our calculations for spin-lattice relaxation at very high temperatures, shown in Fig. 4. The most natural way to think about these results is in terms of the Gaussian approximation.anderson53 ; moriya56a ; moriya56b  The basic idea of this approach is to assume that $`𝑑𝒒H_n(𝒒)S(𝒒,\omega )`$ is a Gaussian function of $`\omega `$, and then evaluate the frequency cumulants of this function by means of short-time expansions of time-dependent correlation functions. At $`T=\mathrm{}`$ the calculations are especially simple, because the expectation values of correlators $`𝑺_i𝑺_j`$ vanish for sites $`ij`$. For three sites of interest in Heisenberg ladders, the Gaussian approximation yields the following exchange dependences of the spin-lattice relaxation rates at $`T=\mathrm{}`$:

$$1/{}_{}{}^{\mathrm{Cu}}T_{1}^{}1/\sqrt{1+\frac{1}{2}J_{}^2},$$
(4)
$$1/{}_{}{}^{\mathrm{O}(1)}T_{1}^{}1/\sqrt{1+J_{}^2},$$
(5)
and $`1/{}_{}{}^{\mathrm{O}(2)}T_{1}^{}`$ does not have any $`J_{}`$ dependence at all. (Recall that $`J_{}1`$; in all of these results there is an overall factor of $`1/J_{}`$.) If this last result seems peculiar, let us note that it can be derived in another way, by considering the strong-$`J_{}`$ limit. Then one most naturally thinks about the states in terms of singlets and triplets on the rungs. The relevant energy scale for the dynamics of the total spin on a rung, which is relevant to $`1/{}_{}{}^{\mathrm{O}(2)}T_{1}^{}`$, would seem to be proportional to $`J_{}`$ (that is, the bandwidth in lowest-order perturbation theory for a triplet excitation in a single backgroundbarnes93 ), and with the hypothesis of a single energy scale in $`𝑑𝒒H_{\mathrm{O}(2)}S(𝒒,\omega )`$ one reproduces the Gaussian approximation result.

We see in Fig. 4 that $`1/T_1`$ for the copper and ladder oxygen sites decreases with increasing $`J_{}`$, qualitatively in agreement with the Gaussian approximation, although the dependence on $`J_{}`$ is not as strong as that approximation suggests. Furthermore, $`1/{}_{}{}^{\mathrm{O}(2)}T_{1}^{}`$ exhibits an increase with $`J_{}`$. The rather poor performance of the Gaussian approximation is somewhat disappointing, considering how well it works for estimating spin-lattice relaxation rates in square-lattice Heisenberg antiferromagnets.gelfand93 ; sokol93  It is not too surprising, perhaps, given that the dynamic correlations in the $`S=1/2`$ Heisenberg chain are far from Gaussian at $`T=\mathrm{}`$.roldan86  So, there is yet another open problem in the area of low-energy spin dynamics of Heisenberg ladders.

###### Acknowledgements.

This work was supported by the US National Science Foundation through grant DMR 94–57928. We thank T. Imai for several stimulating discussions and also for communicating the results of his group’s experiments prior to publication.


format: text

text:
# 1 Introduction

## 1 Introduction

As yet, little is known about the structure of the inner regions of powerful radio galaxies and the impact of the activity on circumnuclear regions. The high resolution and sensitivity afforded by imaging observations using the Hubble Space Telescope (HST) has revealed a wealth of complex structures in such objects, with dust lanes, jets and large regions of scattered emission from hidden nuclei. One of the key targets in studies which aim to elucidate the interplay between active nuclei and their host galaxies, and between galactic activity of different types, is the archetypal powerful radio galaxy Cygnus A (see Carilli & Barthel 1996 for a review).

HST infrared imaging observations of Cygnus A by Tadhunter et al. (1999) revealed an edge-brightened bi-conical structure centred on the nuclear point source, which is strikingly similar to those observed around young stellar objects. The edge-brightening of this structure provides evidence that the bicone is defined as much by outflows in the nuclear regions as by the polar diagram of the illuminating quasar radiation field. The HST observations also show an unresolved nuclear source at 2.0 and 2.25$`\mu `$m. However, from the imaging observations alone it is unclear whether this unresolved source represents the highly extinguished quasar nucleus seen directly through the obscuring torus, or emission from a less-highly-extinguished extended region around the nucleus.

Near-IR polarization observations have the potential to remove the uncertainties surrounding the nature of the unresolved nuclear sources and, in addition, to provide further important information about the obscuration and anisotropy in the near-nuclear regions of powerful radio galaxies. Previous ground-based polarimetric observations of Cygnus A by Packham et al. (1998) demonstrate that the nuclear regions are highly polarized in the K-band, with a measured polarization of $`P_k4`$% for a 1 arcsecond diameter aperture centred on the compact IR nucleus. However, the resolution of the ground-based observations is insufficient to resolve the structures and determine the polarization mechanism unambiguously. In this letter we present new diffraction-limited infrared imaging polarimetry observations of Cygnus A made with the Near Infrared Imaging Camera and Multi-Object Spectrometer (NICMOS: MacKenty et al. 1997) on the HST. These observations resolve the polarized structures, and raise new questions about the nature of the anisotropy in the near-nuclear regions of this key source.

## 2 Observations and data reduction

NICMOS Camera 2 ‘long wavelength’ infrared imaging polarization observations were taken December 1997 and August 1998, giving a pixel scale of 0.075 arcseconds and a total field of 19.4$`\times `$19.4 arcseconds. The NICMOS polarizers are self-contained spectral elements, with the three long wavelength polarizers effective from $`\lambda 1.9`$$`2.1\mu `$m, resulting in an effective central wavelength of $`2.0\mu `$m (MacKenty et al. 1997). The polarizers are oriented at approximately 60 intervals and have characteristics as presented by Hines (1998). Each exposure consisted of a number of non-destructive reads of the detector which were optimally combined in the reduction software to remove cosmic rays. Regular chops were made to offset fields in order to facilitate accurate background subtraction. The total integration time was 2400 seconds per polarizer.

The reduction of the data used standard IRAF/STSDAS pipeline processing together with pedestal removal as given by van de Marel (1998). The three final, clean polarization images were combined following the prescription of Sparks & Axon (1999), using the pipeline produced variance data. The output data comprised a set of images containing each of the Stokes parameters $`I,Q,U`$, their variances and covariances, a debiassed estimate of the polarization intensity and polarization degree using the method of Serkowski (1958), and also position angle and uncertainty estimates on each of those images, as described in detail in Sparks & Axon (1999).

The two epochs of observation were acquired with different instrumental orientation, and were analysed independently to provide a robust check on our estimated uncertainties and on the possibility of systematic errors. A variety of spatial resolutions were used, from full diffraction-limited NICMOS resolution down to $`1`$ arcsec by smoothing the input three images prior to polarization analysis. Both epochs were fully consistent within the estimated statistical errors, implying that there are no significant sources of systematic uncertainty. Measurements of several Galactic stars in the field of Cygnus A give 2.0$`\mu `$m polarizations of $`P_{2.0\mu m}<2.0`$% for 5 pixel diameter apertures. This can be regarded as an upper limit on the level of systematic error in estimates of P. The typical statistical uncertainties for the most highly polarized regions measured in our full resolution images are $`\pm `$1.5% for $`P`$, and $`\pm `$3 for the polarization angle. In the following, we will present the data for the December 1997 observations only, since these have a higher S/N as a consequence of longer exposure times.

## 3 Results

Figure 1 shows the first epoch polarization results. As expected, the total intensity image (Stokes $`I`$) is very similar to the direct images published previously by Tadhunter et al. (1999). In particular, it shows an apparent edge-brightened, reasonably symmetric, bi-conical structure centred on the nucleus with an opening angle of 116 degrees and whose axis is closely aligned with the large scale radio jet.

The image of polarized intensity, however, reveals intriguing differences compared to total intensity. The only regions of strongly polarized emission (apart from the nucleus discussed below) are confined to a quasi-linear structure running along the NW-SE limb of the bi-conical structure. This feature shows approximate reflection symmetry about the nucleus, as opposed to the axial symmetry about the radio axis in the total intensity image. The brightness ratio of the two limbs of the cone to the east of the nucleus is approximately 2:1 in total intensity image, while in polarized intensity image it is $`>`$12:1. Note that the polarization structure visible in our 2.0$`\mu `$m image is strikingly different from that of the optical V- and B-band polarization images (Tadhunter et al. 1990, Ogle et al. 1997), in which the polarized emission appears uniformly distributed across the kpc-scale ionization cones and shows no clear preference for the NW-SE limb of the bi-cone.

Typical measured degrees of polarization are in excess of $`10`$% up to a maximum of $``$25% in the polarized region to the south east of the nucleus. However, these measures underestimate the true degree of intrinsic polarization in the extended structures, because starlight from the host galaxy makes a substantial contribution to the total flux. For example, using the azimuthal intensity profile measured in annulus with inner radius 4 pixels and outer radius 12 pixels, we estimate that the diffuse starlight from the host galaxy contributes 50 - 70% of the total flux in the south east arm of the bicone. Assuming that the starlight is unpolarized, the degree of intrinsic polarization in the south east arm is $`P_{2.0\mu m}^{intr}`$50 – 70%. Such high degrees of measured and intrinsic IR polarization are unprecedented in observations of active galaxies in which the synchrotron emitting jets are not observed directly at infrared wavelengths.

### 3.1 The nucleus

The nuclear point source, discussed in detail by Tadhunter et al. (1999), is also highly polarized. In the polarized intensity image the main nuclear component appears unresolved ($`FWHM=2.24`$ pixels) and its position agrees with that of the nucleus in the total intensity image to within 0.5 pixels (0.04 arcseconds). Thus, it appears likely that the bulk of the polarization is associated with the compact nucleus rather than a more extended region around the nucleus. The core does, however, show a faint entension to the NW in the polarized intensity image. This extension is aligned with the larger scale polarization structures, and its polarization E-vector is close to perpendicular to the radius vector from the nucleus.

From our full resolution polarization images, the measured degree of polarization at the peak flux of the nuclear point source is $`P_{2.0\mu m}^m20`$%. However, spurious polarization can arise because of small mis-alignments between the polarization images, especially in the nuclear regions where there are sharp gradients in the light distribution. To guard against such effects we have smoothed the polarization data using a 5$`\times `$5 pixel boxcar filter (0.375$`\times `$0.375 arcseconds) and re-measured the polarization in the nuclear regions. As expected, the measured degree of polarization in the nucleus in the smoothed image is less ($`P_{2.0\mu m}^m10`$%) than in the full resolution image, because of the greater degree of contamination by unpolarized starlight and extended structures around the nucleus. In order to determine the intrinsic polarization of the point source it is necessary to first determine the proportion of flux contributed by the point source to the total flux in the nuclear regions. Experiments involving the subtraction of a Tiny Tim generated point spread function (Krist & Hook 1997) suggest that an upper limit on the fractional contribution of the nucleus to the total flux in a 5$`\times `$5 pixel pixel box centred on the nucleus is $`f_{nuc}<35`$%. Thus, assuming that all the polarization in the near-nuclear regions is due to the unresolved compact core, and that the remainder of the light is unpolarized, the intrinsic polarization of the unresolved core source is $`P_{2.0\mu m}^{intr}=P_{2.0\mu m}^m/f_{nuc}>28`$%.

The position angle of polarization E-vector of the core source measured in our full resolution polarization map ($`PA=201\pm 3`$), is close to perpendicular to the radio jet axis ($`PA=105\pm 5`$). This is similar to the situation seen in other AGN, and in particular Cen A, where, towards longer wavelengths and smaller apertures, the infrared polarization becomes more and more closely perpendicular to the radio jet (Bailey et al. 1986).

## 4 Discussion

### 4.1 The nature of the unresolved core source

A major motivation for the HST observations was to investigate the nature of the compact core source and the cause of the relatively large polarization measured in the core by Packham et al. (1998). The explanation favoured by Packham et al. is that the compact core source represents transmitted quasar light, while the high polarization is due to dichroic absorption by aligned dust grains in the central obscuring torus. Because the dichroic mechanism is relatively inefficient, a high polarization implies a large extinction: from observations of Galactic stars it is known that at least 55 magnitudes of visual extinction is required to produce a K-band polarization of 28% for optimum grain alignment (Jones 1989). More typically, the correlation between K-band polarization and extinction deduced for Galactic stars by Jones (1989) implies that an extinction of $`A_v350`$ magnitudes would be required for $`P_k=28`$%. For comparison, an upper limit on the K-band extinction in Cygnus A, estimated by comparing the 2.25$`\mu `$m core flux with mid-IR and X-ray fluxes, is $`A_v<94`$ magnitudes (see Tadhunter et al. 1999 for details). Thus, the dichroic mechanism is only feasible if the efficiency of the mechanism in Cygnus A is greater than it is along most lines of sight in our Galaxy. Such enhanced efficiency cannot be entirely ruled out, given that the Galactic dichroic polarization involves a randomly oriented magnetic field component, whereas the magnetic fields in the central obscuring regions of AGN may be more coherent. In this context it is notable that near- and mid-IR polarization measurements of the central regions of the nearby Seyfert galaxy NGC1068 provide evidence for a greater dichroic efficiency than predicted by the Jones (1989) correlation, with $`P_k=5`$% produced by $`A_v`$20 – 40 magnitudes (Lumsden et al. 1999). However, even the greater dichroic efficiency deduced for NGC1068 would not be sufficient to produce the high polarization measured in the core of Cygnus A if $`A_v<94`$ magnitudes.

The efficiency problem might be resolved if the extinction to the core source in the K-band is higher than the $`A_v<94`$ estimated on the basis of comparisons of the K-band flux with the mid-IR and X-ray fluxes. Indeed, substantially higher extinctions have been deduced for Cygnus A, both from modelling the X-ray spectrum of the core ($`A_v=170\pm 30`$: Ueno et al. 1994) and from comparisons between hard-X-ray continuum, \[OIII\] emission line and mid-infrared continuum fluxes ($`A_v=143\pm 35`$: Ward 1996, Simpson 1995). If such high extinctions also apply to the quasar nucleus in the K-band, the low efficiency of the dichroic mechanism would be less of a problem. However, for any reasonable quasar SED, the contribution of such a highly obscured quasar nucleus to the flux and the polarization of the detected 2.0$`\mu `$m core source would be negligible (i.e. we would not expect to detect the quasar nucleus directly in the K-band). Thus, it is more likely that the relatively low extinction deduced from the K-band flux measurements reflects contamination of the K-band core by emission from a less-highly-obscured region, which is close enough to the central AGN to remain unresolved at the resolution of our HST observations. Although it has been proposed that the contaminating radiation in the K-band may include hot dust emission and/or line emission from quasar-illuminated regions close to the nucleus (e.g. Stockton & Ridgway 1996), such emission would have a low intrinsic polarization, and a large dichroic efficiency would still be required in order to produce the polarization of this component by dichroism.

The most plausible alternative to dichroic extinction is that the K-band core source represents scattered- rather than transmitted quasar light. In this case, the polarization is a consequence of scattering in an unresolved region close to the illuminating quasar; we do not detect the quasar nucleus directly in the K-band; and previous extinction estimates based on the K-band fluxes substantially underestimate the true nuclear extinction. Note that the presence of such a scattered component would resolve the discrepancy between the extinction estimates based on K-band flux measurements, and those based on fluxes measured at other wavelengths.

Finally, we must also consider the possibility that the core polarization is due to synchrotron radiation associated with the pc-scale jet visible in VLBI radio images (Krichbaum et al. 1996). Although the integrated polarization of the radio core is small even at high radio frequencies ($`P_{22GHz}<5`$%: Dreher 1979), we cannot entirely rule out the possibility that we are observing a highly polarized sub-component of the jet which suffers a relatively low extinction, or alternatively that the radio core source as a whole suffers large Faraday depolarization at radio wavelengths, and would appear more highly polarized at infrared wavelengths. Polarization observations at sub-mm wavelengths will be required to investigate this latter possibility.

### 4.2 The extended polarization structures

An intriguing feature of our HST observations is the high degree of polarization measured along, and only along, the NW-SE limb of the bicone. The orientation of the polarization measured along the limb is consistent with the scattering of light from a compact illuminating source in the nucleus, while the high degree of polarization is consistent with the edge-brightened bi-cone geometry of Tadhunter et al. (1999), in the sense that the scattering angle for the edge-brightened region will be close to the optimal 90 required for maximal polarization. However, the fact that the polarization is measured along only one limb is difficult to reconcile with the simplest bicone model in which the illuminating IR radiation field is azimuthally isotropic, and the scattering medium is uniformly distributed around the walls of the funnels hollowed out by the circum-nuclear outflows. In this simplest model both limbs would be highly polarized in the direction perpendicular to the radius vector of the source, and this is clearly inconsistent with the observations.

Our observation require that one or more of the assumptions implicit in the simple model must be relaxed. In general terms this requires either invoking specific matter distributions within the cone and/or an anisotropic illumination pattern of the central source itself.

Perhaps the simplest way of reconciling the polarization characteristics with the bi-cone geometry is to adjust the relative importance of scattering and intrinsic emmission with azimuth around the cone, so that one limb of the cone is dominated by scattered radiation, while the other is predominatly intrinsic radiation. Since the band-pass of the NICMOS polarizers contains the Paschen alpha line, this provides an obvious potential source of the diluting radiation for the unpolarized regions. Unfortunately there is no obvious reason why such an asymmetry should exist. Furthermore, direct images with the F222M filter show that there is no radical change in the relative brightness of the two limbs of the eastern cone between 2.0$`\mu `$m and 2.25$`\mu `$m. This is an argument against the Paschen alpha model for the intrinsic emission, since the F222M filter admits no emission lines as strong as Paschen alpha.

An alternative possibility is that the NW-SE limb is brighter because the near-IR radiation field of the illuminating AGN is more intense in that direction (i.e. the illuminating radiation field is azimuthally anisotropic within the cones). In this case, the clear difference in structure between the optical and near-IR reflection nebulae suggests that the sources of illumination at the two wavelengths are different: while the source of the shorter wavelength continuum must produce a radiation field which is azimuthally isotropic within cones defined by the obscuring torus, the source of illumination at the longer wavelengths is required to display considerable degree of azimuthal anisotropy. The near-IR continuum source must also have a relatively red spectrum, in order to avoid producing similar structures at optical and infrared wavelengths. The near-IR anisotropy might arise in the following ways.

1. Beamed radiation from the inner radio jet. The near-IR continuum is emitted by a component of the inner synchrotron jet which has a direction of bulk relativistic motion siginificantly displaced from the axis of the large-scale radio jet, such that the radiation is beamed towards the NW-SE limb of the bicone. However, given the remarkable degree of collimation, and the lack of bending, observed in the Cygnus A jet on scales betweem 1pc and 100kpc, a major difficulty with this model is that the the jet would have to bend through a large angle on a scale smaller than $``$1pc (the resolution of the VLBI maps), whilst retaining the rotational symmetry in the jet structure about the nucleus. A further requirement of this model is that, if the inner jet is precessing, the precession timescale must be greater than the light travel time ($``$5000 years) across the bicone structure.
2. Anisotropic hot dust emission. The near-IR continuum is emitted by hot dust in the inner regions of the galaxy, with a larger projected area of the emitting region visible from one limb of the bicone than from the other. For example, if the near-IR radiation is emitted by dust in a warped disk close to the central AGN — perhaps an outer part of the accretion disk — the warp could be oriented such that the NW-SE limb has an almost face-on view of the emitting region, whereas the the NE-SW limb has a more oblique view. There is already direct observational evidence for a warped outer accretion disk in at least one active galaxy (Miyoshi et al. 1995). Given that this mechanism would produce a relatively mild, broad-beam anisotropy, an optically thick torus on scales larger than the hot dust emitting region would still be required in order to produce the sharp-edges to the illuminated bicone structure.

Note that regardless of how any anisotropy in the IR radiation field might be produced, such anisotropy would not by itself explain the nature of the unpolarized emission along the SW-NE limb of the bicone, and the lack of variation in the brightness ratio of the two limbs of the eastern cone between 2.0$`\mu `$m and 2.25$`\mu `$m. It may also be difficult to reconcile the anisotropic illumination model with the polarization properties of the unresolved core source: if the core polarization is due to scattering, the orientation of the core polarization vector implies that a substantial flux of illuminating photons must escape at large angles to the NW-SE limb of the bicone.

We expect future spectropolarimetry observations to resolve the uncertainties concerning the origin of the near-IR polarization structures in Cygnus A. For example, in the case of the anisotropic illumination mechanisms considered above, the anisotropy is in the continuum flux rather than the broad lines associated with the AGN. Thus, if this model is correct, the broad lines will be relatively weak or absent in the polarized spectrum of the extended structures. In contrast, for a non-uniform distribution of scattering material but isotropic illumination within the cones, the broad lines and continuum will be scattered equally, and the equivalent widths of the broad lines in the polarized spectrum should fall within the the range measured for steep spectrum radio quasars.

## 5 Conclusions and Future Work

Our NICMOS polarimetry observations of Cygnus A have demonstrated the existence of a compact reflection nebula around the hidden core, but one whose polarization properties are inconsistent with the simplest illumination model suggested by the imaging data. The predominantly axial symmetry of the total intensity imaging is replaced by axial asymmetry and reflection symmetry about the nucleus in polarized light.

We have discussed several mechanisms to explain the near-IR polarization structures. While none of these is entirely satisfactory, it is clear that the near-IR polarization properties have the potential to provide key information about the geometries of the central emitting regions in AGN, and the near-IR continuum emission mechanism(s). In this context, it will be interesting in future to make similar observations of a large sample of powerful radio galaxies in order to determine whether the extraordinary IR polarization properties of Cygnus A are a common feature of the general population of such objects.

Acknowledgments. Based on Observations made with the ESA/NASA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. We thank the referee — Stuart Lumsden — for useful comments. A. Robinson acknowledges support from the Royal Society. References

Bailey, J.A., Sparks, W.B., Hough, J.H., Axon, D.J., 1986, Nat, 332, 150

Carilli C., Barthel, P.D., 1996, A&ARev, 7, 1

Dreher, J.W., 1979, ApJ, 230, 687

Hines, D.C., 1998, NICMOS & VLT, ESO Workshop and Conference Proceedings, 55, Wolfram Freudling and Richard Hook (eds), p63

Jones, T.J., 1989, ApJ, 346, 728

Krichbaum, T.P., Alef, W., Witzel, A., 1996, in Cygnus A — Study of a Radio Galaxy, ed. C.L. Carilli & D.E., Harris (Cambridge: Cambridge University Press), 93.

Krist J.E., Hook, R., 1997, TinyTim User Guide, Version 4.4 (Baltimore:STScI)

Lumsden, S.L., Moore, T.J.T., Smith, C., Fujiyoshi, T., Bland-Hawthorn, J., Ward, M.J., 1999, MNRAS, 303, 209

MacKenty J.W., et al. 1997, NICMOS Instrument Handbook, Version 2.0 (Baltimore: STScI)

Miyoshi, M., Moran, J., Herrnstein, J., Greenhill, L., Nakai, N., Diamond, P., Inoue, M., 1995, Nature, 373, 127

Ogle P.M., Cohen, M.H., Miller, J.S., Tran, H.D., Fosbury, R.A.E., Goodrich, R.W., 1997, ApJ, 482, L37

Packham, C., Hough, J.H., Young, S., Chrysostomou, A., Bailey, J.A., Axon, D.J., Ward, M.J., 1996, MNRAS, 278, 406

Packham, C., Young, S., Hough, J.H., Tadhunter, C.N., Axon., 1998, MNRAS, 297, 939

Serkowski, 1958, Acta.Astron., 8, 135

Sparks, W.B., Axon, D.J., 1999, PASP, in press

Stockton A., Ridgway S.E., Lilly, S., 1994, AJ, 108, 414

Stockton, A., Ridgway, S.E., 1996, in Cygnus A — Study of a Radio Galaxy, ed. C.L. Carilli & D.E., Harris (Cambridge: Cambridge University Press), 1

Tadhunter C.N., Scarrott S.M., Rolph C.D., 1990, MNRAS, 246, 163

Tadhunter C.N., Metz, S., Robinson, A., 1994, MNRAS, 268, 989

Tadhunter C.N., Packham, C., Axon, D.J., Jackson, N.J., Hough, J.H., Robinson, A., Young, S., Sparks, W., 1999, ApJ, 512, L91

Ueno S., Katsuji K., Minoru N., Yamauchi S., Ward M.J., 1994, ApJ, 431, L1

Ward M.J., Blanco P.R., Wilson A.S., Nishida M., 1991, ApJ, 382, 115

Ward M.J., 1996, in Cygnus A — Study of a Radio Galaxy, ed. C.L. Carilli & D.E., Harris (Cambridge: Cambridge University Press), 43

van der Marel, R., 1998: http://sol.stsci.edu/ marel/software.html

Figure 1. Infrared (2.0$`\mu `$m) polarization images of Cygnus A. Top left – total intensity (Stokes I) at full resolution; top right – polarization degree at full resolution; bottom left – polarized intensity at full resolution; and bottom right – polarization vectors plotted on a contour map of the polarized intensity image derived from the data smoothed with a 5$`\times `$5 pixel box filter, with length of vectors proportional to the percentage polarization. The line segment in the intensity image shows the direction of the radio axis. At the redshift of Cygnus A, 1.0 arcsecond corresponds to 1.0 kpc for $`H_0=75`$ km s<sup>-1</sup> Mpc<sup>-1</sup> and $`q_0=0.0`$.


text:
# The Detection of the Diffuse Interstellar Bands in Dusty Starburst Galaxies

## 1 Introduction

The Diffuse Interstellar Bands (“$`DIBs`$”) have been studied for over 60 years, since Merrill (1934) first established their origin in the interstellar medium. Despite decades of intensive investigation, the identity of the carrier or carriers of the $`DIBs`$ has not been established (see the comprehensive review by Herbig 1995). The most likely candidates are large carbon-rich molecules (e.g. Sonnentrucker et al 1997), perhaps Polycyclic Aromatic Hydrocarbons (PAH’s - Salama et al 1999). The strongest and best-studied $`DIBs`$ in the optical spectrum are known empirically to trace the $`HI`$ phase of the ISM, with strengths that correlate well with the line-of-sight color excess $`E(BV)`$, the $`HI`$ column density, and the $`NaI`$ column density as probed with the $`NaI\lambda \lambda `$5890,5896 (“$`NaD`$”) doublet (see Herbig 1993).

To date, $`DIBs`$ have been observed almost exclusively in our own Galaxy, and to a limited extent in the Magellanic Clouds (Morgan 1987). Supernova 1986G in NGC 5128 (Cen A) allowed the detection of $`DIBs`$ produced within the famous dusty gas disk in this elliptical galaxy (di Serego Alighieri & Ponz 1987). Most recently, Gallagher & Smith (1999) have reported the possible discovery of the $`DIB`$ at $`\lambda `$6283.9 Å in the spectra of two “super starclusters” near the nucleus of the prototypical starburst galaxy M 82. This is intriguing, since it suggests that the $`DIB`$ carriers are present at a normal level even in the ISM of an intense starburst, in which the ambient radiation intensity and gas pressure are orders-of-magnitude higher than in the diffuse gas in the Milky Way disk (e.g. Colbert et al 1999).

We have recently analyzed the properties of the interstellar $`NaD`$ absorption- line in the spectra of 18 high-luminosity, infrared-selected (dusty) starbursts (Heckman et al. 2000 - hereafter HLSA). In the course of this analysis, we examined the 7 starbursts with the highest quality spectra for the presence of $`DIBs`$ (section 2). As we report in section 3 below, we have detected one or both of the $`\lambda `$6283.9 Å and $`\lambda `$5780.5 Å $`DIB`$ features (normally the two strongest $`DIBs`$ in the optical spectral region) in all seven cases. We have also been able to map the spatial distribution of the $`DIBs`$. These data allow us to directly compare the properties of the $`DIBs`$ in the ISM of these extreme starbursts to sight-lines in the Galaxy having similar gas column density and reddening (section 4).

## 2 Observations & Data Analysis

Details concerning the following are given in HLSA, so we only summarize the most salient points here.

The starburst sample presented here is a subset of the 32 galaxies observed by HLSA. The HLSA sample itself was selected from two far-infrared-bright samples: the Armus, Heckman, & Miley (1989) sample of galaxies with very warm far-IR colors and the Lehnert & Heckman (1995) sample of far-IR-bright disk galaxies seen at high inclination. The combined sample is representative of the far-IR-galaxy phenomenon, but is not complete.

HLSA found that the $`NaD`$ line was of predominantly interstellar origin in 18 of the 32 galaxies, while cool stars contributed significantly to the line in the other 14 cases. The plethora of weak absorption features in the spectra of cool stars greatly complicate the detection of the $`DIBs`$, while the strength of $`DIBs`$ in our Galaxy correlate strongly with the ISM $`NaI`$ column density. Thus, the galaxies in the present paper were drawn exclusively from the 18 “interstellar-dominated” objects in HLSA. We then selected the objects in HLSA having the highest signal-to-noise spectra obtained with a resolution better than $``$ 100 km s-1 (see below). This results in a sample of 7 galaxies, as listed in Table 1.

The observations were undertaken in 1993 and 1994 using two different facilities: the 4-meter Blanco Telescope with the Cassegrain Spectrograph at $`CTIO`$ and the 4-meter Mayall Telescope with the RC Spectrograph at $`KPNO`$. The spectral resolution ranged from 1.1 Å FWHM in the $`KPNO`$ data to 1.8 Å FWHM in the $`CTIO`$ data. Details regarding spectrograph configurations are listed in Table 2 of HLSA.

The spectra were all processed using the standard LONGSLIT package in $`IRAF`$ (bias-subtracted, flat-fielded using spectra of a quartz-lamp, geometrically-rectified and wavelength-calibrated using a $`HeNeAr`$ arc lamp, and then sky-subtracted). See HLSA for details. No explicit correction was made for the presence of weak telluric absorption-features, but these are not a problem for our analysis. The strongest feature of relevance is the O2 band from $``$ 6276 to 6284 Å (e.g. Figure 2a in Benvenuti & Porceddu 1989). Fortunately, the redshifts of our galaxies are sufficient to move the $`\lambda `$6283.9 Å $`DIB`$ out from under this feature.

The spectra were analyzed using the interactive SPLOT spectral fitting package in $`IRAF`$. In all cases, a one-dimensional “nuclear” spectrum was extracted, covering a region with a size set by the slit width and summed over 5 pixels in the spatial direction (the resulting aperture is typically 2 by 4 arcsec). The corresponding linear size of the projected aperture is generally a few hundred parsecs to a few kpc in these galaxies (median diameter 600 pc). This is a reasonable match to the typical sizes of powerful starbursts like these (e.g. Meurer et al 1997; Lehnert & Heckman 1996). Prior to further analysis, each 1-D spectrum was normalized to unit intensity by fitting it with, and then dividing it by, a low-order polynomial. Similar one-dimensional spectra for off- nuclear regions were extracted over the spatial region with adequate signal- to-noise in the continuum for each galaxy.

It is essential to remove the myriad absorption features due to cool stars from the spectra before searching for the relatively weak $`DIB`$ features. We have therefore used the average spectrum of several Galactic K giant stars as a template. After redshifting the normalized stellar template to the galaxy rest-frame, we have iteratively scaled and subtracted the template from the normalized galaxy spectrum until the residuals in the difference spectrum were minimized in the spectral regions that exclude potentially detectable interstellar features. The scale factors found for the stellar template imply that cool stars typically contribute 20 to 30% of the continuum light at $``$ 6000Å. This is consistent with both the less rigorous estimates reported in HLSA for these galaxies, and with theoretical expectations for red supergiants in a mature metal-rich starburst (Bruzual & Charlot 1993; Leitherer et al 1999). To compensate for the effects of the continuum subtraction, we added back an equivalent amount of featureless continuum. Thus, the depths and equivalent widths of the $`DIBs`$ in the original data are preserved by our analysis. As an example, we show the spectrum of the nucleus of M82 before and after the subtraction of a suitably-scaled K-star spectrum in Figure 1. These final processed spectra are shown in Figures 2 and 3.

We have estimated the uncertainties in our measurements in two ways. First, we compared the measurements for the four galaxies in the sample for which we have more than one independent spectrum (taken at a different position angle). Second, we have calculated the rms noise in the cool-star-subtracted spectra and used this to calculate the implied uncertainties (assuming standard error propogation for Poissonian noise). We report these uncertainties in Tables 1 through 3.

## 3 Results

### 3.1 The Identified Features

The two most conspicuous $`DIBs`$ along typical sight-lines in the ISM of our Galaxy are the strong, relatively narrow features at $`\lambda `$6283.9 Å and $`\lambda `$5780.5 Å (Herbig 1995). In all but one case, both features are within the wavelength coverage of our spectra. The next strongest $`DIBs`$ in the Milky Way in the relevant spectral region are at 5797.0 Å, 6010.1 Å, 6203.1 Å, and 6613.6 Å. We have searched for all these features in our spectra.

We turn our attention first to the $`\lambda `$6283.9 Å feature, which is the strongest feature in the Milky Way, and is not seriously confused by stellar photospheric lines in the starburst spectra. As can be seen in Figure 2, the $`\lambda `$6283.9 $`DIB`$ is detected in 6 of the 7 starburst nuclei (Table 1). The only exception is NGC 6240, where the feature would lie within the blue shoulder of the very strong and broad \[OI\]$`\lambda `$6300 nebular emission-line, making it very difficult to detect. In the other six cases, the equivalent width of this $`DIB`$ ranges from $``$ 0.4 to 0.9 Å with a normalized residual intensity at line-center of 0.83 to 0.94. These values correspond to some of the strongest features seen along sight-lines in the ISM of the Milky Way (e.g. Chlewicki et al 1986; Benvenuti & Porceddu 1989).

Weaker absorption due to the $`\lambda `$5780.5 $`DIB`$ is definitely present in three of the seven members of our sample (NGC2146, M82, and NGC6240), and possibly present in three more (NGC1614, NGC1808, and NGC3256). No measurement can be made in IRAS10565+2448, since the feature lies just outside our spectral passband. In the five cases in which both features are detected, the $`\lambda `$5780.5 $`DIB`$ is typically about 25% as strong as the $`\lambda `$6283.9 feature, compared to a mean value of about 45% along comparably-reddened lines-of-sight in the Milky Way (Chlewicki et al 1986; Benvenuti & Porceddu 1989). We emphasize that the measurement of the $`\lambda `$5780.5 $`DIB`$ is difficult in our spectra owing to its proximity to the comparably strong stellar photospheric CrI+CuI$`\lambda `$5782 feature (with which it is badly blended). We estimate that this introduces an uncertainty of $`\pm `$ 50 mÅ in the quoted equivalent widths (which is generally larger than the formal measurement uncertainties estimated above).

The two starburst nuclei with the strongest $`\lambda `$6283.9 Å $`DIB`$ feature are NGC 2146 and M82. These two spectra also have the highest signal-to-noise and (along with IRAS10565+2448) have the best spectral resolution and broadest wavelength coverage in our sample. In these two spectra, several other weaker $`DIB`$ features can be identified, namely those at 5797.0 Å, 6010.1 Å, 6203.1 Å, and 6613.6 Å (Figure 3). The equivalent widths of these features are $``$100 mÅ or typically 10 to 15% as large as those of the $`\lambda `$6283.9 Å feature. These relative strengths agree reasonably well with Galactic $`DIBs`$ (Chlewicki et al 1986; Benvenuti & Porceddu 1989; Herbig 1995). We summarize this information in Table 2, and note that similarly-weak features could be present in the noisier spectra of the other five members of our sample.

### 3.2 Kinematics

We have measured the width and centroid of the $`\lambda `$6283.9 Å $`DIB`$ feature in all cases but NGC 6240 (where we have instead used the $`\lambda `$5780.5 Å $`DIB`$). The measured line widths (Table 3) range from $``$ 5 to 9 Å. The intrinsic width of the $`\lambda `$6283.9 ($`\lambda `$5780.5) $`DIB`$ in the Milky Way is $``$ 4 (2) Å (Herbig 1995). Taking our instrumental resolution into account, the implied Doppler broadening of the $`DIBs`$ due to macroscopic motions in the starburst ISM ranges from FWHM 160 to 430 km s-1. In four of the seven cases, these Doppler widths are smaller than the widths of the $`NaD`$ doublet (by 25 to 60%). In NGC1808, NGC2146 and M82, the $`DIB`$ $`NaD`$ lines have the roughly the same Doppler widths. Interestingly, HLSA find that these are the three cases in the present sample in which the nuclear $`NaD`$ lines do not show significant blueshifts with respect to the galaxy systemic velocity ($`v_{sys}`$).

The centroids of the $`DIBs`$ are within $``$ 100 km s-1 of $`v_{sys}`$. However, in all four cases with strongly blueshifted $`NaD`$ lines (NGC1614, NGC3256, IRAS10565+2448, and NGC6240), the $`DIBs`$ are mildly blueshifted (by $``$ 50 to 110 km s-1) with velocities that are intermediate between $`v_{NaD}`$ and $`v_{sys}`$. The velocities of the $`DIB`$ and $`NaD`$ absorbers roughly agree with one another (and lie close to $`v_{sys}`$) in the other three cases. This kinematic information is summarized in Table 3.

Taken together, these results suggest that the $`DIBs`$ trace gas that is more quiescent on-average than that probed by the $`NaD`$ line. That is, the $`NaD`$ absorption in the four “outflow” nuclei is probably produced by a combination of quiescent material ($`vv_{sys}`$ with smaller Doppler width) and disturbed, outflowing material. The bulk of the $`DIB`$ absorption would be associated with the former, and this component would dominate both the $`NaD`$ and $`DIB`$ absorption in the three other cases in our sample.

### 3.3 Spatial Extent

We have used our long-slit data to map the extra-nuclear spatial extent of the $`DIBs`$ in these galaxies. As listed in Table 1, these sizes range from $``$ 1 to 6 kpc. The absorbing region is larger (3 to 6 kpc) in the more powerful starbursts (NGC1614, NGC3256, IRAS10565+2448, and NGC6240, with $`logL_{bol}`$ = 11.3 to 12.0 $`L_{}`$), and smaller (0.9 to 1.8 kpc) in the less powerful cases (NGC1808, NGC2146, and M82, with $`logL_{bol}`$ = 10.5 to 10.7 $`L_{}`$). In the nearby (less powerful) starbursts, these sizes reflect the extent of the absorbing material. In the more distant (more powerful) starbursts, these sizes are lower limits set by the region with adequate signal-to-noise in the stellar continuum.

## 4 Discussion

### 4.1 Comparison to Galactic DIBs

The strengths of the prominent Galactic $`DIBs`$ correlate well with the column densities of both $`HI`$ and $`NaI`$ and with the reddening along the line-of-sight (e.g. Chlewicki et al 1986; Herbig 1993). This implies that the $`DIB`$ carrier is most plausibly associated with the cool atomic phase of the ISM. We can use the data discussed in HLSA to estimate the values for $`N_{NaI}`$ and $`E(BV)`$ in our sample of seven starbursts, to see if the $`DIBs`$ in our starburst sample obey the same empirical relations defined by the ISM of the Milky Way.

We follow HLSA and derive estimates for $`N_{NaI}`$ using the average of the values obtained from the classical “doublet ratio” method (Spitzer 1968) and the variant described by Hammann et al (1997). We estimate the line-of-sight reddening to the stellar continuum using the observed colors compared to theoretical models for a starburst stellar population (Leitherer et al 1999; see HLSA for details). We list the results in Table 1.

Our best-measured $`DIB`$ by-far is the strong $`\lambda `$6283.9 feature. The data compiled by Chlewicki et al (1986) and Benvenuti & Porceddu (1989) show that the mean ratio of the equivalent width of this feature and the color excess is $`<W_{6284}/E(BV)>`$ = 1.2 Å for heavily-reddened Galactic sight-lines. For our small starburst sample we find a similar result: $`<W_{6284}/E(BV)>`$ = 0.8 Å. This is shown in Figure 4 where we have plotted $`W_{6284}`$ vs. $`E(BV)`$ for a large sample of Galactic sight-lines using the extensive data compled by Herbig (1993). To compare our starburst data directly to this Galactic data we have converted the values given by Herbig (1993) for the equivalent width of the $`\lambda `$5780.5 Å $`DIB`$ into estimated values for $`W_{6284}`$ assuming that the mean ratio measured by Chlewicki et al (1986) and Benvenuti & Porceddu (1989) applies ($`W_{6284}/W_{5780}`$ = 2.2). In Figure 5 we have likewise plotted $`W_{6284}`$ vs. $`N_{NaI}`$ for both the Galactic data and our starburst data. The starbursts lie at the high-end of the relationship defined by $`DIBs`$ in the Milky Way.

### 4.2 Relationship to the $`\lambda `$2175 Å Dust Feature

Over the years, there has been considerable speculation as to a possible connection between the $`DIBs`$ and the strong and broad feature at $`\lambda `$2175 Å in the Galactic extinction curve (see Benvenuti & Porceddu 1989). In this context, the detection of strong $`DIBs`$ in starburst spectra is noteworthy. As shown by Calzetti et al (1994), the $`\lambda `$2175 feature is extremely (undetectably) weak in the UV spectra of starbursts. This implies that the carriers of the $`DIBs`$ and the $`\lambda `$2175 feature must be quite distinct (in agreement with the conclusions of Benvenuti & Porceddu (1989) for the Galactic ISM).

### 4.3 Speculations

On the face of it, the above results may seem surprising given the extreme differences between the physical conditions in the ISM of intense starbursts and our own Galactic disk. The strong starbursts in our sample have bolometric surface brightnesses of $`\mathrm{\Sigma }_{bol}10^{10}`$ to $`10^{11}`$ L kpc-2 (e.g. Meurer et al. 1997), typical star-formation rates per unit area of $`\mathrm{\Sigma }_{SFR}`$ 10 M year-1 kpc-2, and surface mass densities in gas and stars of $`\mathrm{\Sigma }_{gas}\mathrm{\Sigma }_{stars}`$ 109 M kpc-2 (e.g. Kennicutt 1998). These are roughly 103 ($`\mathrm{\Sigma }_{SFR}`$), 102 ($`\mathrm{\Sigma }_{gas}`$) and 101 ($`\mathrm{\Sigma }_{stars}`$) times larger than the corresponding values in the disks of normal galaxies. These values for $`\mathrm{\Sigma }_{bol}`$ correspond to a radiant energy density inside the star-forming region that is roughly 103 times the value in the ISM of the Milky Way (and see Colbert et al 1999 for direct measurements of this quantity). The rate of mechanical energy deposition (supernova heating) per unit volume in these starbursts is of-order 103 times higher than in the ISM of our Galaxy (e.g. Heckman, Armus, & Miley 1990), as is the cosmic ray heating rate (Suchkov, Allen, & Heckman 1993). Finally, simple considerations of hydrostatic equilibrium imply correspondingly high pressures in the ISM: $`PG\mathrm{\Sigma }_g\mathrm{\Sigma }_{tot}`$ few $`\times `$ 10-9 dyne cm-2 (P/k $``$ few $`\times `$ 107 K cm-3, or several thousand times the value in the local ISM in the Milky Way). These high pressures have been confirmed observationally (e.g. Heckman, Armus, & Miley 1990; Colbert et al 1999).

The interesting result of the above is that despite the extreme conditions prevailing inside these starbursts, the dimensionless ratio of the ISM pressure to the energy density in UV photons (or cosmic rays) is quite similar in starbursts and the disk of the Milky Way. This would in turn imply that (for a given ISM temperature) the ratio of the number densities of the gas particles and UV photons (or cosmic rays) would also be similar to their values in the local ISM. Wang, Heckman, & Lehnert (1998) have discussed the evidence that this analysis is correct for the diffuse ionized medium in starbursts and the disks of normal late-type galaxies.

This “homologous” behavior of the ISM in regions spanning over three orders-of-magnitude in heating and cooling rates per particle may help to explain why the ratio of the column density of $`DIB`$ carriers to that of both $`Na`$ atoms (Figure 5) and dust grains (Figure 4) appears so similar in extreme starbursts and the ISM of our own Galaxy. In the absence of a well-understood origin for the $`DIBs`$, further speculation seems premature.

## 5 Summary

Despite over six decades of investigation, the nature and origin of the Diffuse Interstellar Bands remain a mystery (Herbig 1995). We have presented evidence that - far from being a possibly pathological property of the local ISM in our Galaxy - $`DIBs`$ are probably ubiquitous in the spectra of far-infrared-bright (dusty) starbursts.

In our own Galaxy, the two most conspicuous $`DIBs`$ are the features at $`\lambda `$6283.9 Å and $`\lambda `$5780.5 Å. We have detected one or both of these two $`DIBs`$ in all seven starbursts selected on the basis of strong interstellar $`NaI\lambda \lambda `$5790,5796 ($`NaD`$) absorption from the larger starburst sample studied by Heckman et al (2000 - HLSA). The equivalent widths of these features are $``$ 400 to 900 mÅ and $``$ 100 to 400 mÅ for the $`\lambda `$6280.9 and $`\lambda `$5780.5 features respectively. These roughly correspond to the greatest $`DIB`$ strengths observed in the Milky Way (Herbig 1993; Chlewicki et al 1986). In two members of our sample (M82 and NGC2146) the spectra are of high enough signal-to-noise to detect four other weaker $`DIBs`$ (at 5797.0 Å, 6010.1 Å, 6203.1 Å, and 6613.6 Å). These have typical equivalent widths of $``$ 100 mÅ. The relative strengths of these $`DIBs`$ are rather similar to those in the Milky Way (Herbig 1995; Chlewicki et al 1986; Benvenuti & Porceddu 1989).

The $`DIBs`$ can be mapped over an extensive region in and around the nuclear starbursts. In the moderately powerful starbursts ($`L_{bol}`$ = few $`\times `$ 1010 L), this region is $``$ 1 kpc in size vs. several kpc in the more powerful starbursts ($`L_{bol}`$ = few $`\times `$ 1011 L). The kinematics of the gas producing the $`DIBs`$ is evidently more quiescent than that producing the $`NaD`$ absorption studied by HLSA. In the four starbursts with broad and strongly blueshifted $`NaD`$ lines, the $`DIBs`$ are less Doppler-broadened and much less blueshifted ($`v_{DIB}`$ \- $`v_{sys}`$ $``$ -100 km s-1).

In the Milky Way, the $`DIBs`$ are known to trace a dusty atomic phase of the ISM, since their equivalent widths correlate strongly with the $`HI`$ column density, the $`NaI`$ column density, and the reddening parameter $`E(BV)`$ (Herbig 1995 and references therein). We show that these starburst $`DIBs`$ obey the same trends with $`N_{NaI}`$ and $`E(BV)`$ (e.g. $`W_{6284}`$ 1.2 $`E(BV)`$ Å at log$`N_{NaI}`$ 14 cm-2). Thus, the abundance of the $`DIB`$ carrier(s) relative to $`Na`$ atoms and dust grains appears to be very similar in intense starbursts and the diffuse ISM of our own Galaxy.

This seems surprising, given the thousand-fold greater energy density in photons and cosmic rays in the ISM of an intense starburst (e.g. Colbert et al 1999; Suchkov, Allen, & Heckman 1993). However, the gas pressures and densities in the starburst ISM are correspondingly larger as well (e.g. Heckman, Armus, & Miley 1990). Thus, such key dimensionless ratios as gas/photon density and gas-pressure/radiant-energy-density are similar in the ISM of starbursts and the disks of normal spiral galaxies (Wang, Heckman, & Lehnert 1998). This apparent “homology” may help explain the strikingly similar $`DIB`$ properties.

Finally, we point out that starbursts apparently produce strong $`DIBs`$ without producing a detectable $`\lambda `$2175 Å dust feature in their UV spectra (Calzetti et al 1994). This underscores the quite distinct origin of the two types of features.

We thank David Neufeld, Ken Sembach, and Don York for useful conversations at various stages of this project. The partial support of this project by NASA grant NAGW-3138 is acknowledged.

Note. Col. (2) The equivalent width of the $`\lambda `$6283.9 $`DIB`$ in mÅ. Col. (3) The equivalent width of the $`\lambda `$5780.5 $`DIB`$ in mÅ. The uncertainty is due primarily to the accuracy with which contamination by the stellar photospheric CrI+CuI$`\lambda `$5782 can be removed. We estimate this leads to an uncertainty of $`\pm `$50mÅ. The detection of this $`DIB`$ is therefore only tentative in NGC1614, NGC1808, and NGC3256 (indicated by a colon). Col. (4) The angular size (in arcsec) over which the $`\lambda `$6283.9 $`DIB`$ is detectable (the $`\lambda `$5780.5 $`DIB`$ was used in NGC6240). Col. (5) The corresponding physical size (in kpc), for our adopted $`H_0`$ = 70 km s-1 Mpc-1. Col. (6) The estimated color excess along the line-of-sight to the stellar continuum, based on the observed continuum color and a model starburst spectral energy distribution (Leitherer et al 1999; see HLSA for details). Col. (7) The logarithm of the estimated column density of $`NaI`$ atoms (cm-2). These were derived using the standard doublet ratio technique (Spitzer 1968) and its variant in Hammann et al (1997). See HLSA for details. Based on an intercomparison of the values obtained by different techniques, we estimate the uncertainty to be $`\pm `$0.2 dex.

Note. Col. (2) The Doppler broadening (full-width-at-half maximum) in km s-1 for the $`\lambda `$6283.9 $`DIB`$ (the $`\lambda `$5780.5 $`DIB`$ was used in NGC6240). These widths have been corrected for the intrinsic width of the DIB feature (see text) and for the instrumental resolution of the spectrograph (see HLSA). The raw, measured line widths and their associated uncertainties (in Å) are given in parantheses. Col. (3) The full-width-at-half-maximum in km s-1 of the members of the $`NaI\lambda \lambda `$5890,5896 doublet ($`NaD`$). Uncertainties are $`\pm `$20 km s-1. Taken from HLSA. Col. (4) The heliocentric galaxy systemic velocity. Approximate uncertainties range from $`\pm `$10 km s-1 for NGC1808, NGC2146, and M82 to $`\pm `$50 km s-1 for NGC1614 and NGC3256, to $`\pm `$100 km s-1 for NGC6240 and IRAS10565+2448. See HLSA and references therein. Col. (5) The heliocentric velocity of the $`\lambda `$6283.9 $`DIB`$ (the $`\lambda `$5780.5 $`DIB`$ was used in NGC6240). The measurement uncertainties are $`\pm `$30 km s-1 for NGC2146 and M82, $`\pm `$50 km s-1 for NGC1614, NGC1808, and NGC3256, and $`\pm `$80 km s-1 for IRAS10565+2448 and NGC6240. These do not include any uncertainties in the true value of the rest wavelength for the $`DIB`$. Col. (6) The heliocentric velocity of the $`NaD`$ doublet taken from HLSA. Uncertainties are $`\pm `$20 km s-1.


text:
# Study of Multi-muon Events from EAS with the L3 Detector at Shallow Depth Underground

## 1 THE L3+COSMICS EXPERIMENT

The muon component of extensive air showers (EAS), due to the long muon range in the Earth’s atmosphere, carries a wealth of information about the shower development. Study of multi-muon events gives an insight into the primary cosmic ray composition and the physics of high energy hadronic interactions. The L3 detector, situated 30 m underground, offers interesting possibilities to detect and study such events , which are complementary to the data collected in traditional cosmic ray experiments. The hadron component of EAS is absorbed, while the muon component is detected with low threshold (typically, if we exclude access shafts, 15 GeV) and high momentum and spatial resolution by the sophisticated tracking system of the L3 detector. The muon spectrum can be measured up to 2 TeV with high precision. The multi-muon event rate is high enough to make studies of the knee region possible with one year of data taking.

This year, 5 billion triggers were collected with the full L3+Cosmics setup. The independent readout and data acquisition system allows us to take data in parallel with L3. The acceptance of the setup is 200 $`\mathrm{m}^2\mathrm{sr}`$. The angular resolution is better than 3.5 mrad for muons above 100 GeV and zenith angles from 0 to $`50^{}`$. The momentum resolution is 5.0 % at 45 GeV. It is calibrated with $`\mathrm{Z}\mu ^+\mu ^{}`$ events from the LEP calibration runs, where the muon momentum is known exactly.

## 2 MONTE CARLO SIMULATIONS

The simulation program ARROW  is used to calculate the hadron, muon and neutrino flux at the detector level. The method combines simulations for fixed energies and different primary nuclei with a parametrization of the energy dependence and allows to do fast calculations for different geometries and energy thresholds. Results for the L3+C setup with sensitive area 200 $`\mathrm{m}^2\mathrm{sr}`$

are shown in Figure 1.

The primary composition is divided in heavy (Fe) and light (p) components in two limiting hypotheses: Fe S for constant heavy contribution $`30`$% and Fe H for heavy contribution rising from 30% below to 70% above the knee. The events with up to 6 muons are dominated by proton induced showers and above 10 muons the iron takes over. To distinguish between the two hypotheses with this method we need to detect events with $`50`$ muons and more.

## 3 DATA AND OUTLOOK

A first small subset of our data is shown in Figure 2. Only part of the events with up to 6 muons are reconstructed currently. The observed charge ratio $`\mu ^+/\mu ^{}`$ from the raw data is flat in the region between 50 and 500 GeV.

In the year 2000 an EAS array will be mounted above L3+C in order to detect the primary energy and core position. The experimental program includes studies of muon families and the primary composition, sidereal anisotropies, high multiplicity events in coincidences with other experiments, the moon shadow, searches for point sources, gamma ray bursts and exotic events.


text:
# Rashba spin splitting in two-dimensional electron and hole systems

## Abstract

In two-dimensional (2D) hole systems the inversion asymmetry induced spin splitting differs remarkably from its familiar counterpart in the conduction band. While the so-called Rashba spin splitting of electron states increases linearly with in-plane wave vector $`k_{}`$ the spin splitting of heavy hole states can be of third order in $`k_{}`$ so that spin splitting becomes negligible in the limit of small 2D hole densities. We discuss consequences of this behavior in the context of recent arguments on the origin of the metal-insulator transition observed in 2D systems.

At zero magnetic field $`B`$ spin splitting in quasi two-dimensional (2D) semiconductor quantum wells (QW’s) can be a consequence of the bulk inversion asymmetry (BIA) of the underlying crystal (e.g. a zinc blende structure) and of the structure inversion asymmetry (SIA) of the confinement potential. This $`B=0`$ spin splitting  is the subject of considerable interest because it concerns details of energy band structure that are important in both fundamental research and electronic device applications (Refs.  and references therein).

Here we want to focus on the SIA spin splitting which is usually the dominant part of $`B=0`$ spin splitting in 2D systems.  To lowest order in $`k_{}`$ SIA spin splitting in 2D electron systems is given by the so-called Rashba model,  which predicts a spin splitting linear in $`k_{}`$. For small in-plane wave vector $`k_{}`$ this is in good agreement with more accurate numerical computations.  For 2D hole systems, on the other hand, the situation is more complicated because of the fourfold degeneracy of the topmost valence band $`\mathrm{\Gamma }_8^v`$, and so far only numerical computations on hole spin splitting have been performed.  In the present paper we will develop an analytical model for the SIA spin splitting of 2D hole systems. We will show that in contrast to the familiar Rashba model the spin splitting of heavy hole (HH) states is basically proportional to $`k_{}^3`$. This result was already implicitly contained in several numerical computations.  But a clear analytical framework was missing. We will discuss consequences of this behavior in the context of recent arguments on the origin of the metal-insulator transition observed in 2D systems.

First we want to review the major properties of the Rashba model

$$H_{6c}^{\mathrm{SO}}=\alpha 𝐤\times 𝐄𝝈.$$
(1)
In this equation $`𝝈=(\sigma _x,\sigma _y,\sigma _z)`$ denotes the Pauli spin matrices, $`\alpha `$ is a material-specific prefactor,  and $`𝐄`$ is an effective electric field that results from the built-in or external potential $`V`$ as well as from the position dependent valence band edge. For $`𝐄=(0,0,E_z)`$ Eq. (1) becomes (using explicit matrix notation)

$$H_{6c}^{\mathrm{SO}}=\alpha E_z\left(\begin{array}{cc}0& k_{}\\ k_+& 0\end{array}\right)$$
(2)
with $`k_\pm =k_x\pm ik_y`$. By means of perturbation theory we obtain for the spin splitting of the energy dispersion

$$_{6c}^{\mathrm{SO}}(𝐤_{})=\pm \alpha E_zk_{}$$
(3)
where $`𝐤_{}=(k_x,k_y,0)`$. Using this simple formula several groups determined the prefactor $`\alpha E_z`$ by analyzing Shubnikov-de Haas (SdH) oscillations.

Equation (3) predicts an SIA spin splitting which is linear in $`k_{}`$. For small $`k_{}`$ Eq. (3) thus becomes the dominant term in the energy dispersion $`_\pm (𝐤_{})`$, i.e., SIA spin splitting of electron states is most important for small 2D densities. In particular, we get a divergent van Hove singularity of the density-of-states (DOS) at the bottom of the subband  which is characteristic for a $`k`$ linear spin splitting. As an example, we show in Fig. 1 the self-consistently calculated  subband dispersion $`_\pm (k_{})`$, DOS effective mass $`m^{}/m_0`$, and spin splitting $`_+(k_{})_{}(k_{})`$ for an MOS inversion layer on InSb. For small $`k_{}`$ the spin splitting increases linearly as a function of $`k_{}`$, in agreement with Eq. (3). Due to nonparabolicity the spin splitting for larger $`k_{}`$ converges toward a constant.

The spin splitting results in unequal populations $`N_\pm `$ of the two branches $`_\pm (k_{})`$. For a given total density $`N_s=N_++N_{}`$ and a subband dispersion $`_\pm (k_{})=\mu k_{}^2\pm \alpha E_zk_{}`$ with $`\mu =\mathrm{}^2/2m^{}`$ we obtain

$$N_\pm =\frac{1}{2}N_s\pm \frac{\alpha E_z}{8\pi \mu ^2}\sqrt{8\pi \mu ^2N_s\alpha E_z^2}.$$
(4)
This equation can be directly compared with, e.g., the results of SdH experiments.

The Rashba model (1) can be derived by purely group-theoretical means. The electron states in the lowest conduction band are $`s`$ like (orbital angular momentum $`l=0`$). With spin-orbit (SO) interaction we have total angular momentum $`j=1/2`$. Both $`𝐤`$ and $`𝐄`$ are polar vectors and $`𝐤\times 𝐄`$ is an axial vector (transforming according to the irreducible representation $`\mathrm{\Gamma }_4`$ of $`T_d`$).  Likewise, the spin matrices $`\sigma _x`$, $`\sigma _y`$, and $`\sigma _z`$ form an axial vector $`𝝈`$. The dot product (1) of $`𝐤\times 𝐄`$ and $`𝝈`$ therefore transforms according to the identity representation $`\mathrm{\Gamma }_1`$, in accordance with the theory of invariants of Bir and Pikus.  In the $`\mathrm{\Gamma }_6^c`$ conduction band the scalar triple product (1) is the only term of first order in $`𝐤`$ and $`𝐄`$ that is compatible with the symmetry of the band.

Now we want to compare the Rashba model (1) with the SIA spin splitting of hole states. The topmost valence band is $`p`$ like ($`l=1`$). With SO interaction we have $`j=3/2`$ for the HH/LH states ($`\mathrm{\Gamma }_8^v`$) and $`j=1/2`$ for the SO states ($`\mathrm{\Gamma }_7^v`$). For the $`\mathrm{\Gamma }_8^v`$ valence band there are two sets of matrices which transform like an axial vector, namely $`𝐉=(J_x,J_y,J_z)`$ and $`𝓙=(J_x^3,J_y^3,J_z^3)`$ (Refs. ). Here $`J_x`$, $`J_y`$ and $`J_z`$ are the angular momentum matrices for $`j=3/2`$. Thus we get

$$H_{8v}^{\mathrm{SO}}=\beta _1𝐤\times 𝐄𝐉+\beta _2𝐤\times 𝐄𝓙.$$
(5)
Similar to the Rashba model the first term has axial symmetry with the symmetry axis being the direction of the electric field $`𝐄`$. The second term is anisotropic, i.e., it depends on both the crystallographic orientation of $`𝐄`$ and $`𝐤`$. Using $`𝐤𝐩`$ theory we find that the prefactor $`\beta _2`$ is always much smaller than $`\beta _1`$, i.e., the dominant term in Eq. (5) is the first term. This can be easily understood by noting that the $`𝐤𝐩`$ coupling between $`\mathrm{\Gamma }_8^v`$ and $`\mathrm{\Gamma }_6^c`$ is isotropic, so that it contributes to $`\beta _1`$ but not to $`\beta _2`$. The prefactor $`\beta _2`$ stems from $`𝐤𝐩`$ coupling to more remote bands such as the $`p`$ antibonding conduction bands $`\mathrm{\Gamma }_8^c`$ and $`\mathrm{\Gamma }_7^c`$.

For $`𝐄=(0,0,E_z)`$ Eq. (5) becomes (using explicit matrix notation with $`j=3/2`$ eigenstates in the order $`j_z=+3/2,+1/2,1/2,3/2`$)

$$H_{8v}^{\mathrm{SO}}=\beta _1E_z\left(\begin{array}{cccc}0& \frac{1}{2}\sqrt{3}k_{}& 0& 0\\ \frac{1}{2}\sqrt{3}k_+& 0& k_{}& 0\\ 0& k_+& 0& \frac{1}{2}\sqrt{3}k_{}\\ 0& 0& \frac{1}{2}\sqrt{3}k_+& 0\end{array}\right)+\beta _2E_z\left(\begin{array}{cccc}0& \frac{7}{8}\sqrt{3}k_{}& 0& 3/4k_+\\ \frac{7}{8}\sqrt{3}k_+& 0& 5/2k_{}& 0\\ 0& 5/2k_+& 0& \frac{7}{8}\sqrt{3}k_{}\\ 3/4k_{}& 0& \frac{7}{8}\sqrt{3}k_+& 0\end{array}\right).$$
(6)
Here the first term couples the two LH states ($`j_z=\pm 1/2`$), and it couples the HH states ($`j_z=\pm 3/2`$) with the LH states. But there is no $`k`$ linear splitting of the HH states proportional to $`\beta _1`$. The second matrix in Eq. (6) contains a $`k`$ linear coupling of the HH states.

We want to emphasize that $`H_{6c}^{\mathrm{SO}}`$ and $`H_{8v}^{\mathrm{SO}}`$ are effective Hamiltonians for the spin splitting of electron and hole subbands, which are implicitly contained in the full multiband Hamiltonian for the subband problem

$$H=H_{𝐤𝐩}(𝐤_{},k_z=i_z)+eE_zz𝟙.$$
(7)
Here $`H_{𝐤𝐩}`$ is a $`𝐤𝐩`$ Hamiltonian for the bulk band structure (i.e., $`H_{𝐤𝐩}`$ does not contain $`H_{6c}^{\mathrm{SO}}`$ or $`H_{8v}^{\mathrm{SO}}`$) and we have restricted ourselves to the lowest order term in a Taylor expansion of the confining potential $`V(z)=V_0+eE_zz+𝒪(z^2)`$ which reflects the inversion asymmetry of $`V(z)`$. The effective Hamiltonians (2) and (6) stem from the combined effect of $`H_{𝐤𝐩}`$ and the term $`eE_zz`$. For a systematic investigation of the importance of the different terms in $`H`$ we have developed a novel, analytical approach based on a perturbative diagonalization of $`H`$ using a suitable set of trial functions and using Löwdin partitioning.  Though we cannot expect accurate numerical results from such an approach it is an instructive complement for numerical methods, as we can clearly identify in the subband dispersion $`(𝐤_{})`$ the terms proportional to $`E_z`$ which are breaking the spin degeneracy. Neglecting in $`H_{𝐤𝐩}`$ remote bands like $`\mathrm{\Gamma }_8^c`$ and $`\mathrm{\Gamma }_7^c`$ we obtain for the SIA spin splitting of the HH states

$$_{\mathrm{HH}}^{\mathrm{SO}}(k_{})\pm \beta _1E_zk_{}^3.$$
(9)
In particular, we have no $`k`$ linear splitting (and $`\beta _20`$) if we restrict ourselves to the Luttinger Hamiltonian  which includes $`\mathrm{\Gamma }_8^c`$ and $`\mathrm{\Gamma }_7^c`$ by means of second order perturbation theory.  Accurate numerical computations  show that the dominant part of the $`k`$ linear splitting of the HH states is due to BIA. However, for typical densities this $`k`$ linear splitting is rather small. For the LH states we have

$$_{\mathrm{LH}}^{\mathrm{SO}}(k_{})\pm \beta _1E_zk_{}.$$
(10)
Thus we have a qualitative difference between the spin splitting of electron and LH states which is proportional to $`k_{}`$ and the splitting of HH states which essentially is proportional to $`k_{}^3`$. The former is most important in the low-density regime whereas the latter becomes negligible for small densities. Note that for 2D hole systems the first subband is HH like so that for low densities the SIA spin splitting is given by Eq. (9). In Eq. (7) the lengthy prefactors depend on the details of the geometry of the QW. Moreover, we have omitted a weak dependence on the direction of $`𝐤_{}`$. But the order of the terms with respect to $`k_{}`$ is independent of these details. It is crucial that, basically, we have

$$\alpha ,\beta _1,\beta _2\mathrm{\Delta }_0$$
(11)
with $`\mathrm{\Delta }_0`$ the SO gap between the bulk valence bands $`\mathrm{\Gamma }_8^v`$ and $`\mathrm{\Gamma }_7^v`$, i.e., we have no SIA spin splitting for $`\mathrm{\Delta }_0=0`$.  This can be most easily seen if we express $`H_{𝐤𝐩}`$ in a basis of orbital angular momentum eigenstates.

A more detailed analysis of our analytical model shows that both $`H_{6c}^{\mathrm{SO}}`$ and $`H_{8v}^{\mathrm{SO}}`$ stem from a third order perturbation theory for $`k_\pm `$, $`k_z=i_z`$, and $`eE_zz`$. This seems to be a rather high order. Nevertheless, the resulting terms are fairly large.  In agreement with Refs.  this is a simple argument to resolve the old controversy based on an argument by Ando  that spin splitting in 2D systems ought to be negligibly small because for bound states in first order we have $`E_z=0`$. We note that the present ansatz for the prefactors $`\alpha `$ and $`\beta _1,\beta _2`$ is quite different from the ansatz in Ref. . We obtain $`H_{6c}^{\mathrm{SO}}`$ and $`H_{8v}^{\mathrm{SO}}`$ by means of Löwdin partitioning of the Hamiltonian (7) whereas in Ref.  the authors explicitly introduced $`H_{6c}^{\mathrm{SO}}`$ into their model. Moreover, we evaluate the matrix elements of $`eE_zz`$ with respect to envelope functions for the bound states whereas in Ref.  the authors considered matrix elements of $`eE_zz`$ with respect to bulk Bloch functions. The latter quantities are problematic because they depend on the origin of the coordinate frame.

As an example, we show in Fig. 2 the self-consistently calculated  anisotropic subband dispersion $`_\pm (𝐤_{})`$, DOS effective mass $`m^{}/m_0`$, and spin splitting $`_+(𝐤_{})_{}(𝐤_{})`$ for a  grown GaAs/Al0.5Ga0.5As heterostructure. The calculation was based on a $`14\times 14`$ Hamiltonian ($`\mathrm{\Gamma }_8^c`$, $`\mathrm{\Gamma }_7^c`$, $`\mathrm{\Gamma }_6^c`$, $`\mathrm{\Gamma }_8^v`$, and $`\mathrm{\Gamma }_7^v`$). It fully took into account both SIA and BIA. The weakly divergent van Hove singularity of the DOS effective mass at the subband edge indicates that the $`k`$ linear splitting is rather small. (Its dominant part is due to BIA. ) Basically, the spin splitting in Fig. 2 is proportional to $`k_{}^3`$.

Only for the crystallographic growth directions  and  the hole subband states at $`k_{}=0`$ are pure HH and LH states. For low-symmetry growth directions like  and  we have mixed HH-LH eigenstates even at $`k_{}=0`$, though often the eigenstates can be labeled by their dominant spinor components.  The HH-LH mixing adds a $`k`$ linear term to the splitting (9) of the HH states, which often exceeds $`\beta _2E_zk_{}`$. However, this effect is still small when compared with the cubic splitting.

For a HH subband dispersion $`_\pm (k_{})=\mu k_{}^2\pm \beta _1E_zk_{}^3`$ we obtain for the densities $`N_\pm `$ in the spin-split subbands

$$N_\pm =\frac{1}{2}N_s\pm \frac{\beta _1E_zN_s}{\sqrt{2}\mu X}\sqrt{\pi N_s(64/X)}$$
(13)
with

$$X=1+\sqrt{14\pi N_s\left(\frac{\beta _1E_z}{\mu }\right)^2}.$$
(14)
The spin splitting according to Eq. (Rashba spin splitting in two-dimensional electron and hole systems) is substantially different from Eq. (4). For electrons and a fixed electric field $`E_z`$ but varying $`N_s`$ the difference $`\mathrm{\Delta }N=N_+N_{}`$ increases like $`N_s^{1/2}`$ whereas for HH subbands it increases like $`N_s^{3/2}`$. Using a fixed density $`N_s`$ but varying $`E_z`$ it is more difficult to detect the difference between Eqs. (4) and (Rashba spin splitting in two-dimensional electron and hole systems). In both cases a power expansion of $`\mathrm{\Delta }N`$ gives $`\mathrm{\Delta }N=a_1|E_z|+a_3|E_z|^3+𝒪(|E_z|^5)`$ with $`a_3<0`$ for electrons and $`a_3>0`$ for HH subbands.

The proportionality (11) is completely analogous to the effective $`g`$ factor in bulk semiconductors.  Lassnig  pointed out that the $`B=0`$ spin splitting of electrons can be expressed in terms of a position dependent effective $`g`$ factor $`g^{}(z)`$. In the following we want to discuss the close relationship between Zeeman splitting and $`B=0`$ spin splitting from a more general point of view. Note that in the presence of an external magnetic field $`𝐁`$ we have $`𝐤\times 𝐤=(ie/\mathrm{})𝐁`$ and the Zeeman splitting in the $`\mathrm{\Gamma }_6^c`$ conduction band can be expressed as

$$H_{6c}^Z=\frac{i\mathrm{}}{e}\frac{g^{}}{2}\mu _B𝐤\times 𝐤𝝈=\frac{g^{}}{2}\mu _B𝐁𝝈$$
(15)
with $`\mu _B`$ the Bohr magneton. Thus apart from a prefactor we obtain the Rashba term (1) from Eq. (15) by replacing one of the $`𝐤`$’s with the electric field $`𝐄`$. In the $`\mathrm{\Gamma }_8^v`$ valence band we have two invariants for the Zeeman splitting

$$H_{8v}^Z=2\kappa \mu _B𝐁𝐉+2q\mu _B𝐁𝓙.$$
(16)
Here, the first term is the isotropic contribution, and the second term is the anisotropic part. It is well-known that in all common semiconductors for which Eq. (16) is applicable the dominant contribution to $`H_{8v}^Z`$ is given by the first term proportional to $`\kappa `$ whereas the second term is rather small.  Analogous to $`\beta _1`$ and $`\beta _2`$ the isotropic $`𝐤𝐩`$ coupling between $`\mathrm{\Gamma }_8^v`$ and $`\mathrm{\Gamma }_6^c`$ contributes to $`\kappa `$ but not to $`q`$. The latter stems from $`𝐤𝐩`$ coupling to more remote bands such as $`\mathrm{\Gamma }_8^c`$ and $`\mathrm{\Gamma }_7^c`$.

Several authors  used an apparently closely related intuitive picture for the $`B=0`$ spin splitting which was based on the idea that the velocity $`v_{}=\mathrm{}k_{}/m^{}`$ of the 2D electrons is perpendicular to the electric field $`E_z`$. In the electron’s rest frame $`E_z`$ is Lorentz transformed into a magnetic field $`B`$ so that the $`B=0`$ spin splitting becomes a Zeeman splitting in the electron’s rest frame. However, this magnetic field is given by $`B=(v_{}/c^2)E_z`$ (SI units) and for typical values of $`E_z`$ and $`v_{}`$ we have $`B2\mathrm{}20\times 10^7`$ T which would result in a spin splitting of the order of $`5\times 10^9\mathrm{}5\times 10^5`$ meV. On the other hand, the experimentally observed spin splitting is of the order of $`0.1\mathrm{}10`$ meV. The $`B=0`$ spin splitting requires the SO interaction caused by the atomic cores. In bulk semiconductors this interaction is responsible for the SO gap $`\mathrm{\Delta }_0`$ between the valence bands $`\mathrm{\Gamma }_8^v`$ and $`\mathrm{\Gamma }_7^v`$ which appears in Eq. (11). The SO interaction is the larger the larger the atomic number of the constituting atoms. In Si we have $`\mathrm{\Delta }_0=44`$ meV whereas in Ge we have $`\mathrm{\Delta }_0=296`$ meV. Therefore, SIA spin splitting in Si quantum structures is rather small.

Recently, spin splitting in 2D systems has gained renewed interest because of an argument by Pudalov  which relates the metal-insulator transition (MIT) in low-density 2D systems with the SIA spin splitting. Based on the Rashba model  it was argued that the SIA spin splitting “results in a drastic change of the internal properties of the system even without allowing for the Coulomb interaction.” . However, as we have shown above, this argument is applicable only to electron and LH states. The MIT has been observed also in pure HH systems in, e.g., Si/SiGe QW’s.  As noted above, SO interaction and spin splitting in these systems are rather small,  so that it appears unlikely that here the broken inversion symmetry of the confining potential is responsible for the MIT. We note that in Si 2D electron systems the effective $`g`$ factor is enhanced due to many body effects.  It can be expected that similar effects are also relevant for the $`B=0`$ spin splitting.

In conclusion, we have analyzed the SIA spin splitting in 2D electron and hole systems. In 2D hole systems the splitting differs remarkably from its familiar counterpart in the conduction band. For electron states it increases linearly with in-plane wave vector $`k_{}`$ whereas the spin splitting of heavy hole states can be of third order in $`k_{}`$. We have discussed consequences of this behavior in the context of recent arguments on the origin of the metal-insulator transition observed in 2D systems.

The author wants to thank O. Pankratov, S. J. Papadakis, and M. Shayegan for stimulating discussions and suggestions.