text:
# Infrared Signatures of Protoplanetary Disk Evolution

## 1. Introduction

Protoplanetary disks are primarily detected through their signature infrared excess emission relative to a stellar photosphere (e.g., Mendoza 1968; Rydgren, Strom, & Strom 1976; Cohen & Kuhi 1979; Rucinski 1985). The IR excesses arise from the reprocessing of stellar radiation and the liberation of accretion luminosity within the disk (e.g., Lynden-Bell & Pringle 1974; Adams & Shu 1986; Adams, Lada, & Shu 1987; Kenyon & Hartmann 1987; Adams, Emerson, & Fuller 1990; Lada & Adams 1992; Hillenbrand et al. 1992). At other wavelengths, disk masses and velocity structures have been measured with radio and mm interferometers (Beckwith et al. 1990; Beckwith & Sargent 1993; Koerner, Sargent, & Beckwith 1993; Koerner & Sargent 1995; Dutery et al. 1996; Wilner & Lay 2000) and more recently, high resolution HST observations and ground based adaptive optics and speckle imaging techniques have imaged disks via scattered light (O’Dell, Wen, & Hu 1993; Burrows et al. 1996; Stapelfeldt et al. 1998; Padgett et al. 1999; Krist et al. 2000; Cotera et al. 2001; Grady et al. 2000; Roddier et al. 1996; Koresko 1998).

The disks that are detected range in mass, structure, and formation mechanisms. Massive, optically thick, flaring circumstellar disks, $`M_{\mathrm{disk}}10^3M_{}`$, are a natural product of the star formation process (e.g., Shu, Adams, & Lizano 1987) and are detected around pre-main sequence Classical T-Tauri stars. Low mass “debris disks” ($`M_{\mathrm{disk}}10^5M_{}`$) are detected around older main sequence stars (Backman & Paresce 1993). The dust that gives rise to the IR excesses and scattered light images of debris disks is believed to form from collisions of planetessimals within the disks.

The advent of large IR detector arrays now allows for simultaneous multiwavelength observations of large numbers of disks in star clusters (e.g., Haisch, Lada, & Lada 2000). The upcoming SIRTF mission will greatly increase the wavelength coverage for studying clusters, thereby enabling detailed studies of disks over a wide range in age and mass. We therefore need to know what spectral indicators are best suited for detecting disks and determining their mass evolution.

In this paper our working hypothesis is that through time the disk mass decreases but the disk structure and stellar luminosity do not change. In our models the dominant opacity source for $`\lambda 100\mu `$m is small ($`50\mu \mathrm{m}`$) grains, so another way of stating our model is that we assume the small particle mass traces the total mass of gas + dust. The small grain mass may decrease by the accretion of disk material or coagulation and growth of large grains and rocks while the gas is either depleted or accreted onto the star keeping the gas/dust ratio constant. The recent detection of $`H_2`$ in $`\beta `$ Pictoris’ disk and the derived gas/dust $`100`$ is consistent with our assumption (Thi et al. 2001).

We investigate the observational signatures of this simple model for disk evolution and present radiative equilibrium models for a range of disk masses and viewing angles. Model spectra are compared with semi-analytic approximations for optically thick flat and flared disks and predictions are made as to what spectral regions are most sensitive to the evolution of disk mass and at what wavelengths disk detections are robust. We perform simulations for a fiducial disk structure and dust size distribution that fits the SED of the edge-on disk system HH30 IRS. (Wood et al. 2001). § 2 describes the ingredients of our models, § 3 presents our SED simulations and shows how SIRTF colors are sensitive to disk mass. § 4 presents models that explore variations in parameters other than disk mass and discusses alternative disk evolution models. § 5 compares our models with color-color diagrams of Taurus-Auriga sources and we summarize our findings in § 6.

## 2. Models

### 2.1. Radiative Equilibrium Calculation

Theoretical temperature distributions and emergent spectra have been calculated for T Tauri disks using a variety of techniques. These include approximations for opaque flat disks (Adams & Shu 1986; Adams, Lada, & Shu 1987; Lada & Adams 1992) and extensions of this approach to flared disks (Kenyon & Hartmann 1987; Chiang & Goldreich 1997, 1999); two dimensional radiation transfer techniques (Efstathiou & Rowan-Robinson 1991); diffusion approximations (e.g., Sonnhalter, Preibisch, & Yorke 1995; Boss & Yorke 1996); disks which are modeled as “spherical sectors” (Men’shchikov & Henning 1997); vertical structure calculations performed by dividing the disk into plane parallel annuli (Calvet et al. 1992; Bell et al. 1997; Bell 1999; D’Alessio et al. 1998, 1999); and Monte Carlo techniques (Wolf, Henning, & Stecklum 1999). Common approximations in many of the “traditional” codes are to conduct the radiation transfer in only one direction in the disk and assume the dust grains scatter radiation isotropically. Also many techniques are limited to the study of optically thick and hence massive disks.

Monte Carlo techniques are straightforward to adapt to any geometry, mass, and can accurately include polarization and multiple, anisotropic scattering. We use the Monte Carlo radiative equilibrium technique of Bjorkman & Wood (2001) which conserves energy exactly. For very optically thick disks we use our Monte Carlo technique for the upper layers of the disk (i.e., the disk “atmosphere”) and a diffusion approximation for the densest midplane regions. For the most massive disks we simulate, typically less than 1% of the flux requires a diffusion treatment. More details of our adaptation of the Bjorkman & Wood (2001) technique for simulating T Tauri disks will be presented in Bjorkman, Wood, & Whitney (2001, in preparation). An advantage of using Monte Carlo techniques for studying a range of disk masses, is that we are not restricted to one-dimensional radiation transfer and can therefore simulate SEDs of low mass optically thin disks in which radial transport of photons is important. The output of our code is the disk temperature structure (due to heating by stellar photons and accretion luminosity) and the emergent SED and polarization spectrum at a range of viewing angles. A calculation of the hydrostatic disk structure (e.g., Chiang & Goldreich 1997; D’Alessio et al. 1999) can be included in the Monte Carlo technique, but this would require an iterative scheme. At present we have not implemented such a scheme and instead perform the radiative equilibrium calculation for a fixed disk geometry.

### 2.2. Disk Structure

The determination of disk structure from fitting SEDs and scattered light images does not yield a single structure and dust size distribution that applies to all disks. Some systems are fit with passive flat disks (e.g., Adams et al. 1987; Adams et al. 1990; Miyake & Nakagwa 1995), while others require flared disks heated by starlight and accretion luminosity (e.g., Kenyon & Hartmann 1987; Burrows et al. 1996; Stapelfeldt et al. 1998). The scattered light image of the edge-on disk of HH30 IRS (Burrows et al. 1996) has for the first time allowed the vertical structure of a protoplanetary disk to be studied directly. We therefore adopt the HH30 IRS disk as our fiducial model for our SED models. A fixed disk density structure which fits the scattered light images (Burrows et al. 1996; Cotera et al. 2001) and SED (Stapelfeldt & Monetti 1999; Wood et al. 2001) of HH30 IRS is

$$\rho =\rho _0\left(\frac{R_{}}{\varpi }\right)^\alpha \mathrm{exp}\frac{1}{2}[z/h(\varpi )]^2,$$
(1)
where $`\varpi `$ is the radial coordinate in the disk midplane and the scale height increases with radius, $`h=h_0\left(\varpi /R_{}\right)^\beta `$. For the HH30 IRS disk we adopt $`\beta =1.25`$, $`\alpha =2.25`$, and $`h_0=0.017R_{}`$, giving $`h(100\mathrm{A}\mathrm{U})=17`$AU.

In our simulations the inner edge of the disk is truncated at the dust destruction radius, $`R_{\mathrm{dust}}`$. Assuming $`T_{}=4000`$K, $`R_{}=2R_{}`$, and circumstellar dust sublimates at 1600K, then $`R_{\mathrm{dust}}8R_{}`$. This is larger than the dust destruction radius for optically thin dust because the reprocessed emission from the disk provides additional heating over and above the direct stellar radiation, increasing the size of the dust destruction zone (a detailed discussion of the shape of the dust destruction region will be presented in Bjorkman et al. 2001, in preparation). In currently popular magnetic accretion models the disk is truncated at a radius $`R_0`$, which may not be equal to $`R_{\mathrm{dust}}`$. If $`R_0<R_{\mathrm{dust}}`$ there will be a gas disk extending from $`R_0`$ to $`R_{\mathrm{dust}}`$ which may give rise to additional IR emission. Our models therefore assume that any material within $`R_{\mathrm{dust}}`$ is optically thin, which is a good approximation for low mass disks. For high mass disks the gas may be optically thick producing larger near-IR excesses.

### 2.3. Adopted Circumstellar Dust Properties

Recent modeling of HST images (Cotera et al. 2001) and the SED of HH30 IRS (Wood et al. 2001) indicates that the circumstellar dust size distribution extends to larger grain radii than typical ISM grains. This is in agreement with many other observations indicating grain growth within protoplanetary disks (e.g., Beckwith et al. 1990; Beckwith & Sargent 1991). This paper primarily investigates the effects of disk mass on the SED and adopts circumstellar dust properties that reproduce the HH30 IRS SED. The dust model (chemical composition, mathematical form for the size distribution, calculation of opacity and scattering parameters, etc) is described in Wood et al. (2001) and we only summarize the main features here. Specifically, we adopt, a size distribution

$$n(a)da=C_ia^p\mathrm{exp}([a/a_c]^q)da,$$
(2)
with $`p=3.5`$, $`q=0.6`$, $`a_c=50\mu \mathrm{m}`$, $`a_{\mathrm{min}}=0.01\mu \mathrm{m}`$, and $`a_{\mathrm{max}}=1`$mm. The exponential scalelength, $`a_c`$, yields dust particle sizes extending up to and in excess of $`50\mu \mathrm{m}`$. Figure 1 shows the wavelength dependence of the opacity, scattering albedo, and Heyney-Greenstein phase function asymmetry parameter (Heyney & Greenstein 1941) for this size distribution.

Recent Infrared Space Observatory spectra of Herbig Ae/Be stars (e.g., Meuss et al. 2000; Chiang et al. 2001; van den Ancker et al. 2000) are now allowing the circumstellar dust chemistry to be studied. Because we adopt the dust properties of Fig. 1 for our simulations, we have not investigated the effects of different chemical compositions on the the resulting disk SEDs. How the circumstellar chemistry effects the SED is an interesting problem, but it is beyond the scope of this paper and we present models for different disk masses and the dust properties shown in Figure 1. In addition, our models do not include additional heating from transiently heated small grains (see the new radiation transfer code of Misselt et al. 2001) which is not important in the cooler Classical T Tauri stars considered in this paper.

### 2.4. Energy Sources

The energy input to the disk is from stellar photons and accretion luminosity liberated in the disk. As discussed in the previous section, we fix the disk structure for our radiation transfer simulations. Given the disk structure (Eq. 1), $`\alpha `$ disk theory determines the accretion rate for a given disk mass. Our parameterization of the disk density and accretion follows that presented in the review by Bjorkman (1997), apart from the term $`\pm \sqrt{R_0/\varpi }`$. The accretion rate and viscosity parameter, $`\dot{M}`$ and $`\alpha _{\mathrm{disk}}`$, are related to the disk parameters by

$$\dot{M}=\sqrt{18\pi ^3}\alpha _{\mathrm{disk}}V_c\rho _0h_0^3/R_{},$$
(3)
where the critical velocity $`V_c=\sqrt{GM_{}/R_{}}`$. The flux due to viscous disk accretion, $`GM_{}\dot{M}/2R_{}`$, is generated throughout the disk midplane region according to (Shakura & Sunyaev 1973, Lynden-Bell & Pringle 1974)

$$\frac{dE}{dAdt}=\frac{3GM_{}\dot{M}}{4\pi \varpi ^3}\left[1\sqrt{\frac{R_{}}{\varpi }}\right].$$
(4)
For low mass disks, the heating due to accretion luminosity is negligible and stellar irradiation dominates the disk heating. In general we choose $`\alpha _{\mathrm{disk}}=0.01`$ (Hartmann et al. 1998), but for the most massive disk we simulate this results in a very large accretion luminosity, $`L_A>0.8L_{}`$. For this case $`\alpha _{\mathrm{disk}}`$ is adjusted so that $`L_A<0.2L_{}`$ in line with recent observational determinations of accretion luminosities in classical T Tauri stars (Hartmann et al. 1998).

## 3. Model Results: SEDs, Colors

The following models use $`T_{}=4000`$K, $`R_{}=2R_{}`$, $`R_{\mathrm{disk}}=100`$AU, and a distance of 500pc to the system.

### 3.1. SED Evolution with Disk Mass

Figure 2 shows the effect on the SED of changing the disk mass, but keeping the disk structure fixed to that of our fiducial model. The massive optically thick disks produce SEDs that resemble that of the Kenyon & Hartmann (1987) flared disk model, aside from differences (addition of scattered light to pole-on views and silicate features) due to our inclusion of a finite albedo and non-gray opacity. The most massive disk has a very large near-IR excess due to the large accretion luminosity present in this model. The model SEDs display the characteristic features present in other simulations: large infrared excess emission, flat spectrum sources at intermediate inclinations, and double peaked spectra (optical and far-IR peaks) for very large inclinations. At large inclinations, the optical peak is due to scattered starlight as the dense disk totally obscures the star and disk emission at short wavelengths.

Reducing the disk mass has the most dramatic effect at long wavelengths with the SED rapidly declining with decreasing mass. At short wavelengths the disk remains optically thick to stellar photons over a wide range of mass so the near-IR excess is not as sensitive to mass. Figure 4 presents our models again, with the three panels showing SEDs for the range of disk masses at a given inclination. These results show similar features to the models of Men’shchikov & Henning (1997, their Fig. 12) who presented SEDs for a range of optical depths in a spherical geometry with evacuated bipolar cones.

The SIRTF sensitivity limits show that at a distance of 500pc very low mass disks ($`M_{\mathrm{disk}}10^6M_{}`$) are detectable out to a wavelength of $`70\mu \mathrm{m}`$. At $`25\mu \mathrm{m}`$ SIRTF will be sensitive to photospheric flux levels at 500pc, thus allowing for the detection of even lower mass disks. Our overall conclusion from these simulations is that near-IR excesses detect disks while far-IR excesses can be used to study their mass.

### 3.2. Color Evolution with Disk Mass

Figure 4 shows the variation of colors with disk mass for $`i<60^{}`$. At long wavelengths we use the simulated flux at $`70\mu \mathrm{m}`$ and $`160\mu \mathrm{m}`$ in forming the color and have not adopted any particular color system. The colors are fairly insensitive to inclination for $`i<60^{}`$, but see §3.3 for color-color and color-magnitude diagrams that include highly inclined disks. For a passive disk with $`M_{\mathrm{disk}}=10^1M_{}`$, $`\mathrm{\Delta }(KL)0.7`$, decreasing to a fairly constant $`\mathrm{\Delta }(KL)0.4`$ for $`10^7M_{}M_{\mathrm{disk}}10^3M_{}`$. Larger $`KL`$ colors arise in massive disks where accretion luminosity is included, $`\mathrm{\Delta }(KL)1`$ for $`M_{\mathrm{disk}}=10^1M_{}`$. In our models, less massive disks do not sustain large accretion rates (see §2.3) and the $`KL`$ excess is due to reprocessing of starlight. For these disks, $`KL`$ is fairly insensitive to mass as the disk remains optically thick in the near-IR for masses $`M_{\mathrm{disk}}10^7M_{}`$. Therefore ground based near-IR observations are capable of detecting very low mass disks. Again, we emphasize that this is the mass of small ($`50\mu \mathrm{m}`$) particles which dominate the near-IR opacity.

At longer wavelengths the SED decreases with decreasing disk mass and this is reflected in the other color indexes remaining relatively flat with disk mass until the disk becomes optically thin at the waveband under study. At the longest MIPS wavelength, $`160\mu `$m, the $`K160`$ color shows a clear progression from the most massive to least massive disk.

### 3.3. Color-Color and Color-Magnitude Diagrams

Figure 5 shows color-color diagrams for our model disks assuming intrinsic stellar colors from Kenyon & Hartmann (1995, Table A5). When we include all disk inclinations there are degeneracies in the color-mass parameter space, with large viewing angles generally leading to redder colors. This result differs from other investigations (e.g., Lada & Adams 1992; Kenyon, Yi, & Hartmann 1996; Meyer, Calvet, & Hillenbrand 1997) which assumed the disk emission was proportional to $`\mathrm{cos}i`$, giving the bluest colors for edge-on viewing. More detailed radiation transfer modeling shows that a simple $`\mathrm{cos}i`$ scaling does not give the correct inclination dependence for flared disk models.

The effect of changing the disk mass is not very strong in the $`(KL)/(KN)`$ diagram for $`M_{\mathrm{disk}}10^6M_{}`$, reflecting that at short wavelengths the disks remain optically thick over a wide range of mass. The disk masses separate out more clearly in the $`(KN)/(K70)`$ diagram. Therefore color-color diagrams of large numbers of sources that compare near and far infrared colors provide a means of determining the range of disk masses within a cluster.

Figure 5 showed that inclination effects lead to a large spread in the location of different disk masses in the color-color diagrams. This leads to degeneracies in the location of different disk masses. However, Fig. 6 shows how color-magnitude diagrams can help in breaking the mass-inclination degeneracy. Placing all our disk models on color-magnitude diagrams shows that disks viewed at large inclinations occupy the lower righthand corners, i.e., highly inclined disks are faint and red.

## 4. SED and Color Dependence on Other Disk Properties

### 4.1. Flat Disks

The previous sections presented SEDs for disks which have the flaring parameters that fit HH30 IRS images and SED. However, other investigations indicate a variety of disk structures in Classical T Tauri stars (e.g., Kenyon & Hartmann 1987; Miyake & Nakagawa 1995; Chiang et al. 2001). To investigate the effects of flatter disks on SEDs and colors we repeated the simulations of Fig. 2 with a reduced scaleheight $`h_0=0.003R_{}`$ giving $`h(100\mathrm{A}\mathrm{U})=3`$AU. The smaller scaleheight results in the disk intercepting less stellar radiation (e.g., Kenyon & Hartmann 1987) and consequently smaller excesses and bluer colors than our models in § 3. For a given disk mass, the mm flux is unaltered, but the mid to far-IR emission is sensitive to the disk structure (Kenyon & Hartmann 1987; Chiang et al. 2001). The qualitative variation of the SEDs and colors with disk mass remains the same.

### 4.2. Grain Growth

In our models the dust size distribution does not change as the disk mass decreases. However, small dust grains may coagulate to form larger grains and rocks thereby altering the grain size distribution in the disk (e.g., Beckwith, Henning, & Nakagawa 2000). This will lead to changes in the wavelength dependence of the opacity. D’Alessio et al. (2001) investigated grain growth by keeping the disk mass constant and increasing the maximum grain size (effectively reducing the population of small grains) and changing the slope of the power law size distribution. For a slope $`p=3.5`$, their SED models show the largest effects occur for $`\lambda 20\mu `$m, while for $`p=2.5`$ (which puts more mass into larger particles) the SEDs are similar to ours in which the total disk mass decreases. SEDs arising from either a low mass disk or a massive disk with large grains may be distinguished with detailed SED modeling.

### 4.3. Disk Radius

Decreasing the disk radius, but keeping the mass constant results in a disk that is denser and optically thicker than our models of § 3. A smaller disk effectively removes cool material at large radii that provides the bulk of the long wavelength emission, yielding smaller far-IR and mm fluxes (e.g., Beckwith et al. 1990). The other effect of squeezing the disk into a smaller volume is an increase in the height above the midplane at which stellar photons are absorbed. This gives a geometrically thicker dust disk that can intercept more stellar photons, raising the near and mid-IR emission. The overall effects of a very small disk are to increase the wavelength at which the disk becomes optically thin yielding larger $`KL`$ colors and a slower decline in the $`K180`$ color with decreasing mass. The larger optical depth for small disks will lower the minimum disk mass that may be detected from $`KL`$ colors. Large radius disks give the opposite of these effects: less near and mid-IR, and more far-IR and mm emission.

### 4.4. Inner Disk Holes and Non-homologous Disk Evolution

In order to investigate general trends, we have presented a simple model for disk evolution of protoplanetary disks in which the mass decreases homologously with time. Real disks will be more complicated than this, with the disk radius growing, changes in the flaring parameters due to dust settling, and the possibility of accretion being terminated through the opening up of large inner disk holes and gaps, resulting in “inside-out” evolution. Current data does not provide much support for the existence of large holes in protoplanetary disks: if they do not show near-IR excesses, mid-IR and mm excesses are usually absent also (Stassun et al. 2001, Haisch, Lada & Lada 2001a), but there are some exceptions notably GM Aur (Koerner, Sargent, & Beckwith 1993). SIRTF should be able to distinguish between low mass protoplanetary disks, whose mass has evolved homologously (without creating large inner holes), and disks which have evolved by creating large inner holes and thus directly test the inside-out disk clearing scenario.

Debris disks often exhibit large inner holes or a ring-like structure (e.g., Koerner et al. 1998; Jayawardhana et al. 1998; Schneider et al. 1999). However, debris disks are not likely to be protoplanetary disks that have evolved large inner holes. Debris disks are sufficiently old that the dust in them is not remnant protoplanetary material but instead is created and continually replenished by collisions between planetessimals formed previously in the disk. Moreover, protoplanetary disks which simply evolve from the inside out by clearing large inner holes (but otherwise maintaining their initial structure), will still contain large reservoirs of material in their outer regions and thus should be considerably more massive than debris disks.

## 5. Comparison to Observations: Disk Masses in Taurus-Auriga

Observations at mm and radio wavelengths, where the disk is optically thin, provide the best probe of disk mass (e.g., Beckwith et al. 1990). When such data is not available our models show that far-IR, and in certain circumstances even mid-IR data, can provide another means of probing disk mass (if the disks are not too massive and optically thick). For example, Kenyon & Hartmann (1995) presented a compilation of optical through far-IR observations of Taurus-Auriga sources and noted a pronounced gap in the $`KN`$ distribution between the bluest Class II and the reddest Class III sources. Figure 7 shows the $`(KL)/(KN)`$ color-color diagram for the Taurus-Auriga sources along with our face-on model colors for various disk masses. In the context of our models, the gap in the $`KN`$ distribution corresponds to disks with $`M_{\mathrm{disk}}10^6M_{}`$ and indicates that there are few circumstellar disks with masses between $`10^6`$ and $`10^8`$ $`M_{}`$ in the Taurus-Auriga cloud (of course, this disk mass does not include rocks and planets). This, in turn, could be interpreted to indicate that once the disk mass falls below $`10^6M_{}`$ the timescale for clearing the remaining material in the inner ($`510`$AU) disk is very rapid.

## 6. Summary

We have investigated the observational signatures of an evolutionary model in which the disk mass decreases homologously, but the disk structure and stellar parameters remain constant with time. Our main conclusions are:

* Near-IR observations detect disks. Disks remain optically thick in the near-IR over a wide range of disk mass resulting in measurable $`KL`$ excesses for disks down to $`M_{\mathrm{disk}}10^7M_{}`$. This corresponds to a dust mass of $`M_{\mathrm{dust}}10^9M_{}`$.
* Mid-IR observations probe disk structure. Observations in the $`20\mu \mathrm{m}\lambda 100\mu \mathrm{m}`$ range are sensitive to the disk flaring parameters and are crucial for determining the degree of dust settling.
* Far-IR most sensitive to disk mass. The far-IR emission decreases rapidly with disk mass resulting in a strong correlation between $`K160`$ and disk mass. Our simulations indicate that at a distance of 500pc, SIRTF will be able to detect $`M_{\mathrm{disk}}10^6M_{}`$ with a 500s exposure at $`70\mu \mathrm{m}`$ and even lower masses at $`25\mu \mathrm{m}`$.
* General trends are reproduced for a range of disk parameters. Despite degeneracies in fitting SEDs of particular sources, we find that varying the disk parameters within generally accepted limits introduces a spread in the color-mass relationship, but does not affect the general trends.

The degeneracies inherent in fitting disk models to wavelength restricted datasets illustrates that it is necessary to have full SED coverage to constrain the disk structure and dust size distribution (e.g., Men’shchikov & Henning 1997; Chiang et al. 2001; D’Alessio et al. 2001). SIRTF, SMA, and ALMA will provide detailed SEDs for nearby sources enabling detailed modeling of the disk mass and structure. For studies of stellar clusters, such data will not be available and mm observations, which are most sensitive to disk mass, cannot as yet achieve the required sensitivity and resolution to study low mass disks and/or distant clusters. Therefore, we must appeal to statistics and determine average properties of sources in different clusters for studying disk evolution. SIRTF will provide broad band colors for hundreds of sources within a given cluster, from which we can construct median colors and the spread around the median for each cluster (e.g., Haisch, Lada, & Lada 2001b). From our models it is apparent that NIR excesses detect disks, and therefore the absence of a $`KL`$ excess within a cluster of a given age can determine the maximum disk lifetime. By comparing the median colors from different clusters with a spread in age (as a whole and as functions of position within the cluster), we can search for trends as a function of cluster age. Therefore, the far-IR probed by SIRTF may allow us to determine a timescale for disk evolution — i.e., in our homologous evolutionary model we would be able to convert the color-mass correlation to one which tracks disk mass in clusters as a function of age, analogous to the recent ISO analysis of debris disks by Spangler et al. (2001).

We acknowledge financial support from NASA’s Long Term Space Astrophysics Research Program, NAG5 6039 (KW), NAG5 8412 (BW), NAG5 7993 (MW), NAG5 3248 (JEB); the National Science Foundation, AST 9909966 (BW and KW), AST 9819928 (JEB), and a PPARC Advanced Fellowship (KW).


format: text

text:
# Non-topological solitons in brane world models

## Abstract

We examine some general properties of a certain class of scalar filed theory models containing non-topological soliton solutions in the context of brane world models with compact large extra dimensions. If a scalar field is allowed to propagate in extra space, then, beside standard Kaluza-Klein type excitations, a whole new class of very massive soliton-type states can exist. Depending on their abundance, they can be important dark matter candidates or give significant contribution to entropy and energy density in our universe.

It is well known that non-linear scalar filed theory models have a very rich structure. The scalar filed theory, with only linear couplings, admits only plane wave solutions. Wave packets can be constructed by superposition of plane waves but they are always dispersive. If non-linear couplings are present, then the theory can contain spatially confined nondispersive solutions — solitons. In renormalizable relativistic local field theories, all the soliton solutions belong to two groups — topological and non-topological solitons .

The necessary condition for the existence of a topological soliton is that there should exist degenerate vacuum states so that the boundary conditions at infinity for a soliton state are topologically different from that of a physical vacuum state. There are well known examples of such solutions in field theory — monopoles, strings, domain walls etc. Their stability is insured against decay into plain waves because of topologically non-trivial boundary conditions at infinity.

The necessary condition for the existence of a non-topological soliton is that there should exist an additive conservation law. The boundary conditions at infinity are the same as for the vacuum state. Because of the absence of non-trivial topology, the stability of non-topological solitons depends on which type of solution is of the lowest energy***It is also possible to achieve the stability by including some additional fields. For example, the stability of electroweak strings can be provided by fermionic zero modes in theory .

There is a simple result called Derrick’s theorem  that imposes severe restrictions on the number of scalar field models which can contain solitons. In principle, if the number of space-like dimensions in theory is greater than one, then the only time independent solutions of finite energy are the vacuum states which are constant everywhere and which minimize the potential energy. There are several ways to go around this result. In order to have soliton solutions, we can (1) include fields with non-zero spin (2) consider scalar models with constraints (3) consider time-dependent but non-dispersive solutions.

In literature, a great part of attention was devoted to topological solitons in arbitrary number of space-time dimensions with both compact and non-compact extra dimensions. However, analysis of non-topological soliton solutions in literature was mainly limited to calculations based on assumption that the space-time is $`(3+1)`$-dimensional on all scales and that the strength of the gravitational interaction is given by the Plank mass scale, $`M_{Pl}`$ . The goal of this paper is to present some general features of non-topological solitons in the context of recently proposed theories with large compact extra dimensions. We will show that new assumptions about the dimensionality of our space-time and the gravitational strength at small distances allows the existence of two new types of particles which can possibly play a significant role in cosmology.

Since we are mainly interested in questions of importance for cosmology, we analyze non-linear scalar field theory models in the presence of gravity. Consider gravitational Lagrangian in $`(4+d)`$ dimensions with a complex massive scalar field $`\mathrm{\Phi }`$:

$$L=d^{4+d}x\sqrt{g^{(4+d)}}\left[^{(4+d)}+g^{\mu \nu }_\mu \mathrm{\Phi }_\nu \mathrm{\Phi }+U(\mathrm{\Phi })\right]$$
(1)
where $``$ is a Ricci scalar and $`U(\mathrm{\Phi })`$ is a potential for the field $`\mathrm{\Phi }`$. The field $`\mathrm{\Phi }`$ propagates in all $`(4+d)`$ dimensions. There are two distinct classes of non-topological soliton solutions in these models. The necessary condition for both of them is that the theory is invariant under space independent phase rotations:

$$\mathrm{\Phi }\mathrm{\Phi }^{}=e^{i\alpha }\mathrm{\Phi }$$
(2)
This provides conservation of a quantity $`N`$ which we identify with the particle number. Class I comprise models where $`U(\mathrm{\Phi })`$ is such that the theory does not contain non-topological soliton solutions in the absence of the gravitational field. Class II comprise models where $`U(\mathrm{\Phi })`$ is such that the theory does contain non-topological soliton solutions in the absence of the gravitational field .

Analysis of these models has been done only for the case of $`(3+1)`$ space-time dimensionsAfter this paper appeared on the web we noticed the complementary work in .. Here, we examine models of class I and II in the context of recently proposed ”world as a brane scenario”. We consider an original  of large extra dimensions (for other types of compactification of extra dimensions see ). In this proposal, our universe is a direct product of an ordinary $`(3+1)`$-dimensional space and an extra compact $`d`$-dimensional space. All the standard model fields are trapped on the brane which represents an ordinary $`(3+1)`$-dimensional space, while gravity can propagate everywhere. This scenario allows a solution of so-called hierarchy problem by bringing down a fundamental quantum gravity energy scale — $`M_F`$ — down to the electroweak energy scale. For example, in the case of toroidal compactification of the extra space, $`M_F`$ can be as low as $`1`$TeV if the radius of compactification is $`r_c1`$mm in $`d=2`$. For an observer on the brane, the standard Newton’s law of gravity remains valid on distances larger than $`r_c`$, i.e. for $`r>r_c`$ the gravitational force, $`F_g`$, between two bodies of mass $`m_1`$ and $`m_2`$ is governed by the $`4`$-dimensional gravitational constant, $`G_4=1/M_{Pl}^2`$

$$F_g=\frac{G_4m_1m_2}{r^2}.$$
(3)
However, for $`r<r_c`$ the space is effectively $`(4+d)`$ dimensional and the Newton’s law is modified:

$$F_g=\frac{G_{4+d}m_1m_2}{r^{2+d}},$$
(4)
where the $`(4+d)`$-dimensional gravitational constant $`G_{4+d}=1/M_F^{2+d}`$. It is important to note that the strength of gravitational interaction is now determined by a fundamental energy scale, $`M_F`$, and not by the Plank scale, $`M_{Pl}`$. We will see that this fact will considerably modify the basic properties of soliton states.

For comparison, we give the basic features of $`(3+1)`$ and $`(3+1+d)`$-dimensional models in parallel. The model which gives rise to class I soliton solutions in $`(3+1)`$ dimensions is described by the Lagrangian (1) with $`d=0`$ and the potential given by

$$U(\mathrm{\Phi })=m^2\mathrm{\Phi }^{}\mathrm{\Phi },$$
(5)
where $`m`$ is the mass of the field $`\mathrm{\Phi }`$ . This model supports a classical soliton solution which is regular everywhere and zero at infinity. Its amplitude is determined by an attractive force coupling constant. In this case the amplitude is proportional to $`1/\sqrt{G_4}`$ since the only attractive force in the model is gravity. Ignoring gravity loops, a renormalizable quantum theory can be obtained in analogy with well established procedures for quantum solitons.

Because of Bose-Einstein statistics, it is possible to have very massive objects with the extremely compact characteristic dimensions. The equations of motion can be derived by minimizing the total (gravitational plus scalar field) energy while keeping the number of particles, $`N`$, fixed . The spherically symmetric solution

$$\mathrm{\Phi }=1/\sqrt{2}\sigma (r)e^{i\omega t}$$
(6)
is regular everywhere for $`\omega <m`$ (non-topological soliton can be considered as an analytical continuation of a plane wave solution for which $`\omega m`$ always) and at infinity goes to zero as

$$\sigma (r)e^{\sqrt{m^2\omega ^2}r}0.$$
(7)
Thus, the characteristic size of a class I soliton object, in the lowest energy state, is proportional to $`m^1`$. Here, we assume that $`\mathrm{\Phi }`$ is some Higgs-type scalar field, and thus for the usual range of the scalar field masses the linear size of a class I soliton is microscopic. It can be shown that, because of gravitational interaction, for large $`N`$, these states are stable.

It is important to note that very reliable order of magnitude estimates describing almost all of the peculiar properties of non-topological solitons are possible to get without solving any differential equations. To see this, examine the set of $`N`$ particles (solitons) in the same orbit in a lowest energy state (a zero-node or s-state) of wavelength $`R`$. Define the effective mass of this configuration by $`MNm`$,where $`m`$ is the mass of the free particle. The kinetic energy $`N/R`$ prevents the collapse of the object, while gravitational energy $`G_4M^2/R`$ prevents its decay. The balance between the kinetic and gravitational energies gives :

$$\frac{N}{R}\frac{G_4M^2}{R}.$$
(8)
com From here we get:

$$M\frac{1}{G_4m}\frac{M_{Pl}^2}{m}.$$
(9)
For example, if $`m300`$GeV then $`M10^9`$kg. These states are known in literature as mini-soliton stars.

Consider now class I soliton model in the context of large extra dimensions. Since the linear dimensions of this kind of non-topological solitons are smaller than the compactification radius of the extra dimensions, this state is effectively a spherically symmetric $`(4+d)`$-dimensional object. At these distances, the gravitational strength is determined by a fundamental scale $`M_F`$ rather than $`M_{Pl}`$. For simplicity, we assume that the bulk cosmological constant is zero. The equilibrium condition (8) is now

$$\frac{N}{R}\frac{G_{4+d}M^2}{R^{1+d}}.$$
(10)
From here we get:

$$M\frac{R^dM_F^{d+2}}{m}.$$
(11)
This formula is valid as long as $`R`$ is smaller than the compactification radius of the extra dimensions, $`r_c`$. However, there is an upper limit on the mass of such an object in a given state. The critical mass $`M_c`$ for the formation of a black hole can be estimated by substituting in (11) an expression for the Schwarzschild radius, $`R_s`$, of a $`(4+d)`$-dimensional black hole :

$$R_s\frac{1}{M_F}\left(\frac{M}{M_F}\right)^{\frac{1}{d+1}}.$$
(12)
We get:

$$M_c\left(\frac{M_F}{m}\right)^{d+1}M_F.$$
(13)
Note that, if we take $`d=0`$ and replace $`M_F`$ with $`M_{Pl}`$, we recover $`(3+1)`$-dimensional results. It is interesting that for the generic values of parameters in the model, say $`mM_F1`$TeV, the critical mass is $`M_c1`$TeV and $`R_s1`$Tev-1. However, if the classical analysis of these states is to be believed one need to choose parameters in the model which yield the soliton mass much higher than the fundamental quantum gravity scale $`M_F`$.

Let us now examine $`(3+1)`$-dimensional class II soliton states. A model which gives rise to this class of non-topological solitons solutions contains a complex scalar field $`\mathrm{\Phi }`$ whose action is invariant under (2). Beside this, it has to satisfy an additional condition — that in the absence of the gravitational field, the theory has non-topological soliton solutions. This can be achieved by introducing a new Hermitian scalar field $`\chi `$ to the Lagrangian (1) (with $`d=0`$) with the potential:

$$V(\chi )=\frac{1}{2}\mu ^2\chi ^2\left(1\frac{\chi }{\chi _0}\right)^2,$$
(14)
where $`\mu `$ is the mass of the field $`\chi `$. The normal vacuum state is $`\chi =0`$, while $`\chi =\chi _0`$ is the false vacuum state. Without gravity, the theory contains a non-topological soliton solution. The soliton contains an interior in which $`\chi \chi _0`$, a shell of width $`\mu ^1`$ over which $`\chi `$ changes from $`\chi _0`$ to $`0`$, and an exterior which is essentially the vacuum. The field $`\mathrm{\Phi }`$, with conserved particle number $`N`$ is confined to the interior. The only role of the field $`\chi `$ is to limit the space volume available to the filed $`\mathrm{\Phi }`$. In opposite, the only solutions for the filed $`\mathrm{\Phi }`$ in the absence of gravity would be the plain waves.

Consider a shell of radius $`R`$ with above mentioned configuration. The kinetic energy inside the shell is:

$$E_k\frac{N}{R}.$$
(15)
Since the value of $`\chi `$ from interior to exterior changes from $`\chi _0`$ to $`0`$, the shell contains a surface energy

$$E_ssR^2,$$
(16)
where the surface tension $`s\mu \chi _0^2`$. The equilibrium radius can be found by minimizing the total energy $`E=E_k+E_s`$ with respect to $`R`$. From there we find the minimum of $`E`$ which is the soliton mass, $`M`$. Assuming $`\chi _0\mu `$ this is $`M\mu ^3R^2`$. Since the total conserved quantum number is $`NsR^3`$, the total soliton mass is $`MN^{2/3}`$. The exponent of $`N`$ is less than unity which implies, that for large enough $`N`$, the total soliton mass is always less than the sum of free particles masses. Therefore this soliton configuration is stable.

Gravitational effects become important when $`R`$ approaches the Schwarzschild radius. By setting $`R_sG_4M`$ we estimate the critical mass, $`M_c`$, for the formation of a black hole:

$$M_c\left(\frac{M_{Pl}}{\mu }\right)^4\mu $$
(17)
For $`\mu 300`$GeV we get $`M_c10^{12}M_{Sun}`$ and $`R_s1`$ light month. This kind of solution is described in .

Consider now class II soliton model in the context of large extra dimensions. If the linear dimension, $`R`$, of a composite soliton object is much larger than the compactification radius of the extra dimensions, then the order of magnitude estimate presented above will not change much. However, on distances smaller than the compactification radius, the gravitational interaction and the functional dependence of the surface energy is quite different than in $`(3+1)`$-dimensional case. This will lead to the existence of the equilibrium configuration of the soliton states whose liner dimensions are microscopic.

Consider again a shell of radius $`R`$ which is smaller than the compactification radius of the extra dimensions. The kinetic energy inside the shell remains $`E_kN/R`$. But, the surface energy is $`E_ssR^{2+d}`$, where the surface tension $`s\mu \chi _0^{2+d}`$. By minimizing the total energy $`E=E_k+E_s`$ with respect to $`R`$ and assuming $`\chi _0\mu `$ we get:

$$M\mu ^{3+d}R^{2+d}.$$
(18)
Note that $`MN^{\frac{2+d}{3+d}}`$, and as in the case of $`(3+1)`$-dimensions, the exponent of $`N`$ is less than unity which implies the stability of the composite soliton state for large enough $`N`$.

Gravitational effects become important when $`R`$ approaches the $`(4+d)`$-dimensional Schwarzschild radius (12). Thus, we find the critical mass for the formation of a black hole:

$$M_c\left(\frac{M_F}{\mu }\right)^{(d+2)^2}\mu .$$
(19)
Note again that, if we take $`d=0`$ and replace $`M_F`$ with $`M_{Pl}`$, we recover $`(3+1)`$-dimensional results. Because of the form of the exponent in (19), the critical mass is very sensitive to the choice of the parameters in the model. Unlike in $`(3+1)`$-dimensional models of class II solitons, there is an equilibrium state of a composite soliton state which has microscopic dimensions. Like in case of class I soliton states, it is possible to set parameters close to the generic ones — for example $`\mu M_F1`$TeV — which gives $`M_c1`$TeV and $`R_s1`$TeV-1. However, note again that if the classical analysis of these states is to be believed one has to choose parameters which yield the soliton mass much higher than the fundamental quantum gravity scale $`M_F`$. Also, a large number $`N`$ of individual particles is needed in order to provide the classical stability of a composite soliton state which will result in the large total soliton mass.

Although a careful numerical analysis is required in order to estimate the total cross section production of these states in near future accelerator experiments (like the Large Hadron Collider with the center of mass energy of $`14`$TeV or the Very Large Hadron Collider with the center of mass energy of order $`100`$TeV ) it is conceivable that for some values of the parameters in the model the cross section would not be negligible.

In conclusion, we considered non-linear scalar filed theory models in the presence of gravity in the context of brane world models with large compact extra dimensions. We showed that new exotic non-topological soliton-like particles can exist. One class of these states arises in simple models which contain soliton solutions only in the presence of gravity. The other class arises in a different type of models which contain soliton solutions even in the absence of gravity. In this class, it is possible to have equilibrium states of these solitons which have both macro- and microscopic dimensions. The latter was not possible in $`(3+1)`$-dimensional models.

Depending on the abundance of these new states of matter, they can be important dark matter candidates not only because they are cold and massive but also because of the fact that great part of the dark matter must be of non-baryonic origin. It will be interesting to study the role of these states in cosmology and their contribution to entropy and energy density in our universe.

Acknowledgments: The author is grateful to Valeri Frolov, Glenn Starkman, Tanmay Vachaspati and Bruce Campbell for stimulating discussions. This work was partly supported by the Natural Sciences and Engineering Research Council of Canada. The author is also grateful to the Killam Trust for its financial support.


text:
# Chiral 𝑁𝑁 model and 𝐴_𝑦 puzzle

## Abstract

We analyze the results by chiral $`NN`$ models for the two-nucleon system and calculate the predictions for the nucleon vector analyzing power of elastic nucleon-deuteron ($`Nd`$) scattering, $`A_y`$, by these models. Our conclusion is that a quantitative chiral two-nucleon potential does not resolve the $`Nd`$ $`A_y`$ puzzle (when only two-body forces are included).

The term $`A_y`$ puzzle refers to the inability to explain the nucleon vector analyzing power $`A_y`$ in elastic nucleon-deuteron ($`Nd`$) scattering below 30 MeV laboratory energy for the incident nucleon by means of three-body calculations in which only two-nucleon forces are applied. The problem showed up as soon as it was possible to conduct three-body continuum calculations with realistic $`NN`$ potentials. The first such calculation was performed by Stolk and Tjon  in 1978 using the Reid soft-core potential , and the first calculations with (a separable representation of) the Paris potential  were conducted by the Graz-Osaka group in 1987 ; in both cases, the $`A_y`$ predictions showed the characteristic problem. Finally, the ‘puzzle’ became proverbial when rigorous three-nucleon continuum Faddeev calculations using realistic forces were started on a large scale . Over the years, many measurements and calculations of $`Nd`$ $`A_y`$ were performed (including the $`pd`$ reaction that involves the Coulomb force ) which all confirmed that the problem was real (see Ref.  for a review): For energies below 20 MeV, the $`A_y`$ is predicted about 30% too small in the angular region around 120 deg center-of-mass angle where the maximum occurs.

There have been many attempts to solve the problem. Already in the very early stages of three-body continuum calculations, when only schematic $`NN`$ potentials were applied, it was noticed that the $`Nd`$ $`A_y`$ predictions depend very sensitively on the strength of the input $`NN`$ potential in the triplet $`P`$ waves —a sensitivity that was confirmed in later calculations using realistic forces . Based upon this experience, Witała and Glöckle  showed in 1991 that small changes in those $`{}_{}{}^{3}P`$ wave potentials could remove the discrepancy. This finding gave rise to systematic investigations of the question whether the small variations of the low-energy phase shifts of, particularly, those triplet $`P`$ waves necessary to explain the $`Nd`$ $`A_y`$ are consistent with the $`NN`$ data base. While Tornow and coworkers  suggest that the low-energy $`NN`$ data may leave some lattitude in the $`NN`$ $`{}_{}{}^{3}P`$ waves that could improve the predictions for $`Nd`$ $`A_y`$, Hüber and Friar  find that it is not possible with reasonable changes in the $`NN`$ potential to increase the $`Nd`$ $`A_y`$ and at the same time to keep the two-body observables unchanged.

Another important observation has been that conventional three-nucleon forces (when added to a realistic two-nucleon potential) change the predictions for $`Nd`$ $`A_y`$ only slightly and do not improve them . Therefore, the general perception in the community has shifted towards the believe that the $`A_y`$ puzzle is the ‘smoking gun’ for new types of three-nucleon forces  or new physics .

However, very recently, there has been an apparent indication that the above conclusion/believe may be premature. It was reported  that with a two-nucleon force of a new type, namely, one that is based upon chiral effective field theory, the $`A_y`$ puzzle is resolved.

In recent years, effective field theory methods have become increasingly popular in nuclear physics. The reason for this development is the need to link conventional nuclear physics methods one way or the other to the underlying theory of strong interactions, QCD. After quark cluster models had only a limited success, it was recognized that the symmetries of QCD are more important than the high-energy degrees of freedom of QCD (quarks and gluons). The effective field theory concept distinguishes between different energy scales and assignes appropriate degrees of freedom for each scale while observing the over-all symmetries. For traditional nuclear physics with energies below 1 GeV, the right degrees of freedom are nucleons and pions interacting via a force that is controlled by (broken) chiral symmetry.

The derivation of the nuclear force from chiral effective field theory was initiated by Weinberg  and pioneered by Ordóñez  and van Kolck . Subsequently, many researchers became interested in the subject . As a result, efficient methods for deriving the nuclear force from chiral Lagrangians emerged and the quantitative nature of the chiral $`NN`$ potential improved. This trend shows up, in particular, in the excellent work by Epelbaum et al.  where the chiral $`NN`$ force was constructed using a unitary transformation and applying systematic power counting in next-to-leading order (NLO) and NNLO; and it is this potential in NLO that was applied in Ref. , resulting—seemingly—in a resolution of the long-standing $`A_y`$ puzzle. It is the purpose of this note to critically investigate the predictions by the chiral $`NN`$ model and the implications for the $`Nd`$ $`A_y`$ puzzle.

We start our investigation by taking a close look at important phase shifts of two-nucleon scattering. In Table I, we list $`S`$-wave phase shifts and, in Table II, we show $`{}_{}{}^{3}P`$-wave phase shifts for energies between 1 and 200 MeV. Since charge-dependence of the $`NN`$ interaction is not a crucial factor in the $`A_y`$ puzzle , and since the present chiral $`NN`$ potentials are all charge-independent and adjusted to the neutron-proton ($`np`$) data, we consider $`np`$ phase shifts and $`np`$ data.

It is of interest to compare the phase shifts produced by the chiral NLO (next-to-leading order) model by Epelbaum et al.  (which is the chiral potential applied in Ref. ) with the empirical ones from the Nijmegen multi-energy $`np`$ analysis  (PWA93) and the predictions by one representative of the family of the high-precision potentials constructed in the 1990’s (CD-Bonn ). In $`S`$ waves (Table I), there is, generally, good agreement up to 50 MeV. Above 50 MeV, differences between NLO and PWA93 show up and increase with energy. However, since the $`S`$-waves are not very important for $`Nd`$ $`A_y`$, this may not have much impact on the predictions.

We turn now to the triplet-$`P`$ waves (Table II) which are crucial for $`Nd`$ $`A_y`$, and focus, first, on the energy range below 30 MeV. The NLO $`{}_{}{}^{3}P_{0}^{}`$ phase shifts are about 2% larger than PWA93, which is not significant. For $`{}_{}{}^{3}P_{1}^{}`$, the differences are more drastic: The NLO value at 10 MeV is about 5% smaller and the one at 25 MeV is 3% smaller than PWA93. Finally, the NLO prediction for $`{}_{}{}^{3}P_{2}^{}`$ at 25 MeV is enhanced by 13%.

To understand how variations of the $`{}_{}{}^{3}P`$ phase shifts may effect observables, we consider the spin-orbit phase shift combination,

$$\mathrm{\Delta }_{\mathrm{LS}}=\frac{1}{12}\left(2\delta _{{}_{}{}^{3}P_{0}^{}}3\delta _{{}_{}{}^{3}P_{1}^{}}+5\delta _{{}_{}{}^{3}P_{2}^{}}\right),$$
(1)
which is a measure for the strength of the spin-orbit force. Results are shown in Table III. At 10 MeV, one obtains $`\mathrm{\Delta }_{\mathrm{LS}}^{\mathrm{PWA93}}=0.203`$ deg and $`\mathrm{\Delta }_{\mathrm{LS}}^{\mathrm{NLO}}=0.212`$ deg for PWA93 and NLO, respectively, implying that the NLO value is larger by 4.4% as compared to PWA93. At 25 MeV, the corresponding figures are $`\mathrm{\Delta }_{\mathrm{LS}}^{\mathrm{PWA93}}=0.93`$ deg and $`\mathrm{\Delta }_{\mathrm{LS}}^{\mathrm{NLO}}=1.09`$ deg implying that NLO is 17% larger. These numbers show that the spin-orbit force of NLO is enhanced as compared to the Nijmegen analysis and similar enhancements are obtained when comparing to any of the high-precision potentials, like CD-Bonn (cf. Table III). It is well-known that the two-nucleon spin-orbit force is magnified in the $`Nd`$ system. Therefore, a moderate enhancement of the $`NN`$ $`LS`$-potential leads to a substantial enlargement of the $`Nd`$ $`A_y`$. Thus, the enhanced NLO spin-orbit force as reflected in the larger values for $`\mathrm{\Delta }_{\mathrm{LS}}^{\mathrm{NLO}}`$ could very well be the explanation of the large $`Nd`$ $`A_y`$ predictions by NLO reported in Ref. .

Next, we consider the $`{}_{}{}^{3}P`$ phase shifts above 30 MeV. Here, the differences between PWA93 and CD-Bonn, on the one hand, and NLO, on the other, are in general larger and increase with energy (cf. Table II). This trend is most dramatic in $`{}_{}{}^{3}P_{2}^{}`$ where the discrepancies quickly grow into the hundreds of standard deviations. The sign of these differences are such as to drastically enhance the spin-orbit force, see Table III for the value at 50 MeV. Note that the $`NN`$ $`t`$-matrix, on- and off-shell, is input to the three-body continuum calculations. Thus, the description of the two-nucleon data at energies above 30 MeV has impact on low-energy three-body predictions. Therefore, the unrealistically strong NLO spin-orbit force above 30 MeV may—by means of an off-shell effect—further enhance the $`Nd`$ $`A_y`$ predictions. However, one may not have much confidence in this type of off-shell effect.

The above observations trigger the question if chiral models can also make more accurate predictions for $`NN`$ phase shifts; and if so, what are the implications for the $`Nd`$ $`A_y`$ problem, in such a case. The natural way is to include higher order terms in power counting which should improve not only the quality of the $`NN`$ phase shift reproduction but also extend the energy range in which it works. For that purpose, we pick up the chiral $`NN`$ potential of Ref.  (subsequently denoted by ‘Idaho’) that was recently developed. In the chiral Idaho model , contact terms up to order four are included which introduces more parameters allowing for a better fit of the lower partial waves in a much wider energy range. In Table I and II, it is clearly seen that the chiral Idaho $`NN`$ potential reproduces the empirical $`np`$ phase shifts of PWA93 up to 200 MeV, accurately.

We consider now the observable $`Nd`$ $`A_y`$ which is the focus of this paper. We have calculated the predictions for $`Nd`$ $`A_y`$ at energies 3, 10, and 65 MeV for the incident nucleon. The results are shown in Fig. 1. In this figure, the shaded band represents the prediction by the family of high-precision potentials (using always the $`np`$ version of those models), namely, CD-Bonn , Argonne $`V_{18}`$ , and the Nijmegen potentials Nijm-I, Nijm-II, and Nijm93 . The dashed line is the predicition by the Idaho chiral $`NN`$ potential  and it is clearly seen that this prediction follows accurately the narrow band made up from the variations among those high-precision potentials. In conclusion, at 3 and 10 MeV, we are still having an $`A_y`$ problem if the chiral $`NN`$ potential is a quantitative one.

The evidence presented may be perceived as a convincing proof that a quantitative chiral potential does not resolve the $`Nd`$ $`A_y`$ puzzle. However, there remains one objection that can be raised. In the literature, notably in Ref. , one can find the suggestion that the $`{}_{}{}^{3}P`$ waves at low energy are not as well determined as claimed in PWA93 . If true, then moderate variations of the $`{}_{}{}^{3}P`$ phase shifts at low energy could be consistent with the low-energy $`NN`$ data. This variations could be such as to enhance the low-energy spin-orbit force and, thus, lead to an improved prediction for $`Nd`$ $`A_y`$. For the purpose of seriously checking out this possibility, we have constructed a variation of the Idaho chiral potential with modified $`{}_{}{}^{3}P`$ phase shifts at low energy. The column ‘Modified’ of Table II shows the $`{}_{}{}^{3}P`$ phase shifts and the corresponding column of Table III reveals that for this fit the spin-orbit force is enhanced, similarly to NLO. However, in contrast to NLO, the ‘modified’ model is much more realistic, since the phase shifts do not diverge to unrealistic values at higher energies.

We have calculated the $`Nd`$ $`A_y`$ as predicted by the ‘modified’ chiral model and find, indeed, a considerable improvement (see dotted curve in Fig. 1). Is this the resolution of the $`A_y`$ puzzle by a chiral $`NN`$ potential? To answer this question one needs to know if the modified chiral model is a realistic and quantitative $`NN`$ potential. A precise and reliable answer cannot be given by just looking at phase shifts. As stressed repeatedly by the Nijmegen group in the past, only a direct confrontation with the $`NN`$ data can reveal if a $`NN`$ potential is quantitative or not. For that reason we have calculated the $`\chi ^2`$/datum for the fit of the world $`np`$ data as represented by the 1999 database , see Table IV. For a proper interpretation of the results of Table IV, it is necessary to establish a standard concerning what $`\chi ^2`$/datum represents a ‘quantitative’ reproduction of the data. This issue was debated a lot in the 1990’s, and the consensus that emerged was that only values below 1.1 are acceptable. Deliberatly, we losen this standard and consider a fit with $`\chi ^2/\mathrm{datum}1.2`$ as ‘quantitative’, while we will perceive higher values as not acceptable.

Applying this standard, the modified model produces unacceptable values for $`\chi ^2`$/datum for all energy intervals below 75 MeV (cf. Table IV). Thus, the modified chiral model is not a quantitative one and, consequently, it is not the resolution of the $`Nd`$ $`A_y`$ puzzle. Since it is well known that the off-shell character of the $`NN`$ potential plays essentially no role in three-nucleon scattering, one can further draw the more general conclusion: No model that reproduces the $`NN`$ data correctly can solve that $`Nd`$ $`A_y`$ puzzle.

Table IV shows also the $`\chi ^2`$/datum of the other models discussed in this paper. It is seen that PWA93, Idaho, and CD-Bonn reproduce the $`np`$ data below 210 MeV with the perfect $`\chi ^2/\mathrm{datum}=0.97`$ and 0.98. The chiral NLO potential by Epelbaum et al.  produces $`\chi ^2/\mathrm{datum}=37`$ which is grossly unacceptable. In fact, only for the interval 0-8 MeV is NLO acceptable. This range of validity is so tiny that no serious implications can be drawn from any prediction by this potential.

Finally, we like to take this opportunity to also present an overview of other interesting $`Nd`$ observables, which are shown in Figs. 2-8. Concerning the deuteron vector analyzing power, $`iT_{11}`$, at 3 and 10 MeV (Fig. 2), it should be noted that a seemingly drastic reduction of the discrepancy between theory and $`pd`$ data seen at 3 MeV has its origin in large effects of the long range Coulomb force acting between two protons, which is not taken into account in our calculations . Taking this Coulomb force into account, $`iT_{11}`$ is also underpredicted in the peak region —an equally well-known problem, which is why it would be appropriate to speak more generally of the vector analyzing power puzzle in elastic $`Nd`$ scattering. In almost all cases, the Idaho chiral $`NN`$ potential follows the trend of the high-precision potentials. The only exceptions are $`T_{20}`$ and $`T_{21}`$ at 65 MeV where the chiral potential predictions describe slightly better the data in the minimum region as compared to conventional potentials. But, apart from this, the quantitative chiral $`NN`$ model containing contributions of higher orders in power counting does not produce any new signatures.

In summary, our main conclusions are:

* A quantitative chiral two-nucleon potential does not resolve the $`Nd`$ $`A_y`$ puzzle on the two-body force level.
* Low-energy $`{}_{}{}^{3}P_{J}^{}`$ $`NN`$ phase shifts that “solve” the $`Nd`$ $`A_y`$ puzzle are inconsistent with the low-energy $`NN`$ data.

And, finally, as a consequence of the above two points, one may expect that no quantitative two-nucleon force—no matter what the basis is, pure phenomenology, meson theory, chiral EFT, or anything—will ever solve the $`Nd`$ $`A_y`$ puzzle.

An accurate $`NN`$ model requires to take chiral perturbation theory ($`\chi `$PT) to order four. At that order, also many three-body force (3NF) terms occur. According to the basic rules of $`\chi `$PT, all two- and many-body terms must be included for a complete calculation. Conventional 3NFs were shown to be ineffective for the $`A_y`$ problem . The advantage of $`\chi `$PT is that it provides a systematic scheme to generate all terms at a given order. As it turns out, there are several such chiral 3NF terms that were never considered in few-nucleon physics before . It is natural to expect the resolution of the $`Nd`$ $`A_y`$ puzzle from such new chiral 3NFs which, therefore, should be at the focus of future work in the field.

At next-to-leading order (NLO) in $`\chi `$PT, there are no 3NF contributions. So, a calculation with a NLO two-nucleon potential and no 3NF seems to have formal validity. However, since at NLO the $`NN`$ data can only be reproduced for $`T_{\mathrm{lab}}8`$ MeV, such a calculation is doomed to be inconclusive, from the outset, and higher order terms must be taken into account for any meaningful calculation.

###### Acknowledgements.

This work was supported in part by the U.S. National Science Foundation under Grant No. PHY-0099444, by the Ramón Areces Foundation (Spain), and by the U.S. Department of Energy, Office of High Energy and Nuclear Physics, under Grant No. DE-FG02-97ER41033. The numerical calculations have been performed on the CRAY T90 of the John von Neumann Institute for Computing in Jülich, Germany.


text:
# Evidence for “fragile” glass-forming behavior in the relaxation of Coulomb frustrated three-dimensional systems.

## Abstract

We show by means of a Monte Carlo simulation study that three-dimensional models with long-range frustration display the generic phenomena seen in fragile glassforming liquids. Due to their properties (absence of quenched disorder, physical motivation in terms of structural frustration, and tunable fragility), these systems appear as promising minimal theoretical models for describing the glass transition of supercooled liquids.

The viscous slowing down of supercooled liquids that leads to glass formation when the temperature is sufficiently lowered, remains a challenging problem of liquid-state and condensed-matter physics. The phenomenon seems dramatic enough (a change by $`15`$ orders of magnitude of the viscosity $`\eta `$ and the primary ($`\alpha `$) relaxation time $`\tau _\alpha `$ for a mere decrease of the temperature $`T`$ by a factor of $`2`$) and general enough to call for a universal explanation. Yet, the search for an underlying universality is thwarted by the absence of observed criticality, e.g., by the absence of a diverging, or rapidly growing, static susceptibility or correlation lengthEAN96 .

If one accepts the premise that the distinctive features of most glassforming liquids, namely the rapid super-Arrhenius increase of $`\eta `$ and $`\tau _\alpha `$ with decreasing $`T`$ and the non exponential character of the relaxation functions, result from collective or cooperative behavior, progress could be made by identifying the key physical ingredients and formulating a minimal theoretical model  based on these ingredients, in much the same way as spin glasses have been described in terms of quenched disorder and frustration and modelled by the Edwards-Anderson Hamiltonian. One such attempt has been recently madeKKZNT95 , based on the concept of structural frustrationNelsonS89 . Such frustration arises because a locally preferred structure, i.e., an arrangement of neighboring molecules that minimizes some local free energy, cannot tile the whole space, because of the global constraintNelsonS89 . The canonical example is a $`3`$-dimensional system of spherical particles whose locally preferred structure, an icosahedral cluster of $`13`$ particles, cannot form a pure crystal because of the $`5`$-fold rotational symmetry of the icosahedronNelsonS89 . In the frustration-limited domain theory (FLDT) of supercooled liquidsKKZNT95 , structural frustration is implemented as the competition between effective interactions acting on very different length scales: a short-range ordering term that favors the extension of the locally preferred structure and a weak, but long-range ($`1/r`$) frustrating term that generates a strain free energy growing super-extensively with the linear size of the ordered region. These ingredients are most simply incorporated in three-dimensional Coulomb ($`1/r`$) frustrated lattice models with the following Hamiltonian,

$$H=J\underset{<ij>}{}𝐒_i.𝐒_j+\frac{Q}{2}\underset{ij}{}\frac{𝐒_i.𝐒_j}{|𝐫_{ij}|},$$
(1)
where $`J`$ and $`Q`$ are both positive and denote the strength of the ordering and frustrating interactions, respectively, $`<ij>`$ is a sum over distinct pairs of nearest neighbors and $`|𝐫_{ij}|`$ is the distance between sites $`i`$ and $`j`$ (the lattice spacing is taken as unit length); the tunable ratio $`Q/J`$ characterizes the frustration strength. In such models, glassiness has to be self-induced and does not result from quenched randomness or from specific kinetic constraints.

In this letter, we present an extensive Monte Carlo simulation study of Coulomb frustrated models on a cubic lattice. We focus on the relaxation to equilibrium and we do not consider at this stage out-of-equilibrium dynamics and aging phenomena in a non-equilibrated, glassy state. We show that: (i) the relaxation displays the main features characteristic of the slowing down of “fragile” glassforming liquids (non exponential decay of the relaxation function, super-Arrhenius activated $`T`$dependence of $`\tau _\alpha `$); (ii) the “fragility” of the system, e.g., its degree of departure from Arrhenius behavior, increases as the frustration decreases, thereby providing a description that spans the whole range from strong to very fragile glassformer via tuning of the relative amplitude of the frustrating interaction; (iii) the slowing down occurs without the rapid growth of a correlation length and does not follow the standard critical slowing down pattern; (iv) the results, including the existence of a crossover temperature in the vicinity of the critical point of the unfrustrated system, confirm the scenario proposed by the FLDT.

In using Eq. (1), we have chosen spin variables $`𝐒_i`$ that can take several orientations; most of the results have been obtained for the $`5`$-state clock model ($`𝐒_i=(\mathrm{cos}(2\pi \theta _i/q),\mathrm{sin}(2\pi \theta _i/q))`$, with $`\theta _i`$ the orientation of the planar spin and $`q=5`$) and the Ising model ($`q=2`$), but we have also studied the $`q=11`$ state clock model.

In the absence of frustration ($`Q=0`$), the models have a critical point at a temperature $`T_c^0`$, below which they are ferromagnetic. The long-range, coulombic interaction requires that the total magnetization of the system be zero in order to ensure a proper thermodynamic limit. Therefore, long-range ferromagnetic order is prohibited at all $`T^{}s`$ for any nonzero value of the frustration parameter $`Q/J`$; in the FLDT, this mimics the effect of structural frustration that forbids any long-range order based on the periodic repetition of the locally preferred structure. The models with discrete orientations studied here can still form ordered phases characterized by modulated patterns (e.g., lamellar phases)Dzu93 , but as shown recently by Monte Carlo simulation of the Ising modelGTV01b , the transition between paramagnetic and modulated phases is first-order. Following BrazovskiiB75 , this result can be interpreted on the basis of the self-consistent Hartree approximation which predicts the occurence of a fluctuation-induced first-order transition, a transition with no nearby low-$`T`$ spinodal. All the simulations discussed below are made in the disordered (paramagnetic) phase. Since the transition to the modulated phases at $`T_{DO}(Q)`$ is first-order, one could in principle supercool the paramagnetic phase to lower temperatures; however, we have found that the lattice sizes achievable in practice are too small to allow for a proper supercooling below $`T_{DO}`$.

To study the relaxation of the frustrated systems, we have considered the dynamics associated with the Monte Carlo algorithm (simple or slightly modified Metropolis ruleGrous2001 ), time being the number of sweeps per spin. In all cases, the total magnetization is constrained to be zero, as it should be in the thermodynamic limit. The dynamical quantity that we have monitored is the spin-spin self correlation function, $`C(t)=(1/N)_i<𝐒_i(t^{})𝐒_i(t^{}+t)>`$, where the bracket denotes the thermal average and $`N`$ is the total number of lattice sites. We stress that $`C(t)`$ is an equilibrium correlation function that is computed once the system has reached equilibrium. In practice, the thermal average has been performed as an average over $`20`$ different initial times $`t^{}`$, all chosen longer than the relaxation time. Most simulations have been done on cubic lattices of size $`16^3`$ to $`22^3`$ with periodic boundary conditions, and the Coulomb interaction has been handled via Ewald sums. In one special case (see below), we have performed a finite-size study with lattices of linear size $`L=5,10,15`$ and $`20`$. We have covered two orders of magnitude in the frustration parameter $`Q/J`$, always in the weak-frustration regime $`(Q/J<<1)`$. In what follows, we take $`J=1`$.

The evolution of $`C(t)`$ with $`T`$ is illustrated in Fig. 1 for the $`5`$-state clock ($`Q=0.1`$ and $`Q=0.00625`$), the $`11`$state clock ($`Q=0.05`$), and the Ising ($`Q=0.001`$) models. At the highest frustration, illustrated by Fig. 1a, the decay of $`C(t)`$ appears to proceed in a single step at all $`T`$’s. For lower frustrations, a $`2`$step decay develops as $`T`$ is lowered (see Fig. 1b and d). At high $`T`$’s, typically above the critical point of the coresponding unfrustrated model ($`T_c^02.1`$ and $`T_c^04.51`$ for the $`5`$state clock model and the Ising model, respectively), the whole time dependence of $`C(t)`$ is well fitted by a simple exponential, $`\mathrm{exp}(t/\tau _0(T))`$; this is shown in Fig. 2, but is true for all models and all frustrations. At low $`T`$’s, below $`T_c^0`$, it is impossible to describe the entire decay of $`C(t)`$ by a single exponential; as illustrated in Fig. 2, the emerging second step of the relaxation can be described by a stretched exponential, $`\mathrm{exp}((t/\tau _{KWW}(T))^{\beta (T)})`$, while the first step can still be fitted by a simple exponential, $`\mathrm{exp}(t/\tau _0(T))`$. This non-exponential behavior and emergence of a two-step decay, both of which become more marked as $`T`$ is decreased, are typical of many actual glassforming systems, especially the fragile supercooled liquidsEAN96 . One may however notice that the timescale separation between the two relaxation steps that can be achieved in the simulations is not sufficient to observe the development of a true plateau at intermediate times.

To study in more detail the characteristics of the slowing down of the relaxations as one lowers $`T`$, a phenomenon already visible in Fig. 1, we have monitored as a function of $`T`$ the relaxation times characteristic of the short- and long-time behaviors, $`\tau _0(T)`$ and $`\tau _\alpha (T)`$ respectively. The former is obtained from the fit to a simple exponential (see above), but an equivalent behavior is obtained from the condition $`C(\tau _0)=0.9`$. The latter is defined, as often for spin systems, by the condition $`C(\tau _\alpha )=0.1`$; in this case indeed, there was too much uncertainty in using the three-parameter stretched exponential fit that also requires the introduction of a lower cut-off time. The results are illustrated in Fig. 3 for the $`5`$-state clock model. We have plotted the logarithm of $`\tau _0`$ and $`\tau _\alpha `$ versus $`1/T`$. One can see that over the range of $`T`$ studied, the dependence of $`\tau _0(T)`$ is essentially Arrhenius-like, i.e., linear on the diagram, $`\tau _0(T)\tau _{0,\mathrm{}}\mathrm{exp}(E_{\mathrm{}}/T)`$, whereas the dependence of $`\tau _\alpha (T)`$ shows a marked deviation from Arrhenius behavior below some crossover temperature in the vicinity of $`T_c^0`$. Note that at high $`T`$’s (above $`T_c^0`$), $`C(t)`$ is essentially exponential and the two times $`\tau _0(T)`$ and $`\tau _\alpha (T)`$ differ by a trivial constant, so that in this $`T`$ range the Arrhenius dependences of $`\tau _0`$ and $`\tau _\alpha `$ are the same. At lower $`T`$’s, $`\tau _0`$ is characteristic of a “secondary” or “precursor” relaxation; one expects it to be associated with a local, weakly or non cooperative, mechanism of relaxation. On the other hand, the super-Arrhenius rise of $`\tau _\alpha `$ versus $`1/T`$ is indicative of cooperative behavior.

The crossover from Arrhenius to super-Arrhenius behavior of $`\tau _\alpha `$ is shown in Fig. 4 for the $`5`$state clock model for $`5`$ different frustrations from $`Q=0.1`$ down to $`Q=0.00625`$. A similar trend is observed for all other models studied (for the Ising model, see Ref.GTV01a ; Grous2001 )). Several points are worth noting: first, such a super-Arrhenius behavior is typical of the viscous slowing down of fragile glassforming liquids, the more fragile a liquid the more pronounced the super-Arrhenius characterEAN96 ; secondly, the crossover occurs in the vicinity of the critical point of the unfrustrated system; thirdly, the departure from Arrhenius behavior, i.e., the “fragility”, becomes more marked as frustration decreases: see Fig. 5. To our knowledge, such a “tunable” fragility that provides a description over the whole spectrum of glassforming behavior, from strong to very fragile, is shown by no other microscopic model.

The competition of interactions acting on very different length scales can thus generate fragile glassforming behavior similar to that observed in actual supercooled liquids. In the FLDTKKZNT95 , this behavior is attributed to the proximity of the critical point $`T_c^0`$ of the unfrustrated system, a critical point that is avoided in the presence of frustration. Avoided critical behaviorKKZNT95 ; CEKNT96  results in the appearance at $`T`$’s below $`T_c^0`$ of domains whose size is limited by the frustration. These domains are responsible for the non-exponential character of the relaxation and the super-Arrhenius, activated $`T`$dependence of $`\tau _\alpha `$. Scaling arguments predict for instance that the effective activation free energy for the $`\alpha `$ relaxation behaves as KKZNT95

$`E_\alpha (T)`$ $`=`$ $`E_{\mathrm{}},\text{ }T>T^{},`$ (2)
$`=`$ $`E_{\mathrm{}}+BT^{}\left(1{\displaystyle \frac{T}{T^{}}}\right)^{8/3},\text{ }T<T^{},`$ (3)
where $`E_\alpha (T)=T\mathrm{ln}\left(\frac{\tau _\alpha (T)}{\tau _{\mathrm{}}}\right)`$, with $`\tau _{\mathrm{}}`$ a characteristic high-T time, and $`T^{}`$ is close to $`T_c^0`$. The existence of a crossover in the vicinity of $`T_c^0`$ is indeed supported by the Monte Carlo results. More quantitatively, we have compared the above prediction with the simulation data. In doing so, one must remember that the range of $`\tau _\alpha `$’s spanned by our data is limited, as in most simulations of glassforming systems, and is far less than the $`15`$ decades observed experimentally in supercooled liquids. The high-$`T`$ parameters $`\tau _{\mathrm{}}`$ and $`E_{\mathrm{}}`$ are independently obtained by an Arrhenius fit to $`\tau _0(T)`$ (see above, where $`\tau _{\mathrm{}}/\tau _{0,\mathrm{}}=\mathrm{ln}(0.1)`$ due to the definition of $`\tau _\alpha `$ and $`\tau _0`$), and one is left with two adjustable parameters, $`B`$ and $`T^{}`$, for each frustration (and each model). The resulting curves are shown in Fig. 4 for the $`5`$state clock model. For all models the agreement between theory and simulation is good. As predicted by the theory, $`T^{}`$ is found close to $`T_c^0`$ (e.g, $`T^{}/T_c^01.141.20`$ for the $`5`$-state clock and $`1.041.12`$ for the Ising model) and the super-Arrhenius parameter $`B`$ varies roughly as $`Q^1`$. (As already stressed, small frustration $`Q`$ implies large super-Arrhenius character and large fragility.)

We have thus found that our simulation data are compatible with an activated-type expression, $`\tau _\alpha (T)\mathrm{exp}(E_\alpha (T)/T)`$ with $`E_\alpha (T)`$ given by the FLDT expression. (As discussed in SW00 ; GTV01a , the data are also compatible with $`E_\alpha (T)`$ described by the entropic droplet/random-first-order transition approach (EDA)KTW89 ; SW00 , which offers a unique opportunity for comparing these two alternative theories in more detail.) On the other hand, we have checked that they are poorly described by a power-law dependence, $`\tau _\alpha (T)|TT_c|^\gamma `$, except of course when $`Q=0`$Grous2001 . In both the FLDT and the EDA, the slowing down of the relaxation, as well as the non exponential character is attributed to genuinely non-perturbative effects, such as heterogeneities (domains or droplets) and cooperative, activated processes. This is quite different than standard critical slowing down, as could occur if driven by the approach of a spinodal. In particular, the length scale associated with the slowing down is predicted to grow slowly with decreasing $`T`$, at a much slower rate than the relaxation time itself. To check this, we have performed a finite-size study of the $`5`$-state clock model for $`Q=0.00625`$. As seen in Fig. 5a, $`\tau _\alpha (T)`$ is indeed unchanged when one decreases the linear size of the system from $`20`$ to $`10`$; on the other hand, it is strongly altered when passing to $`L=5`$. This seems to confirm that any putative length scale associated with the strong slowing down is rather modest (of the order of $`5`$ or so) over the whole $`T`$ range studied. For comparison, the same finite-size study for the unfrustrated ($`Q=0`$) model, showing the usual pattern associated with critical slowing down, is displayed in Fig. 5b.

In summary, due to their properties (absence of quenched disorder, physical motivation in terms of structural frustration, $`3`$dimensional, finite-range character, tunable fragility), whose combination is found in no other modelsCLMYB , the Coulomb frustrated systems appear as promising minimal theoretical models for describing the slowing down of relaxations in supercooled liquids.

We thank D. Kivelson and L. Berthier for many stimulating discussions.


text:
# Anisotropy of Electrons Accelerated by a High-Intensity Laser Pulse

## Abstract

We describe a realistic model for a focused high-intensity laser pulse in three dimensions. Relativistic dynamics of an electron submitted to such pulse is described by equations of motion with ponderomotive potential depending on a single free parameter in the problem, which we refer to as the ”asymmetry parameter”. It is shown that the asymmetry parameter can be chosen to provide quantitative agreement of the developed theory with experimental results of Malka et al. \[Phys. Rev. Lett. 78, 3314 (1997)\] who detected angular asymmetry in the spatial pattern of electrons accelerated in vacuum by a high-intensity laser pulse.

In their recent paper, Malka et al.  have reported experimental observation of electrons accelerated to relativistic energies by a high-intensity linearly polarized subpicosecond laser pulse in vacuum (see also comments on the paper  and the author’s reply in Ref. ). This effect, known as high-intensity ponderomotive scattering, was discussed in detail in Ref. . It occurs when the quiver amplitude imparted by the laser field to an electron becomes comparable to the focal spot radius of the laser beam. If the beam is Gaussian, the radial restoring force acting on the electron decays exponentially and the electron can be scattered out of the pulse.

The data of Malka et al.  show that the energies gained by the scattered electrons are in good quantitative agreement with calculations of electron trajectories in the polarization plane made with the first-order paraxial model for the laser field . Nevertheless, the first order paraxial model predicts isotropic electron scattering  that is not supported by experimental results. Indeed, accelerated electrons were detected by Malka et al. only in the $`(𝐄,𝐤)`$ plane, while no significant signal was detected after rotating the laser polarization direction by $`90^o`$ . In our opinion, this discrepancy between the theory and experiment is due to the following. In the first-order paraxial model focusing of a plane monochromatic wave only leads to an appearance of nonvanishing longitudinal components of electromagnetic fields in the focal region. However, the focusing is known to affect transverse components of the fields also (see, e.g. Ref. ). In particular , a plane monochromatic wave, polarized linearly along the $`x`$ axis and propagating along the $`z`$ axis, is converted by an aplantic system to a converging spherical wave with non-vanishing $`y`$ and $`x`$ components of electric and magnetic fields, respectively.

In this Letter, we show that experimentally observed  anisotropy of ponderomotive electron scattering can be explained in the framework of a realistic model for the laser field developed in our recent paper . Our model is based on an exact solution of Maxwell equations in three dimensions (3D), which can serve to describe a stationary, focused monochromatic laser beam with characteristic frequency $`\omega `$ and arbitrary intensity. Amplitudes of electric and magnetic fields in the model depend on radial coordinates as well as the coordinate along the direction of the beam propagation. These amplitudes are characterized by parameters $`R`$ and $`L=\omega R^2`$, which can be interpreted as the focal spot radius and the Rayleigh length of the laser beam, respectively. The model admits different field configurations, which are determined by two coordinate functions satisfying certain second-order partial differential equations. Some special choice of these functions describes the Gaussian beams, which are widely used in optics. The model can be generalized by introducing temporal amplitude envelope $`g(\phi /\omega \tau )`$, where $`\phi `$ is the relativistically invariant phase of the traveling wave, to describe a laser pulse with finite duration $`\tau `$. (It is assumed that the function $`g(\phi /\omega \tau )`$ is equal to unity at the point $`\phi =0`$ and decreases exponentially at the periphery of the pulse for $`|\phi |\omega \tau `$.) In this case the electric and magnetic fields of the model constitute an approximate solution of Maxwell equations with the second-order accuracy with respect to small parameters $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }^{}`$, defined as

$$\mathrm{\Delta }^{}=1/\omega \tau \mathrm{\Delta }=1/\omega R1.$$
(1)
For a pulse propagating along the $`z`$ direction, $`x`$ and $`y`$ components of the electric field oscillate with phase difference, which depends on $`z`$ and values of $`x`$ and $`y`$ coordinates of a point in the plain $`z=const`$. Moreover, the aforementioned phase difference depends also on $`\phi `$. Therefore, one cannot ascribe some definite type of polarization to a tightly focused laser pulse. Nevertheless, for a weakly focused pulse ($`\mathrm{\Delta }1`$), there always exists a region near the axis of the beam $`rR`$, where the field properties are very close to those of a plane wave field. This region we call ”the plane wave zone”. It is reasonable to ascribe polarization of the field in this region to the beam as a whole. Hereafter we refer to the field of the pulse as linearly polarized in this sense only. For a tightly focused beam the focal spot radius is of the order of the wavelength and the plane wave zone doesn’t exist. Therefore, only polarization of the parental beam incident on the focusing optical system can be ascribed to the focused beam in this case.

An arbitrary field, linearly polarized along the $`x`$ axis, may be represented  as a superposition of $`E`$\- and $`H`$-polarized waves (i.e. waves with the vectors $`\stackrel{}{E}`$ and $`\stackrel{}{H}`$ being perpendicular to the direction of the pulse propagation, compare with Ref. ). Relative contributions of $`E`$\- and $`H`$-polarized waves to the resulting field are characterized by the ”asymmetry parameter” $`\mu `$

$$\mu =\frac{E_{x0}^hE_{x0}^e}{E_{x0}^h+E_{x0}^e},\mathrm{}<\mu <\mathrm{},$$
(2)
where $`E_{x0}^{e,h}`$ are the $`x`$ components of the electric field for $`E`$\- and $`H`$-polarized waves at the focal point $`𝐫=0`$ for $`\phi =0`$. Note, that in contrast to the amplitude, the quantities $`E_{x0}^{e,h}`$ can take both positive and negative values.

In the lowest approximation in $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }^{}`$, the averaged equations of motion of electrons (ponderomotive equations) in the field of the linearly polarized laser pulse take the form

$$\begin{array}{cccccc}\hfill \frac{d\stackrel{}{q}_{}}{d\phi }& =& \mathrm{\Delta }\frac{m}{q_{}}\frac{U}{\stackrel{}{\rho }_{}},\hfill & \hfill \frac{d\stackrel{}{\rho }_{}}{d\phi }& =& \mathrm{\Delta }\frac{\stackrel{}{q}_{}}{q_{}},\hfill \\ & & & & & \\ \hfill \frac{dq_{}}{d\phi }& =& 0,\hfill & \hfill \frac{d\zeta }{d\phi }& =& \mathrm{\Delta }^2\frac{q_z}{q_{}}.\hfill \end{array}$$
(3)
Here $`q^\mu =p^\mu `$, where $`p^\mu `$ is the 4-momentum of the electron, $`q_{}=q_0q_z`$, $`\stackrel{}{\rho }_{}=\stackrel{}{r}/R`$, $`\zeta =z/R`$, brackets $``$ mean averaging over fast oscillations, and the ponderomotive potential $`U`$ is defined by the expression

$$U=\frac{m\eta _0^2}{2}g^2(\phi /\omega \tau )\left\{\left|F_1\right|^2+\mu ^2\left|F_2\right|^2+\mu \mathrm{cos}2\psi \left(F_1F_2^{}+F_1^{}F_2\right)\right\},$$
(4)
where $`\mathrm{tan}\psi =\varrho _x/\varrho _y`$, and $`\eta _0`$ is the value of the dimensionless field intensity parameter

$$\eta ^2=\frac{e^2\text{E}^2}{m^2\omega ^2},$$
(5)
at the focal point at the moment $`\phi =0`$ ($`\eta _0=a/\sqrt{2}`$, where $`a`$ is the parameter of Malka et al.). Functions $`F_i(\stackrel{}{\varrho }_{},\zeta ;\mathrm{\Delta })`$ are chosen in the form corresponding to the Gaussian beam

$$\begin{array}{c}F_1=(1+2i\zeta )^2\left\{1\frac{\varrho _{}^{}{}_{}{}^{2}}{1+2i\zeta }\right\}\mathrm{exp}\left\{\frac{\varrho _{}^{}{}_{}{}^{2}}{1+2i\zeta }\right\},\\ \\ F_2=\varrho _{}^{}{}_{}{}^{2}(1+2i\zeta )^3\mathrm{exp}\left\{\frac{\varrho _{}^{}{}_{}{}^{2}}{1+2i\zeta }\right\}.\end{array}$$
(6)
The equation (4) shows that the ponderomotive potential $`U`$ depends on the azimuthal angle $`\psi `$, and hence is generally speaking asymmetric. The potential $`U`$ is symmetric only for the case $`\mu =0`$. The shape of the ponderomotive potential in the plane $`z=0`$ for $`\phi =0`$ is shown in Fig.1 for the cases $`\mu =0`$ and $`\mu =1.55`$. Figure 1a represents the ponderomotive potential for the standard case of Gaussian beam commonly used in literature, while Figure 1b illustrates the dramatic difference between the cases $`\mu =0`$ and $`\mu 0`$. For $`\mu 0`$, the ponderomotive potential, possesses (besides the central peak) two extra maxima, which are located in the polarization plane. They arise as a result of the non-uniform intensity distribution in the plane $`z=0`$ for $`\mu 0`$. Locations of the additional maxima, as well as their amplitudes, are determined by the value of $`\mu `$. It is noteworthy, that the case $`\mu =0`$ is the only one when $`E_y`$\- and $`H_x`$-components of the electric and magnetic fields remain to be equal to zero outside the plane wave zone for the pulse polarized along the $`x`$ axis.

We use Eqs. (3) and (4) for our analysis of free electron acceleration by a co-propagating intense laser pulse in vacuum under conditions close to those used in the experiments of Malka et al. . The initial electron energy is taken to be $`\epsilon =10`$ keV ($`v_0=0.2c`$). The laser field parameters are: $`\lambda =1`$ $`\mu `$m, $`\eta _0=2.12`$ (corresponds to the parameter $`a=3`$ of Malka et al.), $`R=10`$ $`\mu `$m and $`\omega \tau =480`$. For the temporal envelope of the pulse $`g(\phi /\omega \tau )`$ a sine-squared shape is taken. The asymmetry parameter is not determined experimentally and remains a free parameter of the problem. Its value, $`\mu =1.55`$, has been chosen for better fitting of our computational results to the results of the experiment.

It is clear that the maximum energy will be gained by electrons that initially propagate exactly along the axis of the laser beam. However, the ponderomotive equations yield zero net energy transfer for such electrons, since they can ”feel” the spatial gradient of the Gaussian laser field, and hence can be scattered out of the pulse, only due to the quiver motion which is absent for the average trajectory described by the ponderomotive equations. Nevertheless, there exist a family of trajectories with nonzero initial distances from the beam axis, for which the gained energy is close to its maximum value. The energy and scattering angle of the electrons depend also on the position at which the particle is overtaken by the pulse . Indeed, the maximum energy is obtained for the electrons that experience the peak field of the laser and therefore meet the laser pulse at some distance before the focus. As a result, we obtain a 3D domain of injection positions of electrons obtaining the final kinetic energy $`W0.9`$ MeV at the scattering angle $`39.5^o`$.

Different cross-sections of this domain for the potential with $`\mu =1.55`$ (Fig. 1b) are shown in Fig. 2. The longitudinal size of the domain is of the order of the Rayleigh length for the laser beam $`L`$, whereas its transverse size is much less than the focal spot radius $`R`$. The cross sections in Fig. 2 display high degree of radial anisotropy. Their shape, of course, is essentially determined by the type of the ponderomotive potential or by the value of $`\mu `$. In particular, the domain of injection positions for the case (not presented here for the sake of compactness) of the potential shown in Fig. 1a with $`\mu =0`$ is purely radial, in agreement with .

To obtain the angular distribution of scattered electrons we have calculated their trajectories numerically, applying the fourth-order Runge-Kutta method to the ponderomotive equations (3). Initial positions of the electrons were taken from the domain shown in Fig. 2. Since the shape of the cross-sections $`z=const`$ varies very slowly at intervals $`\delta zR`$, we have considered equidistant planes $`z=nR`$ with $`n=27,26,\mathrm{},5`$. In each of these cross-sections, the electron injection positions were chosen randomly under condition that their density was constant and equal to $`3\times 10^{15}/R^2`$. Physically, such procedure corresponds to uniformity of the initial electron beam. The total number of electron trajectories considered in such a way was more than $`2\times 10^7`$. We were interested only in those trajectories, which crossed the plane $`z=11.66`$ cm at the points contained inside the ring with radii $`r_1=8.99`$ cm and $`r_2=9.89`$ cm. The latter conditions were determined by the position and angular size of the detector in the experiment .

The results of the calculations are presented in Fig. 3. We plot the normalized number of scattered electrons $`n`$ with final energies $`W0.9`$ MeV as a function of the azimuthal scattering angle $`\alpha `$. The values of $`n`$ at any given $`\alpha `$ were obtained by averaging the number of scattered electrons over the range of $`\alpha `$ equal to the angular size of the detector used in the experiments . It is easily seen, that the angular distribution of scattered electrons essentially depends on the parameter of asymmetry $`\mu `$. At $`\mu =0`$, which corresponds to the symmetrical Gaussian ponderomotive potential shown in Fig. 1a, the distribution is isotropic and purely radial (compare to the result of Ref. ). At the same time, if $`\mu =1.55`$, the number of scattered electrons detected in the $`(𝐄,𝐤)`$ plane ($`\alpha =0`$) is about 30 times higher than that in the $`(𝐇,𝐤)`$ plane ($`\alpha =\pi /2`$). This result is in good quantitative agreement with observations of Malka et al. .

The 30-fold anisotropy of accelerated electrons is clearly explained by asymmetry of the ponderomotive potential (4). The cross term in (4), besides asymmetric corrections to the radial force, gives rise to a tangential force which is responsible for pushing electrons out of the plain perpendicular to the polarization plane. The latter corresponds to minimum of the ponderomotive potential as function of azimuthal angle $`\psi `$, while the perpendicular plane to its maximum. Therefore one could be surprised that the number of electrons scattered at the angle $`\alpha =\pi /2`$ is not equal to zero. Certainly it is explained by complicated structure of the ponderomotive potential (4), namely by the fact that the cross term can change its sign at the periphery of the focus.

The asymmetry of the ponderomotive potential itself is determined by non-zero value of the parameter $`\mu `$ which characterizes relative contributions of E- and H-polarized waves to the resulting field. As far as we know, nobody has never controlled the parameter $`\mu `$ in experiments. The reason is evident. Before the work of Malka et al  there were no experiments where the three-dimensional intensity distribution of a laser pulse influenced physical results. Therefore even qualitative coincidence of our calculations with the results of the experiment  would give a deeper insight into physics of the electron-laser interaction. However, it appeared that the model describes the experiment quantitatively. Of course, the quantitative agreement between the developed theory and the experiment is based on fitting of a single free parameter $`\mu `$ in the problem, which has not been measured experimentally. Therefore, from the standpoint of our model, the experiment of Malka et al.  could be considered as a probe for 3D field distribution in the laser pulse. The correctness of our approach could be verified by an independent experiment for another physical situation performed with the same laser system. Measurements of angular distribution of ATI electrons could serve as a good example of such experiment.

We thank M.V. Fedorov and V.D. Mur for fruitful discussions. This work was supported by the Russian Foundation for Basic Research under projects 00-02-16354 and 00-02-17078.