text:
# A radio-microlensing caustic crossing in B1600+434?

## 1. Introduction

Optical microlensing has unambiguously been detected both in our own Galaxy (e.g. the EROS, MACHO and OGLE collaborations; see these proceedings), as well as in external galaxies (e.g. Q2237+0305; Irwin et al. 1989). However, no convincing detections of microlensing in other wavebands have been claimed thus far. In our Galaxy this is mostly due to the very low surface number density of bright sources, other than stars, that are compact enough (i.e. $`<`$few mas) to be significantly microlensed. For cosmological microlensing there are similar arguments. Because the angular Einstein radius of a point-mass scales approximately as $``$$`D^{1/2}`$, where $`D`$ is the distance to the lens, the lensed source must be $`<`$few $`\mu `$as in angular size to by microlensed by compact stellar-mass lenses at cosmological distances (few Gpc). It was thought until recently that only the optical to $`\gamma `$-ray emitting regions of quasars and AGNs had these extremely small angular scales and would still be observable over cosmological distances.

In this proceeding we will show, however, that ‘radio-microlensing’ has most likely been detected in the CLASS gravitational lens B1600+434 (Jackson et al. 1995) and that it promises to be a new and exciting technique for the study of massive compact objects in galaxies at intermediate redshifts.

## 2. Macro & microlensing of flat-spectrum radio sources

The JVAS/CLASS gravitational-lens survey (e.g. Browne et al. 1997) has discovered at least 17 radio-bright gravitational lens systems. All systems were selected to have a flat spectrum ($`\alpha `$$`<`$0.5 with $`S_\nu `$$``$$`\nu ^\alpha `$) between 5 and 1.4 GHz. This ensures that most of these sources are dominated by a compact radio-bright core. These sources are often variable on short time scales and in many cases show superluminal motion (e.g. Vermeulen & Cohen 1994).

Strong variability seems to imply that these radio sources are very compact. Based on a simple light-travel-time argument, one expects a source at a cosmological distance (few Gpc) to have an angular scale of only a few $`\mu `$as if it varies significantly on a time scale of say one month. However, superluminal motion with a high Doppler-boosting factor ($`𝒟`$) complicates this argument, because the intrinsic variability time scale is reduced by a factor $`𝒟^3`$, whereas the intrinsic angular scale of the source is reduced by only $`𝒟^{1/2}`$ (e.g. Lähteenmäki, Valtaoja & Wiik 1999). Based on variability arguments, the angular size of the source could therefore be severely underestimated if $`𝒟1`$. Furthermore, based on Rayleigh–Jean’s law and a typical surface brightness temperature of $`10^{1112}`$ K (e.g. Lähteenmäki et al. 1999), one expects an angular size as much as several tens of $`\mu `$as for those sources that have been observed in the JVAS/CLASS survey. This is much larger than the typical Einstein radius of a solar-mass object. Overall, it appears that microlensing of these sources is at most marginal, especially for the brighter sources (tens of mJy).

However, matters are not as bad as they seem! First of all, lensed sources are often magnified by the lensing potential. Their intrinsic flux-density is therefore less than the observed flux-density and, because surface brightness (temperature) is conserved in lensing, also their intrinsic angular source size is smaller. Second, flat-spectrum sources at high redshifts often exhibit jet-structures with superluminal motion (e.g. Vermeulen & Cohen 1994). If these jets contain many distinct ‘bullets’ or shock-fronts, their angular sizes can be very small, especially if they are significantly Doppler-boosted ($`𝒟1`$). We expect this to be the case near the core of these flat-spectrum sources, although components further along the jet might grow in size and have lower velocities. Hence, even though the integrated flux-density of the source might be large, suggesting a large angular size, in reality the likely presence of compact substructure with high Doppler-boostings makes this argument significantly weaker. The expected angular size of these flat-spectrum subcomponents is

$$\mathrm{\Delta }\theta \sqrt{\frac{S_{\mathrm{mJy}}(1+z_\mathrm{s})}{\mu 𝒟T_{12}}}\left(\frac{\lambda }{1\mathrm{cm}}\right)\mu \mathrm{as},$$

where $`T_{12}`$ is the surface brightness temperature of the component in units of $`10^{12}`$ K, $`\mu `$ is the magnification by the lens and $`z_\mathrm{s}`$ is the source redshift. For subcomponent flux-densities ($`S_{\mathrm{mJy}}`$) of a few mJy, an angular size of a few $`\mu `$as can be expected. This is small enough to be appreciably microlensed!

Below we will illustrate this by the CLASS gravitational lens system B1600+434, in which believe to have detected radio-microlensing of precisely this type of $`\mu `$as-scale jet-components (Koopmans & de Bruyn 2000; KdB00 hereafter).

## 3. CLASS B1600+434

The CLASS gravitational lens B1600+434 consists of two compact radio components separated by 1.4 arcsec (Fig.1). The lens galaxy is an edge-on spiral at a redshift of 0.41, whereas the source has a redshift of 1.59 (Koopmans, de Bruyn & Jackson 1998). Image A passes predominantly through the dark-matter halo of the lens galaxy ($``$6 kpc above the plane of the galaxy), whereas image B passes through its stellar component (i.e. disk/bulge). The extended optical emission around image A is thought to be associated with its host galaxy at $`z`$=1.59. Both radio components are compact ($`<`$1 mas) and have a flat radio spectrum. The source is highly variable at frequencies between 1.4 and 8.5 GHz (Koopmans et al. 1998).

To determine the time delay between the lensed images, we monitored B1600+434 with the VLA in A- and B-arrays at 8.5 GHz, during a period of about 8 months in 1998 (Koopmans et al. 2000). From these light curves we determined a time delay of 47$`{}_{6}{}^{}{}_{}{}^{+5}`$ days (68%), using the minimum-dispersion method from Pelt et al. (1996). The PRH method (Press, Rybicki & Hewitt 1992) gives the same time delay within the error range. Shifting the light curve of image B back in time over this time-delay and subsequently subtracting it from the light curve of image A, taking the proper flux-density ratio in to account, results in the difference light curve shown in Fig.2. The rms variability of the difference light curve (2.8%) is significantly larger than the expected rms variability due to noise only (1.1%). The difference light curve is consistent with the presence of non-intrinsic variability (i.e. ‘external’ variability) at the 14.6–$`\sigma `$ confidence level (KdB00). Because the time delay was determined using a minimum-dispersion method, time delays other than 47 days would yield an even larger rms variability in the difference light curve.

## 4. ‘External’ variability in B1600+434

What can be the origin of the presence of ‘external’ variability in the image light curves? There are several plausible causes:

* Scintillation by the Galactic ionized interstellar medium.
* Microlensing of the background source by massive compact objects in the lens galaxy.

Before proceeding, however, let us first summarize the results on B1600+434 from our VLA & WSRT monitoring data (until late 1999), that might be relevant in uncovering the true nature of the observed external variability (see KdB00):

1. The short-term rms variability (i.e. the 1st–order modulation index) of the lensed images A and B are 2.8% and 1.6% respectively.
2. The difference light curve has an rms variability of 2.8%.
3. The VLA 8.5-GHz image light curves, as well as the difference light curve, seem to show long-term ($``$1 d) variability, next to some faster variability. The precise time scales are hard to determine, because of the average light-curve sampling rate of once per 3.3 days.
4. Several 5-GHz WSRT 12-h integrated flux-density lightcurves of B1600+434 (June-July 1999), with a sampling rate of once every 5 min, show no evidence for variability $`>`$2% on time scales less than 12 hours.
5. Multi-frequency monitoring with the WSRT shows a decrease in the short-term rms variability from 5 to 1.4 GHz by a factor of $``$3.

Can this considerable body of multi-frequency data be explained in terms of Galactic scintillation (Section 4.1) or radio-microlensing (Section 4.2)?

### 4.1. Galactic Scintillation

From ‘standard’ scintillation theory (e.g. Narayan 1992), we expect a strong anti-correlation between the time-scale of variability and its amplitude (i.e. rms variability). According to the Taylor & Cordes (1993) model for the Galactic ionized medium, B1600+434 should be in the weak-scattering regime and have a variability time scale of at most $``$1 day for an observed rms variability of 1.6–2.8%. However, the WSRT data show no evidence for this short-term variability (point 4). In fact a number of flux-density variations appear to have relatively long time scales (several weeks; see point 3 and also Section 5). These time scales would normally correspond to an rms variability that is well below the noise level ($``$1%). Hence, the longer variability time scales are not compatible with the observed amplitude of variability in the lensed images (KdB00).

Beside this, the images are separated by only 1.4 arcsec but show a factor $``$1.75 difference in their rms variabilities (point 1). This either results from a difference of a factor $``$3.1 in their Galactic scattering measures, or a 90–$`\mu `$as scatter-broadening of image B in the lens galaxy at $`8.5\times (1+z_\mathrm{s})`$$``$12 GHz (KdB00). This requires a scattering measure $``$1 kpc m-20/3 in the lens galaxy, which is very large compared to typical lines-of-sight through our Galaxy. Scatter-broadening at a lower level can not be excluded, however.

Finally, in the case of scintillation the observed rms variability should increase towards longer wavelenghts if the source is extended (i.e. larger than the scattering disk). In the case of B1600+434, however, we find that the rms variability decrease by almost a factor $``$3 between 5 and 1.4 GHz (KdB00; point 5). One might argue that the source contains a compact component that scintillates and a more extended component that does not vary. The compact component would show less variability at 1.4 Ghz than at 5 GHz, because it would be in the strong-scattering regime. However, the time-scale of variability at 5 GHz for this component would remain only a few hours, whereas the 12-h WSRT light curves show no evidence for this (point 4).

Overall, we conclude that the observational evidence does not seem to favor scintillation as the explanation of the observed external variability in B1600+434 (KdB00). To save this hypothesis a considerable rethinking of the properties of the Galactic ionised medium is required, at least in the direction of B1600+434.

### 4.2. Radio-microlensing

How about radio-microlensing? There are different ‘ingredients’ that play a role:

* Source Structure: We assume that the source consist of a stationary core with a jet structure. This jet contains (several) condensations, either physical ‘bullets’ or shock-fronts. These jet-components move with near- or super-luminal motion and are compact ($`\mu `$as scale), especially near the core. Components further along the jet might grow in size.
* Compact Objects in the Lens Galaxy: The lens galaxy consists (partly) of massive compact objects. These objects create a magnification pattern on the source-plane, if they constitute a considerable fraction of the lens-galaxy mass inside a radius of 0.5–1 arcsec,

When these two ‘ingredients’ are combined, one will have jet-components that move over the magnfication pattern, being continuously magnified and/or demagnified, resulting in variability in the integrated flux density of the lensed images. The time-scale for these jet-component to cross the Einstein diameter of a point-mass is

$$\mathrm{\Delta }t4(1+z_\mathrm{s})\beta _{\mathrm{app}}^1\sqrt{(M/M_{})}\sqrt{(D_{\mathrm{ds}}D_\mathrm{s}/D_\mathrm{d}\mathrm{Gpc})}\text{ weeks,}$$

where $`z_\mathrm{s}`$ is the source redshift, $`\beta _{\mathrm{app}}`$ is the apparent jet-component velocity in units of $`c`$, $`M`$ is the lens mass, and $`D_{\mathrm{ds}}`$, $`D_\mathrm{s}`$ and $`D_\mathrm{d}`$ are the angular diameter distances between lens-source, observer-source and observer lens, respectively. More complex microlensing simulations show similar time-scales of the order of weeks to months for typical jet velocities of a few time $`c`$ and compact-object masses of $``$1 M. The longer observed variability time-scales for B1600+434 are therefore consistent with the radio-microlensing hypothesis (KdB00).

The difference in rms variability between the lensed images can be explained either by a difference in the compact-object mass function in their respective lines-of-sight and/or by moderate scatter-broadening ($`>`$few $`\mu `$as) of the image passing through the disk/bulge. If for example the average mass of compact objects in the disk/bulge is much lower than that for the halo line-of-sight, the resulting magnification pattern will show a more dense caustic network. This reduces the expected rms variability for a given angular source size (e.g. KdB00).

The strongest argument in favor of microlensing, however, is the frequency dependence of the observed rms variability. Using the constraints on the jet-structure that we derived from the 8.5-GHz VLA observations in 1998, we predict a decrease in modulation index from 5 to 1.4 GHz by a factor $``$3, if we assume that the jet-components are synchrotron self-absorbed (i.e. grow linearly with wavelength). This decrease has indeed been observed in the independent WSRT observations from 1998/1999 (Section 4; KdB00).

## 5. A caustic crossing in the radio?

During the 1999/2000 VLA A- and B-array configurations we have been monitoring B1600+434 at 8.5, 5 and 1.4 GHz (Koopmans et al. in prep). Some very preliminary results at 5 GHz are shown in Fig.3. This data show several strong ‘events’ (10–30%) at both 5 and 8.5 GHz, whereas they are absent in the 1.4-GHz light curve. What is their origin? All distinct events in image A do not re-occur in image B (after the time-delay of $``$47 days) and are therefore not intrinsic source variability.

Can the strongest feature (around day 80) be an ‘extreme scattering event’ (ESE; Fiedler et al. 1987)? Plasma clouds in our Galaxy move with typically $``$30 km/s. Event durations of several weeks then imply an angular size for these clouds of $``$1 mas, if they are located at $``$0.5 kpc distance from us, thereby covering the entire source. An event of $``$30% at 5 GHz would then be detectable at 1.4 GHz as well, whereas no evidence for this is seen in the VLA 1.4-GHz light curves. Moreover, these very rare ESE’s are not expect to occur several times in the same light curve of image A within a time span of only $``$150 days. Similarly, also scintillation does not show this extreme behavior in amplitude over time scales of several weeks (see Section 4.1).

Rather, we think that this event is caused by a $`\mu `$as-scale jet-component in the lensed source, which has recently been ejected from the core. This component subsequently moves over a single caustic in the magnification pattern that is created by compact objects in the line-of-sight towards lensed image A, causing observable variability in the integrated flux-density of the lensed image.

## 6. Conclusions

The CLASS gravitational lens B1600+434 shows strong ‘external’ variability in the VLA 5 and 8.5-GHz light curves of lensed image A, which passes predominantly through the halo of the lens galaxy. Neither Galactic scintillation nor Extreme Scattering Events (ESE) can satisfactorily explain all the observations. Only ‘radio microlensing’ can explain the observed variability time-scales and amplitudes, as well as its frequency dependence, without invoking extreme assumptions. It also offers a natural explanation for the strong ($``$30%) event that we more recently observed, as being a ‘radio-microlensing caustic crossing’.

Based on the 1998 VLA and 1998/9 WSRT data, we already attempted to place constraints on the mass-function of compact objects in the lens-galaxy halo, finding a lower limit of $`>`$0.5 M (KdB00). However, these constraints are still weak and depend strongly on assumptions about scatter-broadening in image B. The way foreward is to use only the lightcurve of image A to put constraints on the mass function in its line-of-sight. With the detection of distinct radio-microlensing events, this task has become a realistic possibility. We have started to work on this (Koopmans et al. in prep) and hopefully we will soon be able to place significantly better constraints on the mass function of compact objects in the lens-galaxy halo, bringing cosmological radio-microlensing up to Galactic standards!

#### Acknowledgments.

We thank Roger Blandford, Konrad Kuijken, Jean-Pierre Marquart, Penny Sackett and Jane Dennett-Thorpe for critical discussions during the course of this work. We also thank the participants of the ‘Microlensing 2000’ conference for positive feedback. LVEK and AGdeB acknowledge the support from an NWO program subsidy (grant number 781-76-101). This research was supported in part by the European Commission, TMR Program, Research Network Contract ERBFMRXCT96-0034 ‘CERES’.

## References

Browne, et al., 1997, in Cosmology with the new radio surveys, eds. M. Bremer and Jackson, N. J., p305

Fiedler, R. L., Dennison, B., Johnston, K. J. & Hewish, A. 1987, Nature, 326, 675

Irwin, M. J., Hewett, P. C., Corrigan, R. T., Jedrzejewski, R. I. & Webster, R. L. 1989, AJ, 98, 1989

Jackson, N., et al. 1995, MNRAS, 274, L25

Koopmans, L. V. E., de Bruyn, A. G. & Jackson, N. 1998, MNRAS, 295, 534

Koopmans, L. V. E., de Bruyn, A.G., Xanthopoulos, E., Fassnacht, C.D. 2000, A&A, 356, 391

Koopmans, L. V. E. & de Bruyn, A. G. 2000, A&A, in press, astro-ph/0004112 \[KdB00\]

Lähteenmäki, A., Valtaoja, E. & Wiik, K. 1999, ApJ, 511, 112

Narayan, R. 1992, Phil. Trans. Roy. Soc., 341, 151

Pelt, J., Kayser, R., Refsdal, S. & Schramm, T. 1996, A&A, 305, 97

Press, W. H., Rybicki, G. B. & Hewitt, J. N. 1992, ApJ, 385, 416

Taylor, J. H. & Cordes, J. M. 1993, ApJ, 411, 674

Vermeulen, R. C. & Cohen, M. H. 1994, ApJ, 430, 467


format: text

text:
# The Euler Equations and Non-Local Conservative Riccati Equations

## 1 Solving on characteristics

We will solve now the nonlocal Riccati equation on characteristics. We start with an auxiliary problem. Let $`\varphi `$ solve

$$_\tau \varphi +v\varphi =\varphi ^2$$

together with

$$v(x,\tau )=\varphi (x,\tau )+\frac{1}{|Q|}_Q\varphi (x,\tau )𝑑x.$$

We will consider initial data that are smooth, periodic and have zero mean,

$$_Q\varphi _0(x)𝑑x=0.$$

We will also assume that the curl $`\zeta =\frac{v_2}{x}\frac{v_1}{y}`$ obeys

$$_\tau \zeta +v\zeta =\left(\varphi \frac{3}{|Q|}_Q\varphi (x,\tau )𝑑x\right)\zeta .$$

Passing to characteristics

$$\frac{dY}{d\tau }=v(Y,\tau )$$
(15)
we integrate and obtain

$$\varphi (Y(a,\tau ),\tau )=\frac{\varphi _0(a)}{1+\tau \varphi _0(a)}$$

valid as long

$$\underset{aQ}{inf}(1+\tau \varphi _0(a))>0.$$

We need to compute

$$\overline{\varphi }(\tau )=\frac{1}{|Q|}_Q\varphi (x,\tau )𝑑x.$$

The Jacobian

$$J(a,\tau )=Det\left\{\frac{Y}{a}\right\}$$

obeys

$$\frac{dJ}{d\tau }=h(a,\tau )J(a,\tau )$$

where

$$h(a,\tau )=\varphi (Y(a,\tau ),\tau )\overline{\varphi }(\tau ).$$

Initially the Jacobian equals to one, so

$$J(a,\tau )=e^{_0^\tau h(a,s)𝑑s}.$$

So

$$J(a,\tau )=e^{_0^\tau \overline{\varphi }(s)𝑑s}\mathrm{exp}(_0^\tau \frac{d}{ds}\mathrm{log}(1+s\varphi _0(a))𝑑s)$$

and thus

$$J(a,\tau )=e^{_0^\tau \overline{\varphi }(s)𝑑s}\frac{1}{1+\tau \varphi _0(a)}.$$

The map $`aY(a,\tau )`$ is one and onto. The change of variables formula gives

$$_Q\varphi (x,\tau )𝑑x=_Q\varphi (Y(a,\tau ),t)J(a,\tau )𝑑a$$

and therefore

$$\overline{\varphi }(\tau )=e^{_0^\tau \overline{\varphi }(s)𝑑s}\frac{1}{|Q|}_Q\frac{\varphi _0(a)}{(1+\tau \varphi _0(a))^2}𝑑a.$$
(16)
Consequently

$$\frac{d}{d\tau }e^{_0^\tau \overline{\varphi }(s)𝑑s}=\frac{d}{d\tau }\frac{1}{|Q|}_Q\frac{1}{1+\tau \varphi _0(a)}𝑑a.$$

Because both sides at $`\tau =0`$ equal one, we have

$$e^{_0^\tau \overline{\varphi }(s)𝑑s}=\frac{1}{|Q|}_Q\frac{1}{1+\tau \varphi _0(a)}𝑑a$$
(17)
and, using (16),

$$\overline{\varphi }(\tau )=\left\{_Q\frac{\varphi _0(a)}{(1+\tau \varphi _0(a))^2}𝑑a\right\}\left\{_Q\frac{1}{1+\tau \varphi _0(a)}𝑑a\right\}^1$$
(18)
Note that the function $`\delta (x,\tau )=\varphi (x,\tau )\overline{\varphi }(\tau )`$ obeys

$$\frac{\delta }{\tau }+v\delta =\delta ^2+2\frac{1}{|Q|}_Q\delta ^2𝑑x2\overline{\varphi }\delta .$$

We consider now the function

$$\sigma (x,\tau )=e^{2_0^\tau \overline{\varphi }(s)𝑑s}\delta (x,\tau )$$

and the velocity

$$U(x,\tau )=e^{2_0^\tau \overline{\varphi }(s)𝑑s}v(x,\tau ).$$

Multiplying the equation of $`\delta `$ by $`e^{4_0^\tau \overline{\varphi }(s)𝑑s}`$ we obtain

$$e^{2_0^\tau \overline{\varphi }(s)𝑑s}\frac{\sigma }{\tau }+U\sigma =\sigma ^2+\frac{2}{|Q|}\sigma ^2𝑑x.$$

Note that

$$U=\sigma .$$

Now we change the time scale. We define a new time $`t`$ by the equation

$$\frac{dt}{d\tau }=e^{2_0^\tau \overline{\varphi }(s)𝑑s},$$
(19)
$`t(0)=0`$, and new variables

$$\gamma (x,t)=\sigma (x,\tau )$$

and

$$u(x,t)=U(x,\tau ).$$

Now $`\gamma `$ solves the nonlocal conservative Riccati equation

$$\frac{\gamma }{t}+u\gamma =\gamma ^2+\frac{2}{|Q|}\gamma ^2𝑑x$$
(20)
with periodic boundary conditions,

$$u=(\mathrm{\Delta })^1\left[^{}\omega +\gamma \right]$$
(21)
and

$$\frac{\omega }{t}+u\omega =\gamma \omega $$
(22)
The initial data are given by

$$\gamma _0(x)=\delta _0(x)=\varphi _0(x).$$

Using (17) and integrating the equation (19) we see that the time change is given by the formula

$$t=\left(\frac{1}{|Q|}\right)^2_Q_Q\frac{1}{\varphi _0(a)\varphi _0(b)}\mathrm{log}\frac{1+\tau \varphi _0(a)}{1+\tau \varphi _0(b)}dadb$$
(23)
Note that the characteristic system

$$\frac{dX}{dt}=u(X,t)$$

is solved by

$$X(a,t)=Y(a,\tau )$$

where $`Y`$ solves the system (15). This implies the formula

$$\gamma (X(a,t),t)=\alpha (\tau )\left(\frac{\varphi _0(a)}{1+\tau \varphi _0(a)}\overline{\varphi }(\tau )\right)$$
(24)
with

$$\alpha (\tau )=e^{2_0^\tau \overline{\varphi }(s)𝑑s}.$$
(25)
In view of (17), (18), (23), we have obtained a complete description of the general solution in terms of the initial data.

## 2 Blow up

Consider an initial smooth function $`\gamma _0(a)=\varphi _0(a)`$ and assume that it has mean zero and that its minimum is $`m_0<0`$. As it is evident from the explicit formula the blow up time for $`\varphi (Y(a,\tau ),\tau )`$ is

$$\tau _{}=\frac{1}{m_0}$$

and $`\varphi (Y(a,\tau ),\tau )`$ diverges to negative infinity for some $`a`$, and not at all for others. This of course does not necessarily mean that $`\gamma `$ blows up in the same fashion. Let us discuss the simplest case, in which the minimum is attained at a finite number of locations $`a_0`$, and near these locations, the function $`\varphi _0`$ has non-vanishing second derivatives, so that locally

$$\varphi _0(a)m_0+C|aa_0|^2$$

for $`0|aa_0|r`$. Then it follows that the integral

$$\frac{1}{|Q|}\frac{da}{ϵ^2+\varphi _0(a)m_0}$$

behaves like

$$\frac{1}{|Q|}\frac{da}{ϵ^2+\varphi _0(a)m_0}\mathrm{log}\left\{\sqrt{1+\left(\frac{Cr}{ϵ}\right)^2}\right\}$$

for small $`ϵ`$. Taking

$$ϵ^2=\frac{1}{\tau }\frac{1}{\tau _{}}$$

we deduce that, for these kinds of initial data

$$e^{_0^\tau \overline{\varphi }(s)𝑑s}\mathrm{log}\left\{\sqrt{1+\frac{C}{\tau _{}\tau }}\right\}.$$

For the same kind of functions and small $`(\tau _{}\tau )`$, the integral

$$\frac{1}{|Q|}_Q\frac{\varphi _0(a)}{(1+\tau \varphi _0(a))^2}𝑑a\frac{C}{\tau _{}\tau }$$

and $`t(\tau )`$ has a finite limit $`tT_{}`$ as $`\tau \tau _{}`$. The average $`\overline{\varphi }(\tau )`$ diverges to negative infinity,

$$\overline{\varphi }(\tau )\frac{C}{\tau _{}\tau }\left[\mathrm{log}\left\{\sqrt{1+\frac{C}{\tau _{}\tau }}\right\}\right]^1.$$

The prefactor $`\alpha `$ becomes vanishingly small

$$\alpha (\tau )(\mathrm{log}(\tau _{}\tau ))^2$$

and (24) becomes

$$\gamma (X(a,t),t)(\mathrm{log}(\tau _{}\tau ))^2\left(\frac{\varphi _0(a)}{1+\tau \varphi _0(a)}\overline{\varphi }(\tau )\right).$$

If the label is chosen so that $`\varphi _0(a)>0`$ then the first term in the brackets does not blow up and $`\gamma `$ diverges to plus infinity. If the label is chosen at the minimum, or nearby, then the first term in the brackets dominates and the blow up is to negative infinity, as expected from the ODE. From the equation (19)

$$(\alpha (\tau ))^1d\tau =dt$$

it follows that

$$T_{}t(\tau _{}\tau )\left(1+\mathrm{log}\left(\frac{1}{\tau _{}\tau }\right)\right)^2$$

and the asymptotic behavior of the blow up in $`t`$ follows from the one in $`\tau `$. We end by addressing a question that was at some point raised by numerical simulations: can there be a one-sided blow up? From the representation (24) of the solution it follows that

$$M(t)\overline{\varphi }(\tau )e^{2_0^\tau \overline{\varphi }(s)𝑑s}$$

holds for $`M(t)=sup_x\gamma (x,t)`$. If one would assume that, up to the putative blow up

$$M(t)C$$

with some fixed constant $`C`$, then it would follow that

$$\frac{d}{d\tau }e^{2_0^\tau \overline{\varphi }(s)𝑑s}2C$$

and integrating between $`\tau `$ and $`\tau _{}`$ that

$$e^{2_0^\tau \overline{\varphi }(s)𝑑s}2C(\tau _{}\tau ).$$

This in turn would imply that $`T_{}=\mathrm{}`$ and therefore no blow up for $`\gamma `$ can occur in finite $`t`$. So the answer is that for no initial datum can there exist a one sided blow up for $`\gamma `$.

Acknowledgments I thank J. Gibbon and K. Ohkitani for showing me their numerical work prior to its publication. This work was partially supported by NSF-DMS9802611 and by the ASCI Flash Center at the University of Chicago under DOE contract B341495.


text:
# Gluon and quark distributions in large 𝑁_𝑐 QCD: theory vs. phenomenology

## Abstract

We present the QCD analysis of the behavior of quark and gluon distributions in the limit of a large number of colors $`N_c`$. We show that the results agree qualitatively with the phenomenological data. In particular we draw attention to the fact that the polarized gluon distribution is suppressed at large $`N_c`$ compared to the unpolarized one.

aJoint Institute for Nuclear Research, Dubna, 141980, Russia
bInstitut für Theoretische Physik II, Ruhr-Universität Bochum,
D-44780 Bochum, Germany
cPetersburg Nuclear Physics Institute, Gatchina, St. Petersburg, 188350, Russia

The limit of large number of colors of QCD $`N_c\mathrm{}`$ and the $`1/N_c`$ expansion  are rather useful theoretical tools for analyzing the nonperturbative effects of QCD. Although in the real world we have $`N_c=3`$ which cannot be considered as a too large parameter, in many cases the large $`N_c`$ limit is helpful for understanding the dynamics of strong interactions. In particular, it was found that at large $`N_c`$ baryons can be described by a mean field solution of some effective meson action . Although the exact form of this action (depending on bilocal meson fields) is not known, certain conclusions can be made on the large $`N_c`$ behavior of various quantities.

In particular, in the limit of large number of QCD colors $`N_c`$ the twist-two parton distributions in nucleon typically behave as

$$f(x)=N_c^2\varphi (N_cx).$$
(1)
Here the function $`\varphi (z)`$ is stable in the large $`N_c`$ limit (it has the order $`O(N_c^0)`$ and its argument $`z=N_cx`$ is assumed to be fixed in the large $`N_c`$ limit).

One way to derive the large $`N_c`$ behavior (1) is to check that the moments of parton distributions corresponding to (1)

$$M_j=\underset{1}{\overset{1}{}}𝑑xx^{j1}f(x)=N_c^{j+2}\underset{N_c}{\overset{N_c}{}}𝑑yy^{j1}\varphi (y)=O(N_c^{j+2})$$
(2)
behave exactly as the matrix elements of the twist two operators $`O_{\mu _1\mathrm{}\mu _j}`$

$$M_j=n^{\mu _1}\mathrm{}n^{\mu _j}N,p|O_{\mu _1\mathrm{}\mu _j}|N,p=O(N_c^{j+2}).$$
(3)
Here $`n`$ is a light-cone vector $`n_\mu n^\mu =0`$ normalized by condition

$$(np)=1.$$
(4)
Taking nucleon at rest $`p=(m_N,0,0,0)`$ and keeping in mind that the mass of the nucleon at large $`N_c`$ behaves as

$$m_N=O(N_c)$$
(5)
we see from (4) that the nonvanishing components of $`n`$ behave as $`n=O(N_c^1)`$ so that $`j`$ factors $`n^{\mu _1}n^{\mu _2}\mathrm{}n^{\mu _j}`$ contribute explicitly $`O(N_c^j)`$ to the behavior of (3) and in order to derive (3) we must check that

$$N,p|O_{\mu _1\mathrm{}\mu _j}|N,p=O(N_c^2).$$
(6)
This can be done in a straightforward way using the rules of large $`N_c`$ counting , . For example, in the case of the quark operators $`O_{\mu _1\mathrm{}\mu _j}`$ one order of $`N_c`$ comes from the summation over color indices in operator $`O_{\mu _1\mathrm{}\mu _j}`$ and another order of $`N_c`$ is due to the relativistic normalization of the nucleon state

$$N,p^{}|N,p=2p^0\delta ^{(3)}(𝐩𝐩^{}).$$
(7)
Here for a nucleon at rest we have $`p^0=m_N=O(N_c)`$ which gives the second order of $`N_c`$ contributing to (6). Thus (3) is checked.

Alternatively one can derive the large $`N_c`$ behavior (1) of parton distributions directly following the logic of paper  (the derivation of  was formulated in terms of the chiral quark soliton model but its generalization to the “true” effective action of large $`N_c`$ QCD formulated in terms of bilocal meson fields is straightforward).

As an illustration of the large $`N_c`$ behavior (1) one can take the baryon number sum rule

$$\underset{0}{\overset{1}{}}\left[u(x)+d(x)\overline{u}(x)\overline{d}(x)\right]𝑑x=N_c$$
(8)
which obviously agrees with the large $`N_c`$ structure (1) of quark distributions and with the behavior of moments (2).

The physics standing behind the large $`N_c`$ structure (1) becomes rather transparent in the nonrelativistic quark model . Naively the nucleon consists of $`N_c`$ quarks and it is natural to expect that its momentum (in the infinite momentum frame) is more or less proportionally shared by all $`N_c`$ quarks. This explains why we keep $`N_cx`$ fixed in the limit of large $`N_c`$ in (1). The prefactor of $`N_c^2`$ in (1) is fixed e.g. by the normalization condition (8).

As for the spin structures in this model, the proton consists of $`(N_c+1)/2`$ $`u`$-quarks and $`(N_c1)/2`$ $`d`$-quarks. The colored part of the wave function is antisymmetric. Hence the spin-isospin wave function must be symmetric (we assume that all quarks have the same orbital momentum $`l=0`$). Hence the spin wave function of all $`u`$ quarks is symmetric and the $`(N_c+1)/2`$ quarks have the maximal total spin $`S_u`$. The same is valid for the spin of $`(N_c1)/2`$ $`d`$-quarks, i.e.

$$S_u=\frac{1}{4}(N_c+1),S_d=\frac{1}{4}(N_c1).$$
(9)
Interpreting the spin carried by $`u`$ and $`d`$ quarks as the moments of polarized quark distributions we see that

$$\underset{0}{\overset{1}{}}𝑑x\mathrm{\Delta }u(x)\underset{0}{\overset{1}{}}𝑑x\mathrm{\Delta }d(x)𝑑x=O(N_c)$$
(10)
in agreement with the general structure of the parton distributions (1).

In order to obtain the proton state with spin $`S=1/2`$ and isospin $`T=1/2`$, $`S_u`$ and $`S_d`$ must be opposite to each other so that $`O(N_c)`$ part of the spin is cancelled and we are left with

$$\underset{0}{\overset{1}{}}\left[\mathrm{\Delta }u(x)+\mathrm{\Delta }d(x)\right]𝑑x=O(N_c^0).$$
(11)
This means that in the large $`N_c`$ limit function $`\mathrm{\Delta }u+\mathrm{\Delta }d`$ is suppressed compared to (1) and has the form

$$\mathrm{\Delta }u(x)+\mathrm{\Delta }d(x)=N_c\varphi ^{}(N_cx).$$
(12)
The above derivation of this suppression was rather heuristic: we checked only the first moment and even this moment was estimated relying on the naive quark model. Howefer, a straightforward QCD analysis reproduces this result. One way to see it is to note that behavior (12) was derived in  in the context of chiral models but can be trivially generalized to the “true meson effective action” of large $`N_c`$ QCD. Alternatively one can check the behavior of the moments in the spirit of . Note that in both approaches the suppression (12) is essentially due to the so called spin-flavor symmetry  of the mean field solution for the baryon (the mean field approximation is justified at large $`N_c`$).

Actually not only $`\mathrm{\Delta }u+\mathrm{\Delta }d`$ but also some other flavor and polarization combinations of parton distributions have this suppressed large $`N_c`$ behavior:

$$f_{\mathrm{subleading}}(x)=N_c\varphi ^{}(N_cx).$$
(13)
For example, the isovector unpolarized distribution $`ud`$ is also of this type, which can be seen from the isospin sum rule

$$\underset{0}{\overset{1}{}}\left\{\left[u(x)d(x)\right][\overline{u}(x)\overline{d}(x)]\right\}𝑑x=1.$$
(14)
In the table we collect the relevant combinations of twist-two parton distributions. The main result of the large $`N_c`$ counting is that the isoscalar unpolarized and isovector polarized distributions are dominant in the $`1/N_c`$ limit whereas the isovector unpolarized and isoscalar polarized are suppressed (within the chiral quark soliton model this was established in , the generalization to the large $`N_c`$ QCD is straightforward).

| parton distributions | leading  $`N_c^2\varphi (xN_c)`$ | suppressed  $`N_c\varphi ^{}(xN_c)`$ |
| --- | --- | --- |
| unpolarized | $`u+d`$ | $`ud`$ |
| longitudinally  polarized | $`\mathrm{\Delta }u\mathrm{\Delta }d`$ | $`\mathrm{\Delta }u+\mathrm{\Delta }d`$ |
| transversity | $`\delta u\delta d`$ | $`\delta u+\delta d`$ |
| gluon | $`G`$ | $`\mathrm{\Delta }G`$ |

Note that for the antiquark distributions the large $`N_c`$ counting is the same as for the quark ones (at this point instead of the nonrelativistic quark model it is better to rely on the arguments of  or on the direct QCD large $`N_c`$ counting).

As for the gluon distribution, its large $`N_c`$ behavior can be estimated simply by generating the gluons through evolution from quarks. The one-loop DGLAP kernel $`\alpha _sP_{qg}`$ corresponding to this radiation is of order

$$\alpha _sP_{qg}\alpha _sC_F=\alpha _s\frac{N_c^21}{2N_c}\alpha _sN_c=O(N_c^0)$$
(15)
since $`\alpha _s=O(N_c^1)`$. This means that the gluon distributions inherit the large $`N_c`$ behavior of isoscalar quark distributions (1). It should be stressed that here we keep in mind the domain $`x1/N_c`$ of “small but not too small $`x`$”. The region corresponding to the physics of really small $`x`$ is beyond our consideration.

Now let us turn to the phenomenological consequences of the large $`N_c`$ counting. For the parton distributions whose shape is determined with a certain reliability ($`u,d,\mathrm{\Delta }u,\mathrm{\Delta }d`$) the large $`N_c`$ behavior

$$ud\frac{1}{N_c}(u+d),$$
(16)
$$\mathrm{\Delta }u+\mathrm{\Delta }d\frac{1}{N_c}(\mathrm{\Delta }u\mathrm{\Delta }d)$$
(17)
really seems to be in agreement with the experimental data.

Indeed, integrating both sides of (16) (weighted with $`x`$) and (17) over $`x`$ we arrive with the recent global fit of unpolarized data  at $`Q^2=2`$GeV2 at

$$\frac{_0^1𝑑xx\left[(ud)+(\overline{u}\overline{d})\right]}{_0^1𝑑xx\left[(u+d)+(\overline{u}+\overline{d})\right]}\frac{1}{3.5}$$
(18)
and for polarized data at

$$\frac{_0^1𝑑x\left[(\mathrm{\Delta }u+\mathrm{\Delta }d)+(\mathrm{\Delta }\overline{u}+\mathrm{\Delta }\overline{d})\right]}{_0^1𝑑x\left[(\mathrm{\Delta }u\mathrm{\Delta }d)+(\mathrm{\Delta }\overline{u}\mathrm{\Delta }\overline{d})\right]}=\frac{2\mathrm{\Delta }\mathrm{\Sigma }+a_8}{3g_A}\frac{1}{2.8}$$
(19)
Here we used the values $`g_A=1.26`$, $`a_8=0.58`$, $`\mathrm{\Delta }\mathrm{\Sigma }=0.39`$ . The results agree with the $`O(1/N_c)`$ suppression (16), (17).

An interesting result of the above analysis is that $`\mathrm{\Delta }G`$ is $`1/N_c`$ suppressed compared to $`G`$. A naive (but still justified for $`1/N_c`$ counting) way to understand this is again to use the radiative mechanism of the gluon generation: since the polarized $`\mathrm{\Delta }u+\mathrm{\Delta }d`$ distribution is $`1/N_c`$ suppressed in the sense of (13) we expect that $`\mathrm{\Delta }G`$ has the same large $`N_c`$ behavior. A more scientific way to derive the large $`N_c`$ suppression of $`\mathrm{\Delta }G`$ would be to make use of the fact that in the large $`N_c`$ limit the nucleon can be effectively described by a sort of Hartree (mean field) approach in terms of bilocal meson fields (including the composite fields made of gluons). The mean field appearing in this Hartree approach is believed to have mixed spin-flavor symmetry . In the gluon sector this symmetry of the mean field solution simply reduces to the pure rotational symmetry which leads to the suppression of the polarized gluon distribution.

Thus

$$\frac{\mathrm{\Delta }G}{G}\frac{1}{N_c}.$$
(20)
This is consistent with the first direct experimental estimation  of this ratio in the interval $`0.06<x<0.28`$

$$\frac{_{0.06}^{0.28}𝑑x\mathrm{\Delta }G}{_{0.06}^{0.28}𝑑xG}=0.41\pm 0.18$$
(21)
In the case of antiquark polarized distributions we have

$$\mathrm{\Delta }\overline{u}+\mathrm{\Delta }\overline{d}\frac{1}{N_c}(\mathrm{\Delta }\overline{u}\mathrm{\Delta }\overline{d}).$$
(22)
This leads us to the conclusion that in the large $`N_c`$ limit

$$\mathrm{\Delta }\overline{u}\mathrm{\Delta }\overline{d}.$$
(23)
This is supported by calculations in the chiral quark soliton model  where the large flavor asymmetry $`\mathrm{\Delta }\overline{u}\mathrm{\Delta }\overline{d}`$ was obtained (for the numerical parametrization see ). Also note that this result agrees with the instanton model prediction , with the phenomenological approach of  based on the “Pauli blocking principle” and with the statistical model of . Our opinion is that at the current level of the uncertainties the large $`N_c`$ behavior $`\mathrm{\Delta }\overline{u}\mathrm{\Delta }\overline{d}`$ seems more reliable than the ad hoc assumption $`\mathrm{\Delta }\overline{u}=\mathrm{\Delta }\overline{d}`$ which is often used in the analysis of phenomenological data. An indication of a nonzero value for $`\mathrm{\Delta }\overline{u}\mathrm{\Delta }\overline{d}`$ was also observed in .

Concerning the transversity distribution, the large $`N_c`$ QCD predicts $`\delta u=\delta d`$ up to $`1/N_c`$ corrections (see the table). The different signs of $`\delta u`$ and $`\delta d`$ manifest themselves in the opposite signs of $`\pi ^+`$ and $`\pi ^{}`$ azimuthal spin asymmetries in semi-inclusive DIS  and in high $`p_T`$ inclusive pion production  which are observed experimentally .

Certainly, the parameter $`1/N_c=1/3`$ is not too small and looking at a specific quantity one cannot unambiguously conclude whether we really observe the $`1/N_c`$ suppression or not. However, the above examples (see also  for the analysis of the large $`N_c`$ bounds on transversity distributions) give a certain evidence that the experimental data for parton distributions are in a rather reasonable agreement with the large $`N_c`$ counting of QCD. As a particular result of this investigation one can assert that the polarized gluon distribution is suppressed by $`1/N_c`$ compared to the unpolarized one.

Acknowledgment. We are grateful to M. Glück, V.Yu. Petrov, M.V. Polyakov, E. Reya, A. Shuvaev and O. Teryaev for discussions. This work was partially supported by RFBR Grant 00-02-16696, by DFG and by BMBW. A.E. and P.P. are thankful to the Institute for Theoretical Physics of Ruhr University Bochum for warm hospitality.


text:
# Magnetoresistance and electronic structure of asymmetric 𝐆𝐚𝐀𝐬/𝐀𝐥_0.3𝐆𝐚_0.7𝐀𝐬 double quantum wells in the in-plane/tilted magnetic field

## I Introduction

In an idealized, infinitely narrow two-dimensional (2D) system the in-plane component of the magnetic field couples only to the electronic spin degree of freedom. In real samples, however, the finite size of a 2D system in the growth ($`\widehat{z}`$) direction often leads to strong orbital effects of the in-plane field. For example, 2D-2D tunneling studies of a system consisting of two nearby narrow quantum wells, which is the simplest structure with non-trivial growth direction degree of freedom, have revealed dramatic effects of the in-plane field on the Fermi surface topology. In this weakly-coupled double quantum well sample the in-plane field, $`B_{}`$, in effect displaces the origin of the two layer Fermi circles by $`|e|B_{}d/\mathrm{}`$ and allows a sweep of one past the other at a critical field $`B_c=2\mathrm{}k_F/|e|d`$, where $`d`$ is the layer separation and $`k_F`$ is the individual layer Fermi wavevector.

The crossing of Fermi surfaces is replaced with more complicated patterns in samples with strongly coupled quantum wells in which the tunneling rate dominates the electronic scattering rate within a well. In a simple tight-binding model a partial energy gap, equal to the bonding-antibonding gap at zero field, opens at the wavevector corresponding to the crossing point of the displaced dispersion curves for the uncoupled layers. While the upper subband above the energy gap maintains a nearly parabolic shape, a saddle point develops in the lower energy subband. With increasing $`B_{}`$ the bottom of the upper subband moves above the Fermi energy at a critical field $`B_{c,1}`$, resulting in a sudden drop in the density of states. At the second critical field, $`B_{c,2}`$, the saddle point of the lower subband approaches the Fermi energy and the density of states diverges. The magnetoresistance oscillation observed on coupled double and triple quantum wells represents a striking manifestation of the two distinct van Hove singularities in the $`B_{}`$-dependent density of states. For these samples the tight-binding model has provided an accurate quantitative estimate for the critical in-plane fields. Bilayer systems realized in wide single quantum wells also display the magnetoresistance oscillation, however, the orbital effects of the in-plane magnetic field are more complex here and cannot be captured by the simple tight-binding approximation. Instead, the numerical self-consistent-field technique has proven successful for wide quantum wells.

The importance of orbital effects of the in-plane magnetic field has recently been emphasized also in the context of metal-insulator transition studies in Si-MOSFET and GaAs/Al0.3Ga0.7As single heterojuncions. The observed dramatic response to $`B_{}`$ in both metallic and insulating phases has been attributed to the spin coupling to in-plane field as well as to the distortion of carrier Fermi surfaces. A detailed quantitative understanding of these phenomena requires to incorporate the orbital effects of $`B_{}`$ using non-perturbative approaches.

In this paper we present a study of magnetotransport properties of GaAs/Al0.3Ga0.7As heterojunctions with an additional thin Al0.3Ga0.7As barrier introduced into the GaAs buffer. The structures were specially designed to combine properties of single-junction and bilayer systems. A comparison between theory and experiment in our samples constitutes an excellent test for the reliability of the numerical self-consistent-field technique applied to systems with complex growth direction geometries. In Section II we present measured resistance oscillation with in-plane field and Shubnikov-de Haas (SdH) oscillations recorded at magnetic fields slightly tilted from the 2D layer plane. Theoretical calculations of the $`B_{}`$-dependent electronic structure are used, in Section III, for quantitative interpretation of the experimental data. Our conclusions are summarized in Section IV.

## II Experimental results

We studied modulation-doped GaAs/Al0.3Ga0.7As heterojunctions grown by molecular-beam epitaxy on (100) oriented semi-insulating GaAs substrates and patterned in the standard Hall-bar samples. A thin barrier formed by eight monolayers of Al0.3Ga0.7As was grown inside the GaAs buffer producing two coupled quantum wells, one of a distorted rectangular shape and the other one with a nearly triangular geometry, as shown in insets of Fig. 1. We chose the thickness and position of the barrier to achieve high population of the second (antibonding) subband. The particular structure design was based on zero-field numerical self-consistent simulations that predict minimum bonding-antibonding gap for barrier positioned near the node of the first excited state wavefunction (peak of the lowest energy state wavefunction) of the original single-junction. Note that unlike the conventional bilayers realized in double quantum wells or wide single wells our system is intrinsicly, strongly anisotropic. At zero magnetic field the lower (bonding) subband wavefunction has dominant weight in the rectangular well while electrons from the higher (antibonding) subband are more likely to occupy the triangular quantum well. Parameters of the two studied samples are summarized in Table I. The partial occupations of the bonding and antibonding subbands at $`B_{}=0`$ were obtained by Fourier analysis of the low-field SdH oscillations. The Hall measurement at low perpendicular magnetic fields provided an independent check for the total carrier density and was used to determine sample mobilities. The insets in Fig. 1 show the confining potential profiles and wavefunctions of the occupied subbands for two samples.

Magnetotransport data were collected at the temperature 0.45 K and in magnetic fields ranging from 0 to 23 T, using both dc and low-frequency (13 Hz) ac techniques. First, we measured the resistance $`R_{xx}`$ for a magnetic field precisely parallel with the 2D layer plane. The recorded magnetoresitance, $`\mathrm{\Delta }R_{xx}/R_{xx}=(R_{xx}(B_{})R_{xx}(0))/R_{xx}(0)`$, is plotted in Fig. 1. Consistent with the bilayer nature of the studied samples we observe an oscillation on the magnetoresistance trace with the lower critical field $`B_{c,1}`$ and the upper critical field $`B_{c,2}`$. The values of critical fields are $`B_{c,1}=5.8`$ T, $`B_{c,2}=11.7`$ T for sample A and $`B_{c,1}=5.9`$ T, $`B_{c,2}=10.7`$ T for sample B, respectively. The in-plane magnetic field dependence of subband occupations, $`N_a(B_{})`$ and $`N_b(B_{})`$, were obtained from SdH oscillation data measured at different field-tilt angles. Typical SdH traces are shown in Fig. 2. In general, the oscillations are not periodic in $`1/B_{}`$ ($`B_{}`$ is the perpendicular component of the magnetic field) since both $`B_{}`$ and $`B_{}`$ vary when the field is swept at a fixed tilt angle.

Therefore, the standard Fourier technique does not apply. Instead, we use the measured distance between valleys surrounding individual peaks to identify subband densities. For $`B_{}<B_{c,1}`$ one type of oscillations, corresponding to the lower density (antibonding) subband, can be detected at small tilt angles while two subbands are clearly visible at higher angles. With increasing $`B_{}`$ the distribution of electrons between the two subbands changes, as shown in Fig. 3. Eventually, the antibonding subband is depopulated at $`B_{}=B_{c,1}`$. At higher in-plane fields two regimes can be distinguished. For $`B_{c,1}<B_{}<B_{c,2}`$ the SdH oscillations are periodic indicating a single occupied subband. Consistently, the corresponding density is equal to the total 2D density measured at zero in-plane magnetic field. The character of the oscillations changes abruptly at $`B_{}=B_{c,2}`$; the low in-plane field non-periodicity is recovered and the obtained carrier density is significantly reduced (see Fig. 3). These results suggest that the Fermi sea splits into two disconnected parts of densities that vary with $`B_{}`$. Only one Fermi surface, however, can be detected from our SdH data.

## III Theoretical analysis

The single-particle Hamiltonian for an electron confined in a GaAs/Al0.3Ga0.7As heterostructure and subjected to in-plane magnetic field can be written as

$$H=\frac{1}{2m^{}}(𝐩+|e|𝐀)^2+V_{conf}(z),$$
(1)
where $`m^{}`$ is the effective mass in the GaAs conduction band and the vector potential for $`B_{}`$ applied along the $`\widehat{y}`$-direction takes a form $`𝐀=(B_{}z,0,0)`$. The confining potential $`V_{conf}(z)`$ is constructed using the nominal growth parameters and includes Hartree and exchange-correlation potentials generated by the free carriers in the quantum wells. The Hartree potential is derived from the $`z`$-dependent density of electrons, $`\rho (z)`$, by numerical solution of the Poisson equation. The exchange-correlation term is calculated within the local-density approximation. In each loop of the self-consistent procedure we solve numerically the Schrödinger equation with the Hamiltonian (1) to get $`\rho (z)`$. Then, a new $`V_{conf}(z)`$ is constructed, which enters the next loop of the procedure until the self-consistency condition has been achieved. The resulting energy spectra have the form

$$E_i(𝐤)=E_{i,x}(k_x)+\frac{\mathrm{}^2}{2m}k_y^2,$$
(2)
where $`E_{i,x}(k_x)`$ depends on the sample geometry and on the magnitude of the in-plane magnetic field. An index $`i=a,b`$ distinguishes the antibonding and bonding subband states.

2D subband concentrations $`N_i`$ are proportional to the area enclosed by corresponding Fermi contours. The partial densities of states (DOS) $`g_i`$ and the cyclotron effective masses $`m_{c,i}`$ are related to the shape of the Fermi contours by the following expression:

$$g_i=\frac{m_{c,i}}{\pi \mathrm{}^2}=\frac{1}{2\pi ^2}\frac{dk}{|_kE_i|}.$$
(3)
The total density of states $`g=g_a+g_b`$. The distortion of Fermi lines can be probed experimentally by adding a weak perpendicular component of the magnetic field. For $`B_{}B_{}`$, the quantization of the in-plane component of the electron motion can be described in terms of quasiclassical Landau levels with the cyclotron effective mass given by Eq. (3). The degeneracy of the spin-unresolved levels is $`2|e|B_{}/h`$.

The evolution of Fermi surfaces in $`B_{}`$ is illustrated in Fig. 4 for sample B. Fig. 4(a) presents theoretical Fermi contours calculated for several selected values of $`B_{}`$, corresponding cyclotron effective masses are shown in Fig. 4(b). The deviations from the zero-field circles reflect the field-induced changes of $`E_{i,x}(k_x)`$ described for the case of the simple tight-binding model in Section I. The antibonding $`E_{a,x}(k_x)`$ is an asymmetric function of $`k_x`$, narrower than the free-electron parabola, with the minimum (bottom) moving to the Fermi energy which crosses at $`B_{}=B_{c,1}`$. The corresponding Fermi contour acquires the shape of a “lens”, i.e. the oval with the longer axis oriented in the $`\widehat{k}_y`$-direction. Due to this type of deformation the antibonding cyclotron effective mass $`m_{c,a}`$ is a decreasing function of $`B_{}`$.

The field-induced local maximum developing in the $`k_x`$-dependence of $`E_{b,x}(k_x)`$ causes elongation of the bonding Fermi contour in the $`\widehat{k}_x`$-direction. For higher fields, the Fermi line evolves into an asymmetric “peanut” shape. A neck connecting left and right parts of the “peanut” breaks as the saddle point reaches the Fermi energy at $`B_{}=B_{c,2}`$. For $`B_{}<B_{c,2}`$ the bonding effective mass $`m_{c,b}`$ grows with $`B_{}`$ and diverges at the saddle point. Above $`B_{c,2}`$, the Fermi contour splits into two approximately elliptic lines. Very similar results are obtained for sample B.

The position of an electron described by the state $`|i,k_x,k_y`$ in one or the other well is given by $`z_{i,x}=i,k_x,k_y|z|i,k_x,k_y`$. The coordinate $`z_{i,x}`$ is related to the $`k_x`$-component of the wave vector by

$$z_{i,x}=\frac{\mathrm{}k_x}{m\omega }\frac{v_{i,x}}{\omega },$$
(4)
$$v_{i,x}=\frac{1}{\mathrm{}}\frac{E_{i,x}(k_x)}{k_x},$$
(5)
where $`\omega =|e|B_{}/m^{}`$ and $`v_{i,x}`$ is the $`x`$-component of the in-plane group velocity. These relations imply that the larger, almost circular Fermi contour describes the electron layer in the left rectangular well while the smaller elongated oval corresponds to states in the right triangular well. The cyclotron effective masses in two wells are different: above $`B_{c,2}`$ the mass for the left well drops quickly close to the zero-field value $`m^{}`$, the mass for the triangular well also drops, but to substantially larger value.

Both the partial densities of electrons $`\rho _a(z)`$, $`\rho _b(z)`$ and their sum $`\rho (z)=\rho _a(z)+\rho _b(z)`$ are influenced by in-plane magnetic fields. In our structure with a hard-wall barrier the densities have sharp minima in the barrier and peaks inside the wells. The position of peaks is almost the same for all densities and practically field-independent, the magnetic field effect is represented by changes of relative magnitudes of the peaks. We will characterize the transfer of electrons between wells by $`\overline{z}_i(B_{})`$, the field-dependent centroids of electron densities. Results calculated for sample B are shown in Fig. 5. The center of mass of $`\rho (z)`$ is in the middle of the structure at $`B_{}=0`$ and moves closer to the interface as the field increases. This behavior is typical for single-junction structures. Below $`B_{c,1}`$, the centroid $`\overline{z}_a`$ of antibonding electrons is inside the triangular well for all fields. At zero-field, $`\overline{z}_a`$ is far from the barrier. The growing $`B_{}`$ transfers antibonding electrons from the triangular to rectangular well and, consequently, $`\overline{z}_a`$ is shifted towards the barrier. The electrons from the bonding subband exhibit opposite behavior in this range of fields: the magnetic field empties the antibonding subband and corresponding electrons, mostly from the triangular well, become a part of the bonding subband. As result, $`\overline{z}_b`$ moves from position close to the interface to the position near the barrier as $`B_{}`$ increases. Above $`B_{c,1}`$, when only the bonding subband is occupied, the transfer of electrons from the triangular to rectangular well continues. This is the case also above $`B_{c,2}`$ when the bonding Fermi contour is splitted into the left and right parts, as seen in Figs. 3 and 4.

The characteristic of the samples A and B derived from the magnetoresitance experiments are collected in Figs. 1 and 3, together with results of the self-consistent numerical simulation. In the case when theory cannot be directly compared with experiments only the theoretical results for sample B are presented.

We associate the sharp minima in the experimental magnetoresistance traces, presented in Fig. 1 with the sudden decrease of the density of states at $`B_{c,1}`$ that contributes to reduction of the scattering rate at the Fermi energy. The difference between the field-dependence of cyclotron effective masses $`m_{c,a}`$ and $`m_{c,b}`$, calculated for sample B and shown in Fig. 4(b), explains why the SdH oscillations of electrons in the antibonding subband are seen at lower tilt-angles than the oscillations of electrons from the bonding subband.

The increase of theoretical $`m_{c,b}`$ at in-plane fields below $`B_{c,2}`$ is related to the strong distortion of the Fermi contour, as seen in Fig 4. Only the bonding subband is occupied for $`B_{c,1}<B_{}<B_{c,2}`$ and, therefore, the SdH oscillations are periodic in $`1/B_{}`$. The decrease of their amplitude for $`B_{}`$ approaching $`B_{c,2}`$ is an experimental confirmation of calculated sharp increase of $`m_{c,b}`$. The theoretical value for the critical in-plane field precisely matches with $`B_{c,2}`$ obtained from experiment. The peak in the measured in-plane field magnetoresistance at $`B_{c,2}`$ can be understood as a consequence of the divergence of the DOS at Fermi energy or, in other words, of a zero group velocity of electrons at the saddle point.

The experimental magnetoresistance oscillation is clearly associated with van Hove singularities in the $`B_{}`$-dependent DOS. However, apart from these strong features $`R_{xx}`$ does not follow $`g`$ for in-plane fields less than $`B_{c,2}`$ (see Fig. 1). The differences between $`\mathrm{\Delta }R_{xx}/R_{xx}`$ and $`\mathrm{\Delta }g/g`$ curves stem from the in-plane field suppression of the coupling between wells, the transfer of electrons from the triangular to rectangular well, and the changes in the nature of electronic scattering.

At $`B_{}>B_{c,2}`$, the Fermi contour splits into two approximately elliptic lines. The larger, almost circular Fermi contour is for the rectangular well while the smaller elongated oval corresponds to states in the triangular well. Since electrons are localized in one or the other quantum well, the scattering rate becomes nearly independent of $`B_{}`$, and $`\mathrm{\Delta }R_{xx}/R_{xx}\mathrm{\Delta }g/g`$ applies, as seen from Fig. 1. Large cyclotron effective mass and low concentration of electrons in the triangular well explain while these states are not detectable by the SdH measurement; at most two weak oscillations can occur in available magnetic fields. The calculated density of electrons in the rectangular well is again in excellent quantitative agreement with SdH data.

## IV Conclusions

We have performed an experimental and theoretical study of asymmetric GaAs/Al0.3Ga0.7As heterostructure consisting of coupled rectangular and triangular quantum wells. The positions of a strong minimum and maximum in the measured in-plane field dependent resistance of the 2D electron system match precisely with calculated van Hove singularities in the DOS. Theoretical predictions for $`B_{}`$-dependent occupations of electronic subbands and quantum wells are in quantitative agreement with SdH data recorded at magnetic fields slightly tilted from the 2D layer plane. We conclude that the numerical self-consistent-field technique provides a realistic description of orbital effects of the in-plane field on 2D electron systems confined in semiconductor heterostructures with general growth-direction geometry.

To achieve similar level of accuracy for valence-band states with strong spin-orbit coupling represents a formidable challenge for future theoretical work. Particularly intriguing is the derivation of local-spin-density approximation for many-body states with unequal population of different spin-subbands. Recent remarkable studies of the metal-insulator transition in high mobility GaAs hole systems certainly give a strong motivation for pursuing this kind of research.

This work was supported by the Czech-French project Barrande 99011 and by the Grant Agency of the Czech Republic under Contract No. 202/98/0085.


text:
# Generalized Compressibility in a Glass Forming Liquid

## Abstract

We introduce a new quantity to probe the glass transition. This quantity is a linear generalized compressibility which depends solely on the positions of the particles. We have performed a molecular dynamics simulation on a glass forming liquid consisting of a two component mixture of soft spheres in three dimensions. As the temperature is lowered, the generalized compressibility drops sharply at the glass transition. Our results are consistent with the kinetic view of the glass transition, but not with an underlying second order phase transition.

The glass transition is still not well understood despite extensive study. Experimentally the glass transition occurs when the relaxation time exceeds the measurement time and particle motion appears to be arrested. It is characterized by both kinetic and thermodynamic features. In the supercooled liquid kinetic quantities such as the viscosity and relaxation time grow rapidly as the temperature is lowered. The glass transition is also reflected in thermodynamic quantities, e.g., the specific heat at constant pressure has a step–like form and the dielectric constant has a peak at a frequency dependent temperature.

Theoretically there have been two main approaches to the problem : dynamic and thermodynamic. The first category has been dominated by mode coupling theory (MCT) in which ideally the relaxation time diverges at a temperature $`T_C`$ above the experimental glass transition . The kinetic view has produced interesting and fruitful concepts such as the influence of the energy landscape on relaxation processes  and dynamic inhomogeneities . The thermodynamic viewpoint attributes the glass transition to an underlying phase transition hidden from direct experimental observation by extremely long relaxation times . In most scenarios there is an underlying second order phase transition associated with a growing correlation length which produces diverging relaxation times as well as diverging static susceptiblities . More recently Mezard and Parisi  have argued that the underlying transition is actually a random first order transition signaled by a discontinuity in the specific heat.

In an effort to better characterize the glass transition we introduce a novel probe which we call a generalized compressibility . Unlike the specific heat which monitors energy fluctuations, this linear compressibility is a function of the microscopic structure of the system: it depends solely on the positions of the particles and not on their previous history. It is a thermodynamic quantity in the sense that it is purely a function of the microstate of the system dictated by its location in phase space. It is easy to compute numerically, and it is simpler than the dielectric constant which involves both the translation and orientation of electric dipoles. By performing a molecular dynamics simulation of a two component system of soft spheres, we find that the linear generalized compressibility drops sharply as the temperature decreases below the glass transition temperature $`T_g`$. The drop becomes more and more abrupt as the measurement time increases.

We now derive expressions for the linear and nonlinear generalized compressibility. To probe the density fluctuations, we follow the approach of linear response theory and consider applying an external potential $`\frac{\mathrm{\Delta }P}{\rho _o}\varphi (\stackrel{}{r})`$ which couples to the local density $`\rho (\stackrel{}{r})=_{i=1}^N\delta (\stackrel{}{r}\stackrel{}{r}_i)`$ where $`\stackrel{}{r}_i`$ denotes the position of the $`i^{th}`$ particle. $`\rho _o`$ is the average density. $`\mathrm{\Delta }P`$ has units of pressure and sets the magnitude of the perturbation. $`\varphi (\stackrel{}{r})`$ is a dimensionless function of position that must be compatible with the periodic boundary conditions imposed on the system, i.e., it must be continuous across the boundaries, but is otherwise arbitrary. This adds to the Hamiltonian H of the system a term

$$U=\frac{\mathrm{\Delta }P}{\rho _o}_Vd^3r\varphi (\stackrel{}{r})\rho (\stackrel{}{r})=\frac{\mathrm{\Delta }P}{\rho _o}\underset{i}{}\varphi (\stackrel{}{r}_i)\frac{\mathrm{\Delta }P}{\rho _o}\rho _\varphi $$
(1)
where we have defined $`\rho _\varphi =_Vd^3r\varphi (\stackrel{}{r})\rho (\stackrel{}{r})=_i\varphi (\stackrel{}{r}_i)`$. $`\rho _\varphi `$ is the inner product of $`\varphi `$ and $`\rho (\stackrel{}{r})`$, and we can regard it as a projection of the density onto a basis function $`\varphi (r)`$, i.e., $`\rho _\varphi =<\rho |\varphi >`$. It weights the density fluctuations according to their spatial position. The application of the external potential will induce an average change $`\delta \rho _\varphi `$ in $`\rho _\varphi `$:

$$\delta \rho _\varphi =\rho _\varphi _U\rho _\varphi _{U=0}$$
(2)
where the thermal average $`\rho _\varphi _U`$ is given by

$$\rho _\varphi _U=\frac{1}{𝒵}\text{Tr}\left[e^{\beta (H+U)}\rho _\varphi \right]$$
(3)
The partition function $`𝒵=\text{Tr}e^{\beta (H+U)}`$ and $`\beta `$ is the inverse temperature. For small values of $`\mathrm{\Delta }P`$, this change can be calculated using perturbation theory. Up to third order in $`\mathrm{\Delta }P`$, we find

$`\delta \rho _\varphi `$ $`=`$ $`{\displaystyle \frac{\beta \mathrm{\Delta }P}{\rho _o}}\rho _\varphi ^2_c+{\displaystyle \frac{\beta ^2\mathrm{\Delta }P^2}{2\rho _o^2}}\rho _\varphi ^3_c`$ (5)
$`{\displaystyle \frac{\beta ^3\mathrm{\Delta }P^3}{6\rho _o^3}}\rho _\varphi ^4_c,`$
where the cumulant averages are

$`\rho _\varphi ^2_c`$ $`=`$ $`\rho _\varphi ^2\rho _\varphi ^2`$ (6)
$`\rho _\varphi ^3_c`$ $`=`$ $`\rho _\varphi ^33\rho _\varphi \rho _\varphi ^2+2\rho _\varphi ^3`$ (7)
$`\rho _\varphi ^4_c`$ $`=`$ $`\rho _\varphi ^44\rho _\varphi \rho _\varphi ^33\rho _\varphi ^2^2+`$ (9)
$`12\rho _\varphi ^2\rho _\varphi ^26\rho _\varphi ^4`$
with the thermal average $`\rho _\varphi ^n=\rho _\varphi ^n_{U=0}`$. The third order cumulant, eq.(7), is zero in the liquid phase because for every configuration there exists an equivalent configuration with the opposite sign of $`\rho _\varphi \rho _\varphi `$ and so we will not consider this term any further. We can recast eq. (5) as a power series in the perturbation $`\mathrm{\Delta }P`$:

$$\frac{\delta \rho _\varphi }{N}=\frac{1}{6\rho _ok_BT}\chi _l\mathrm{\Delta }P+\frac{1}{6(\rho _ok_BT)^3}\chi _{nl}(\mathrm{\Delta }P)^3$$
(10)
where

$$\chi _l=\frac{6}{N}(\rho _\varphi )^2_c\chi _{nl}=\frac{1}{N}(\rho _\varphi )^4_c.$$
(11)
where $`k_B`$ is Boltzmann’s constant. In the remainder of this paper we will focus our attention on the linear ($`\chi _l`$) and nonlinear ($`\chi _{nl}`$) dimensionless generalized compressibilities defined by the above expressions. We now discuss the choice of the function $`\varphi `$. We consider applying the potential along the direction $`\mu `$ of one of the coordinate axes so that $`\varphi (\stackrel{}{r})=\varphi (r^\mu )`$. A natural candidate for $`\varphi (r^\mu )`$ is $`\mathrm{cos}(k_\mu r^\mu )`$ with $`k=2\pi n/L`$, where $`n=1,2,\mathrm{}`$ and $`L^3`$ is the volume $`V`$. In this case, $`\rho _\varphi `$ is the $`k^{th}`$ mode of the cosine transform of the density. However, it is sufficient to consider the simpler function $`\varphi (r^\mu )=|r^\mu |/L`$. The absolute value means that all the particles feel a force along the $`\mu `$th direction pointing towards the origin. The results are very similar to $`\varphi (r^\mu )=\mathrm{cos}(k_\mu r^\mu )`$ for small $`k`$ at a fraction of the computational cost. So our results in this paper correspond to $`\rho _\varphi =_i|r_i^\mu |/L`$. This is rather like a center of mass. Since the system is isotropic, we average over the direction $`\mu `$.

We have performed a molecular dynamics simulation on a glass forming liquid  consisting of a 50:50 binary mixture of soft spheres in three dimensions. The two types of spheres, labelled A and B, differ only in their sizes. The interaction between two particles a distance $`r`$ apart is given by $`V_{\alpha \beta }(r)=ϵ[(\sigma _{\alpha \beta }/r)^{12}+X_{\alpha \beta }(r)]`$ where the interaction length $`\sigma _{\alpha \beta }=(\sigma _\alpha +\sigma _\beta )/2`$, $`\sigma _B/\sigma _A=1.4`$ ($`\alpha `$, $`\beta =`$ A, B). For numerical efficiency, we set the cutoff function $`X_{\alpha \beta }(r)=r/\sigma _{\alpha \beta }\lambda `$ with $`\lambda =13/12^{12/13}`$. The interaction is cutoff at the minimum of the potential $`V_{\alpha \beta }(r)`$. Energy and length are measured in units of $`ϵ`$ and $`\sigma _A`$, respectively. Temperature is given in units of $`ϵ/k_B`$, and time is in units of $`\sigma _A\sqrt{m/ϵ}`$ where $`m`$, the mass of the particles, is set to unity. The equations of motion were integrated using the leapfrog method  with a time step of 0.005. During each run the temperature was kept constant using a constraint algorithm . $`N=N_A+N_B`$ is the total number of particles. The system occupies a cube with dimensions ($`\pm `$ L/2, $`\pm `$ L/2, $`\pm `$ L/2) and periodic boundary conditions. Since $`N`$ and $`L`$ are fixed in any given run, the density $`\rho _o=N/L^3`$ is also fixed.

As a point of reference we determined the mode coupling $`T_C`$ by fitting the data for the relaxation time $`\tau (T)`$ which is defined as the time when the self part of the intermediate scattering function $`F_{s,\alpha }(\stackrel{}{k},t)`$ falls to $`1/e`$ .

$$F_{s,\alpha }(\stackrel{}{k},t)=\frac{1}{N_\alpha }\underset{i=1}{\overset{N_\alpha }{}}e^{i\stackrel{}{k}(\stackrel{}{r}_i(t)\stackrel{}{r}_i(0))}$$
(12)
where the subscript $`\alpha `$ refers to the particle type, A or B. We choose B particles. $`\stackrel{}{r}_i(t)`$ is the position of particle $`i`$ at time $`t`$, and $`\mathrm{}`$ refers to an average over different configurations. The wavevector $`\stackrel{}{k}=2\pi \stackrel{}{q}/L`$ where $`\stackrel{}{q}`$ is a vector of integers. For an isotropic system $`F_{s,\alpha }(\stackrel{}{k},t)`$ depends only on the magnitude $`k=|\stackrel{}{k}|`$. We choose $`k=k_{max}=2\pi q_{max}/L`$ where $`k_{max}`$ is the position of the first maximum of the partial static structure factor $`S(k)`$. For B type particles we use $`q_{max}=8.3666`$. We fit the ideal MCT form  $`\tau (T)=A^{}(TT_C)^\gamma `$ with $`A^{}=69`$, $`\gamma =1.6`$, $`T_C=0.306\pm 0.005`$. In finding $`T_C`$, we used data from 7 temperatures in the range $`0.33<T0.39`$ below the caging temperature ($`T_{\mathrm{cage}}=0.4`$) where the intermediate scattering function first begins to show a plateau. Our fit is consistent with the mode coupling theory requirement that $`\gamma 1.6`$.

Our procedure for doing runs is as follows. We start each run at a high temperature (T=1.5) and lower the temperature in steps of $`\mathrm{\Delta }T=0.05`$. At each temperature we equilibrate for $`10^4`$ molecular dynamics (md) steps and then measure the quantities of interest for $`N_\tau `$ additional steps where $`N_\tau =10^5`$, $`2\times 10^5`$, $`10^6`$, $`3\times 10^6`$ or $`10^7`$. All the particles move at each md step. The results are then averaged over up to 40 different initial conditions (different initial positions and velocities of the spheres).

As a check on our procedure for measuring $`\chi _l`$ and $`\chi _{nl}`$, we consider first the case of the crystal. To this end, we consider a system of 512 identical ($`\sigma _A=\sigma _B`$) particles at a density $`\rho _o=1.1`$. Figure 1 shows the linear compressibility and the diffusion constant as a function of temperature for the single component liquid. At high temperatures $`\chi _l`$ has a small slope which becomes steeper at low temperature. The salient feature is the very sharp drop around T=0.57 of the linear generalized compressibility and the diffusion constant. The specific heat, which is not shown, has a sharp delta function–like peak around T=0.57. The low temperature phase (T$`<`$0.57) is a crystal with sharp Bragg peaks in the structure factor. Upon heating and cooling, the transition shows hysteresis. All these observations are consistent with the fact that crystallization is a first order transition. Not shown here is the nonlinear compressibility which is zero within our numerical error.

We now examine the response of the two component glass forming liquid. All the runs were done at a density of $`\rho _o=0.6`$ and $`\sigma _B/\sigma _A=1.4`$. For these parameters crystallization is avoided upon cooling.

Figure 2 shows the linear generalized compressibility as a function of temperature for different run times. The compressibility at high temperatures is independent of $`T`$ and about an order of magnitude larger than that of the single component fluid. In the vicinity of the glass transition $`\chi _l`$ drops. Notice that as the measuring time $`t_M`$ increases, the temperature of the drop decreases and becomes more abrupt. The linear compressibility is proportional to the width of the distribution of $`\rho _\varphi `$. If we regard $`\rho _\varphi `$ as a generalized center of mass, then the drop in $`\chi _l`$ corresponds to the sudden narrowing of the distribution $`P(\rho _\varphi )`$ and the sudden arrest in the fluctuations of $`\rho _\varphi `$. This behavior can be quantified using a scaling ansatz: $`\chi _l(t_M,T)=g(\mu =t_M/\tau (T))`$, where the characteristic time has the Vogel–Fulcher form $`\tau (T)=\mathrm{exp}(A/(TT_o))`$. The inset of Figure 2 shows that the data collapse onto a single curve with $`A=0.75`$, $`T_o=0.15`$. (The data could not be fitted using $`\tau (T)=\tau _{MCT}(T)=A^{}(TT_C)^\gamma `$ as suggested by simple mode coupling theories .) This value of $`T_o`$ lies below the mode coupling $`T_C=0.306`$, the upper bound of the glass transition temperature. Notice in Fig. 2 that the drop associated with 107 md steps occurs approximately at the mode coupling $`T_C`$, indicating that our long runs were in equilibrium down to the mode coupling transition temperature. This scaling suggests that $`\chi _l`$ becomes a step function for infinite $`t_M`$ and that the drop in the compressibility would become a discontinuity at infinitely long times. This is consistent with a sudden arrest of the motion of the particles in the liquid which is the kinetic view of the glass transition. The abrupt drop also appears to be in agreement with Mezard and Parisi’s proposal that the glass transition is a first order phase transition . However, the drop indicates that our simulations are falling out of equilibrium and therefore we cannot really tell if there is a true thermodynamic transition.

As an independent check of the glass transition temperature, we have calculated the specific heat which is shown in Figure 3. Notice that the temperature of the peak in the specific heat agrees with the temperature at which $`\chi _l`$ drops.

The behavior seen in Figure 2 is similar to that seen in measurements of the real part of the frequency dependent dielectric function $`\epsilon ^{}(\omega )`$ . In that case as the frequency decreased, the temperature of the peak in $`\epsilon ^{}(\omega )`$ decreased and the drop in $`\epsilon ^{}(\omega )`$ below the peak became more abrupt. By extrapolating their data to $`\omega =0`$, Menon and Nagel  argued that $`\epsilon ^{}(\omega =0)`$ should diverge at the glass transition, signaling a second order phase transition. We have looked for evidence of this divergence by examining samples of different sizes to see if the linear generalized compressibility increased systematically with system size. As shown in Figure 4 we find no size dependence and no indication of a diverging linear generalized compressibility. We also find no size dependence for the specific heat and hence, no evidence of a diverging specific heat (not shown). This is corroborated by recent MCT calculations of a molten salt which find that $`\epsilon ^{}(\omega 0)`$ goes to a finite value as the glass transition is approached .

We have found hysteresis at the glass transition by first cooling a system of 512 particles to our lowest temperature $`T=0.1`$ and then heating in steps of $`\mathrm{\Delta }T=0.05`$. As before we equilibrate at each temperature for $`10^4`$ time steps and then measure quantities for an additional $`10^6`$ time steps. Our results are shown in the inset of Figure 4. Notice the slight hysteresis with the rise in $`\chi _l`$ upon warming being at a slightly higher temperature than the drop in $`\chi _l`$ upon cooling. This hysteresis is consistent with the kinetic arrest of motion.

We now turn to the case of the nonlinear generalized compressibility $`\chi _{nl}`$ given by eq. (11). We are motivated by the case of spin glasses where the nonlinear magnetic compressibility diverges at the spin glass transition while the linear compressibility only has a cusp . There have been a few studies of nonlinear response functions in real glasses , but these have not found any divergences. Our results are consistent with this conclusion. In particular we find that the nonlinear generalized compressibility is zero above and below the glass transition temperature, though it does show a glitch at the glass transition. There is no systematic increase with system size, indicating the absence of a divergence. Because $`\chi _{nl}`$ is sensitive to the tails of the distribution of $`\rho _\varphi `$, one must be careful to obtain a good ensemble average in the liquid above the glass transition temperature. We have done this by doing 32 runs, each involving 200,000 time steps, with different initial conditions, stringing them together as though they were one long run and then taking the appropriate averages. This produces a better ensemble average of $`<\rho _\varphi ^2>^2`$ which enters into $`\chi _{nl}`$ in eq. (7). $`\chi _{nl}`$ also took longer to equilibrate than $`\chi _l`$. A plot $`\chi _{nl}`$ versus run time shows that one needs to run at least $`10^6`$ time steps before $`\chi _{nl}`$ appeared to saturate.

To summarize, we have introduced a new thermodynamic quantity which depends solely on the positions of the particles and not on their histories. This quantity drops abruptly at the glass transition which is compatible with a kinetic arrest of motion, but not with an underlying second order phase transition. This generalized compressibility can be measured experimentally. It can be directly measured in colloidal experiments which monitor the positions of the particles . Measurements of the width of the distribution of $`\rho _q`$, the spatial Fourier transform of the density, would also give the linear generalized compressibility.

We thank Andrea Liu and Sharon Glotzer for helpful discussions. This work was supported in part by CULAR funds provided by the University of California for the conduct of discretionary research by Los Alamos National Laboratory.


text:
# 1 INTRODUCTION

## 1 INTRODUCTION

Gamma-ray bursts (GRBs) are a transient astrophysical phenomenon in which the emission is confined exclusively to high energies: they are detected at gamma-ray bands and last shortly (from a few milliseconds to several hundred seconds). Since the discovery of the events about thirty years ago (Klebesadel et al. 1973), more and more data have been obtained (e.g., Fishman et al. 1994; Meegan et al. 1994, 1996, 1998; Paciesas et al. 1997, 1999). Based on the data available so far, statistical studies have recently unveiled many basic properties of the objects. For example, two separate classes of GRBs with a division at the duration of 2 s were detected (Dezalay et al. 1992; Hurley 1992; Kouveliotou et al. 1993); a correlation between the peak spectral hardness and the peak intensity was discovered (Atteia et al. 1991, 1994, 1997; Belli 1992; Mitrofanov et al. 1992, 1996; Paciesas et al 1992; Nemiroff et al. 1994; Mallozzi et al. 1995; Dezalay et al. 1997).

More recently, the hardness-duration correlation was confirmed by Fishman (1999) with a large number of bursts observed with BATSE. Generally, the hardness ratio is defined as the fluence in channel 3 ($`100`$ to $`300`$ keV) divided by the fluence in channel 2 ($`50`$ to $`100`$ keV) (see, e.g., Fishman et al. 1994). This quantity is more likely to present an intrinsic property of the objects rather than an observational one. It is noticed that, besides the two channel fluences, there are two other channel fluences available in the BATSE burst catalogs, e.g., channels 1 ($`20`$ to $`50`$ keV) and 4 ($`>300`$ keV) (see, e.g., Paciesas et al. 1997). Therefore, one can define more than one hardness ratios. We wonder if there are any meaningful statistical relations between these different hardness ratios. We also want to make clear if there are any differences between the two classes of the objects in the distributions of the ratios.

In Section 2, correlations between various hardness ratios of GRBs will be studied, while in Section 3, differences between the two classes of the objects in the distributions of the ratios will be investigated.

## 2 CORRELATIONS BETWEEN HARDNESS RATIOS

In the 4B catalog (Paciesas et al. 1997), there are 1027 sources with all of their four channel (1-4) fluences available. We then get a rather sizable sample from the catalog (called Sample 1). A hardness ratio $`HR_{ij}`$ is defined as

$$HR_{ij}f_i/f_j,$$
(1)
where $`i`$ and $`j`$ denote the channels $`i`$ and $`j`$ respectively, and $`f_i`$ and $`f_j`$ represent the fluences of the corresponding channels respectively.

The plots of $`\mathrm{log}HR_{ij}\mathrm{log}HR_{i^{}j^{}}`$ for Sample 1 are shown in Figures 1.1 to 1.15. The correlation coefficients for every couple of hardness ratios are presented in Table 1.

Table 1. Correlation coefficients for Sample 1 (N=1027)

|  | $`\mathrm{log}HR_{31}`$ | $`\mathrm{log}HR_{32}`$ | $`\mathrm{log}HR_{41}`$ | $`\mathrm{log}HR_{42}`$ | $`\mathrm{log}HR_{43}`$ |
| --- | --- | --- | --- | --- | --- |
| $`\mathrm{log}HR_{21}`$ | 0.746 | 0.276 | 0.373 | 0.080 | -0.037 |
| $`\mathrm{log}HR_{31}`$ |  | 0.846 | 0.677 | 0.488 | 0.186 |
| $`\mathrm{log}HR_{32}`$ |  |  | 0.678 | 0.640 | 0.298 |
| $`\mathrm{log}HR_{41}`$ |  |  |  | 0.955 | 0.849 |
| $`\mathrm{log}HR_{42}`$ |  |  |  |  | 0.924 |

From Table 1 and Figures 1.1 to 1.15 we find that almost all hardness ratios are mutually correlated, with only $`\mathrm{log}HR_{21}`$ and $`\mathrm{log}HR_{43}`$ being uncorrelated. The most obvious correlation is that between $`\mathrm{log}HR_{41}`$ and $`\mathrm{log}HR_{42}`$. It suggests that, statistically, the slope of the higher part of the spectrum of most GRBs has nothing to do with that of the lower part; emissions at higher energy bands from the objects must be significantly different for different sources, while radiations at lower energy bands must be relatively similar.

## 3 DIFFERENCES BETWEEN THE TWO CLASSES

To study the differences between the two classes of GRBs in the distributions of hardness ratios, we simply concern those sources with their data of the duration as well as the four channel fluences available. The bursts are divided into two classes with a division at the duration of 2 s: those sources with $`T_{90}2`$ s belong to the short duration class and those with $`T_{90}>2`$ s constitute the long duration class, where $`T_{90}`$ is the time during which the burst integrated counts increases from $`5\%`$ to $`95\%`$ of the total counts. We find 249 short duration sources (called Sample 2) and 706 long duration bursts (called Sample 3) from the 4B catalog.

The plots of $`\mathrm{log}HR_{ij}\mathrm{log}HR_{i^{}j^{}}`$ for Samples 2 and 3 are shown in Figures 2.1 to 2.15. The correlation coefficients for every couple of hardness ratios for the samples are presented in Tables 2 and 3.

Table 2. Correlation coefficients for Sample 2 (N=249)

|  | $`\mathrm{log}HR_{31}`$ | $`\mathrm{log}HR_{32}`$ | $`\mathrm{log}HR_{41}`$ | $`\mathrm{log}HR_{42}`$ | $`\mathrm{log}HR_{43}`$ |
| --- | --- | --- | --- | --- | --- |
| $`\mathrm{log}HR_{21}`$ | 0.653 | -0.050 | 0.206 | -0.218 | -0.235 |
| $`\mathrm{log}HR_{31}`$ |  | 0.724 | 0.558 | 0.280 | -0.068 |
| $`\mathrm{log}HR_{32}`$ |  |  | 0.549 | 0.568 | 0.125 |
| $`\mathrm{log}HR_{41}`$ |  |  |  | 0.910 | 0.790 |
| $`\mathrm{log}HR_{42}`$ |  |  |  |  | 0.887 |

Table 3. Correlation coefficients for Sample 3 (N=706)

|  | $`\mathrm{log}HR_{31}`$ | $`\mathrm{log}HR_{32}`$ | $`\mathrm{log}HR_{41}`$ | $`\mathrm{log}HR_{42}`$ | $`\mathrm{log}HR_{43}`$ |
| --- | --- | --- | --- | --- | --- |
| $`\mathrm{log}HR_{21}`$ | 0.757 | 0.308 | 0.324 | 0.059 | -0.060 |
| $`\mathrm{log}HR_{31}`$ |  | 0.855 | 0.606 | 0.425 | 0.142 |
| $`\mathrm{log}HR_{32}`$ |  |  | 0.624 | 0.572 | 0.254 |
| $`\mathrm{log}HR_{41}`$ |  |  |  | 0.964 | 0.874 |
| $`\mathrm{log}HR_{42}`$ |  |  |  |  | 0.939 |

From Tables 2 and 3 we find that, the two classes have similar correlations between hardness ratios: many of them are mutually correlated, with only a few exceptions. For the short duration class, $`\mathrm{log}HR_{21}`$ and $`\mathrm{log}HR_{32}`$ as well as $`\mathrm{log}HR_{31}`$ and $`\mathrm{log}HR_{43}`$ are not correlated, while for the long duration class, $`\mathrm{log}HR_{21}`$ and $`\mathrm{log}HR_{42}`$ as well as $`\mathrm{log}HR_{21}`$ and $`\mathrm{log}HR_{43}`$ are not correlated, suggesting that, statistically, the slope of the higher part of the spectrum of the long duration bursts has nothing to do with that of the lower part. The most obvious correlation for both the short and long duration classes is that between $`\mathrm{log}HR_{41}`$ and $`\mathrm{log}HR_{42}`$, showing that, statistically, emissions at higher energy bands from the bursts of both long and short duration classes must be significantly different for different sources, while radiations at lower energy bands must be relatively similar.

For the most obvious correlations shown in the two tables, those for the long duration class are more significant than that for the short duration class. One can find from Figures 2.1 to 2.15 that the data of hardness ratios for the long duration class are more concentrated, while that for the short duration class are more scattered. This may be the most significant difference between the two classes shown in the figures. It implies statistically that the profiles of the spectra between the long duration bursts must be more similar than that between the short duration bursts. The long duration bursts must share more common properties than the short duration bursts.

Figures 2.1 to 2.15 also show different distributions of the two classes. The short duration bursts tend to occupy the right top part of the figures, showing that, in general, their spectra must be harder than that of the long duration bursts, confirming what the well-known hardness-duration correlation reveals (see, e.g., Fishman 1999). This can be verified by a direct calculation for the two classes. For example, the mean value of $`\mathrm{log}HR_{41}`$ for Sample 2 is 1.5 while that for Sample 3 is 0.67.

## 4 DISCUSSION AND CONCLUSIONS

In the last two sections, we study the correlations between various hardness ratios of GRBs and investigate if there are any differences between the two classes of the objects in the distributions of the ratios with the data of the 4B catalog. The results are: (1) almost all hardness ratios of GRBs are mutually correlated; (2) the most obvious correlation for both the short and long duration classes is that between $`\mathrm{log}HR_{41}`$ and $`\mathrm{log}HR_{42}`$; (3) for the short duration class, $`\mathrm{log}HR_{21}`$ and $`\mathrm{log}HR_{32}`$ as well as $`\mathrm{log}HR_{31}`$ and $`\mathrm{log}HR_{43}`$ are not correlated, while for the long duration class, $`\mathrm{log}HR_{21}`$ and $`\mathrm{log}HR_{42}`$ as well as $`\mathrm{log}HR_{21}`$ and $`\mathrm{log}HR_{43}`$ are not correlated; (4) those most significant correlations for the long duration class are more obvious than that for the short duration class; (5) the data of hardness ratios for the long duration class are more concentrated, while that for the short duration class are more scattered; (6) the sources of the two classes are distributed in different areas in the $`\mathrm{log}HR_{ij}\mathrm{log}HR_{i^{}j^{}}`$ plot. These results suggest that, statistically, the slope of the higher part of the spectrum of the long duration bursts has nothing to do with that of the lower part; emissions at higher energy bands from the bursts of both short and long duration classes must be significantly different for different sources, while radiations at lower energy bands must be relatively similar; the spectrum of the short duration bursts must be harder than that of the long duration bursts, confirming what the well-known hardness-duration correlation reveals; the profiles of the spectra between the long duration bursts must be more similar than that between the short duration bursts. The long duration bursts must share more common properties than the short duration bursts.

According to the most successful model proposed so far, GRBs would result in relativistically expanding fireballs, and the gamma-ray emission would arise after the fireball becomes optically thin, in shocks occurring either because the ejecta run into an external medium or because internal shocks occur in a relativistic wind (see Mészáros & Rees 1997). Several of the features reported in the first GRB detected over timescales of greater than days at X-ray and optical wavelengths, i.e., GB 970228 (Costa et al. 1997), agreed well with the theoretical expectations of the model. In the well-known relativistic beaming model used in AGNs (Rees 1966, 1967), the enhancement of the observed flux for a relativistically moving component is $`\delta ^p`$, with $`p=3+\alpha `$ in the case of a moving sphere and $`p=2+\alpha `$ in the case of a continuous jet, where $`\delta `$ is the Doppler factor and $`\alpha `$ is the spectral index of the component (see, e.g., Ghisellini et al. 1993). Recently, the correlation between polarization and variation for BL Lac objects was found able to be adopted to GRBs and a value as big as 100 of the Doppler factor for the bursts was estimated by the relation (see Fan et al. 1997, 1999; Cheng et al. 1999). As mentioned above, different to that radiated at lower energy bands, emissions of GRBs at higher energy bands are significantly different for different sources. This leads to the prospect that the Doppler boosting for GRBs might play a more important role at higher energy bands than at lower bands. If it is true, higher energies must be emitted by the fastest moving ejecta while lower energies be radiated by the slower moving components in both the isotropic expanding fireball model or the beaming model; or, the range of the spectral index of GRBs at higher energy bands must be significantly wider than that at lower energy bands.

ACKNOWLEDGMENTS

This work was supported by the United Laboratory of Optical Astronomy, CAS, the Natural Science Foundation of China, and the Natural Science Foundation of Yunnan.

REFERENCES

> Atteia, J.-L., et al. 1991, Proc. 22d Int. Cosmic-Ray Conf. (Dublin), 93
> Atteia, J.-L., et al. 1994, A&A, 288, 213
> Atteia, J.-L., et al. 1997, in AIP Conf. Proc. 384, ed. C. Kouveliotou, M. F. Briggs, & G. J. Fishman (New York: AIP), 301
> Belli, B. M. 1992, in AIP Conf. Proc. 265, ed. W. S. Paciesas & G. J. Fishman (New York: AIP), 100
> Cheng, K. S., et al. 1999, astro-ph/9910543
> Costa, E., et al. 1997, Nature, 387, 783
> Dezalay, J.-P., et al. 1992, in AIP Conf. Proc. 265, ed. W. S. Paciesas & G. J. Fishman (New York: AIP), 304
> Dezalay, J.-P., et al. 1997, ApJ, 490, L17
> Fan, J. H., et al. 1997, A&A 327, 947
> Fan, J. H., et al. 1999, astro-ph/9910540
> Fishman, G. J., 1999, A&AS, 138, 395
> Fishman, G. J., et al. 1994, ApJS, 92, 229
> Ghisellini, G., et al. 1993, ApJ, 407, 65
> Hurley, K. C. 1992, in AIP Conf. Proc. 265, ed. W. S. Paciesas & G. J. Fishman (New York: AIP), 3
> Klebesadel, R., Strong, I., Olson, R. 1973, ApJ, 182, L85
> Kouveliotou, C., et al. 1993, ApJ, 413, L101
> Mallozzi, R. S., et al. 1995, ApJ, 454, 597
> Meegan, C. A., et al. 1994, The Second BATSE Burst Catalog, available electronically from the Compton Observatory Science Support Center
> Meegan, C. A., et al. 1996, ApJS, 106, 65
> Meegan, C. A., et al. 1998, in AIP Conf. Proc. 428, Gamma-Ray Bursts: 4th Huntsville Symp., ed. C. A. Meegan, R. D. Preece, & T. M. Koshut (New York: AIP), 3
> Mészáros, P., & Rees, M. J. 1997, ApJ, 476, 232
> Mitrofanov, L., et al. 1992, in Gamma-ray Bursts, ed. C. Ho, R. I. Epstein, & E. E. Fenimore (Cambridge: Cambridge Univ. Press), 203
> Mitrofanov, L., et al. 1996, ApJ, 459, 570
> Nemiroff, R. J., et al. 1994, ApJ, 435, L133
> Paciesas, W. S., et al. 1992, in AIP Conf. Proc. 265, ed. W. S. Paciesas & G. J. Fishman (New York: AIP), 190
> Paciesas, W. S., et al. 1997, The Fourth BATSE Burst Catalog, available electronically at http://cossc.gsfc.nasa.gov/cossc/batse/4Bcatalog
> Paciesas, W. S., et al. 1999, ApJS, 122, 465
> Rees, M. J. 1966, Nature, 211, 468
> Rees, M. J. 1967, MNRAS, 135, 345

FIGURE CAPTION

> Figures 1.1 — 1.15. The plot of $`\mathrm{log}HR_{ij}\mathrm{log}HR_{i^{}j^{}}`$ for Sample 1.
>
> Figures 2.1 — 2.15. The plot of $`\mathrm{log}HR_{ij}\mathrm{log}HR_{i^{}j^{}}`$ for Samples 2 and 3, where an open circle represents a source of Sample 2 and a filled circle stands for a source of Sample 3.